Development of a Single-Phase, Transient, Subchannel Code, within the MOOSE Multi-Physics Computational Framework

Abstract

Subchannel codes have been widely used for thermal-hydraulics analyses in nuclear reactors. This paper details the development of a novel subchannel code within the Idaho National Laboratory’s (INL) Multi-physics Object Oriented Simulation Environment (MOOSE). MOOSE is a parallel computational framework targeted at the solution of systems of coupled, nonlinear partial differential equations, that often arise in the simulation of nuclear processes. As such, it includes codes/modules able to solve the multiple linear and nonlinear physics that describe a nuclear reactor, under normal operation conditions or accidents. This includes thermal-hydraulics, fuel performance, and neutronics codes, between others. A MOOSE-based subchannel code is a new addition to the fleet of INL-developed codes, based on the MOOSE framework. In this work, we present the derivation of the subchannel equations for a single-phase fluid, we proceed with the description of the algorithm that is used to solve these equations and describe how this algorithm was implemented within MOOSE. We also present how this code can be coupled to the BISON fuel performance code. Next, we verify the friction model and the turbulent mixing model. We calibrate the turbulent modeling parameters for momentum mixing and enthalpy mixing, $C_T,\beta$. We validate the code using experimental results and last demonstrate the coupling capabilities using a simple example.

1. Introduction

Generation-IV reactors (Gen-IV) are a set of nuclear reactor designs currently being researched for commercial applications by the Generation-IV International Forum [1]. A plethora of new designs have been proposed, with a wide variety of configurations and physical attributes, e.g., neutron spectra, coolants, moderators, and relevant feedback mechanisms. The novel designs include: pebble bed reactors (PBRs) both of the high-temperature gas-cooled (HTGR) and fluoride-salt-cooled varieties, [2,3], prismatic HTGRs with hexagonal unit cells and cylindrical fuel slugs [4], micro-reactors (MRs) with annular fuel in hexagonal unit cells [5], hexagonal assembly fast [6], and mixed spectrum reactors [7], hexagonal [8,9], and plate fuel [10], nuclear thermal propulsion systems, and thermal [11], and fast open-pool type liquid fuel molten salt reactors (MSRs) [12,13]. The diversity of the reactor Gen-IV designs, necessitates modeling and simulation (M&S) software that permits flexible multi-physics capabilities.

MOOSE [14] is a Multi-physics Object Oriented Simulation Environment, a parallel computational framework targeted at the solution of coupled, nonlinear partial differential equations (PDEs) that often arise in simulation of nuclear processes. The main advantage of the MOOSE framework is that its a flexible finite element tool that leverages the multi-physics capabilities inherent in MOOSE. Gen-IV reactors present a significant challenge in their analysis due to their complexity, advanced performance and new design features. Developing novel nuclear reactor designs such as the Gen-IV and ensuring their safety under normal operating conditions, operational transients, anticipated operational occurrences, design basis accidents (DBA), etc., required the development of novel computational tools. These codes solve the various physics related to nuclear reactors. Neutronics, fuel performance, and thermal-hydraulics, form the primary set of physics that needs to be resolved.

Subchannel codes are thermal-hydraulic codes that offer an efficient compromise for the simulation of a nuclear reactor core, between CFD and system codes. They use a quasi-3D model formulation and allow for a finer mesh than system codes without the high computational cost of CFD. That is why thermal-hydraulic analysis of a nuclear reactor core is mainly performed using the subchannel type of codes. The safety margins and the operating power limits of the nuclear reactor core are considered as the key parameters for subchannel analysis [15]. The subchannel code is used to calculate: critical heat flux ratio (CHFR), critical power ratio (CPR), fuel center line temperature, fuel pin surface temperature, subchannel maximum temperature, bulk coolant outlet temperature [16] along with the distributions of mass flow, cross flow, pressure, enthalpy, temperature, density and viscosity inside the reactor core. CHFR, CPR and fuel center line temperature are the main parameters that limit the maximum operational power of a nuclear reactor [17]. Moorthi et al. [18], presents a comprehensive list of subchannel analysis codes. Among these codes, most prominent are the COBRA [19] type of codes which are the basis of many of the other codes. Other subchannel codes include START [20], MATRA [21,22], VIPRE [23,24,25], etc.

Since most of the nuclear reactors are designed based on the subchannel analysis codes, it is imperative to develop a subchannel code within MOOSE. A MOOSE-based subchannel code extends the area of impact for the fleet of INL developed codes based on the MOOSE framework. The advantage of the code we developed over other subchannel codes is that it is based on the MOOSE framework which means that the code naturally fits into the environment of all other codes and can be easily coupled to them. This subchannel code has been initially developed for water applications in square lattice arrangements, because of the high availability of water benchmarks. This gives us the opportunity to calibrate and validate the solver for various scenarios. Concurrently, it is also being developed for simulating sodium fast reactors (SFR’s) with liquid sodium coolant in a hexagonal lattice.

2. Methods

2.1. Subchannel Governing Equations

The subchannel thermal-hydraulic analysis is based on the conservation equations of mass, linear momentum and energy on the specified control volumes. We present the control volumes in Figure 1 [19,26].

The control volumes communicate in both axial and lateral directions to produce the three dimensional effect of the reactor core. The subchannel equations are derived when we integrate and average the conservation equations over the subchannel control volume $V_i$ excerpt for the lateral momentum equation which is integrated over control volume $V_{ij}$. For a single-phase fluid we have:

• Conservation of mass:

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

  (1)

  where i is the sub-channel index and j the index of the neighbor subchannel. $\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 sub-channel i in the axial direction. ${\rho }_i$ is the density and $w_{ij}$ [kg/sec] is the diversion crossflow in the lateral direction from sub-channel i to neighboring subchannel j, resulting from local pressure differences between the two subchannels.
• Conservation of linear momentum in the axial direction:

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

  (2)
  In the left hand side we retrieve the change of momentum in the axial direction $\Delta (\frac{{\dot{m_i}}^2}{S_i^z{\rho }_i})$ and the inertia transfer between subchannels due to diversion crossflow ${\sum }_j{w}_{ij}{U}^{*}$. $U^{*}$ is the axial velocity of the donor cell and $-g{\rho }_i{S}_i\Delta Z$ represents gravity force, where g is the acceleration of gravity $S_i$ is the flow area and $\Delta Z$ the axial length of the control volume $V_i$. We have assumed that gravity is the only significant body force in the momentum equation. The donor cell is the cell from which crossflow flows out of 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. Henceforward, donor cell quantities will be denoted with the star (${}^{*}$) symbol. $Frictio{n}_i$ is caused by fluid rod interface and possible local form loss coefficient due to spacers/mixing-vanes. $Dra{g}_{ij}$ is caused by viscous stresses within the fluid at the interface between subchannels. $\Delta P_i$ is the axial pressure drop.
• Conservation of linear momentum in the lateral direction:

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

  (3)
  Here, $S_{ij}=g_{ij}\Delta z$ is the lateral flow surface area between the subchannels i,j. $g_{ij}$ is the gap between subchannels. Lateral pressure gradient ($\Delta P_{ij}/L_{ij}$) across the subchannels and/or 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}$ absorbs all the viscous effects and form losses associated with momentum exchange between the fluid and the wall due to the fluid motion through the gap.
• Conservation of enthalpy:

  $$\frac{\mathrm{d}{〈\rho h〉}_i}{\mathrm{d}t}{V}_i+\Delta \left({\dot{m}}_i{h}_i\right)+\sum _j{w}_{ij}{h}^{*}+h_{ij}^{\prime }=q_i^{\prime }\Delta Z.$$

  (4)
  For a single-phase liquid we can neglect dissipation due to viscous stresses and set the total derivative of pressure (work of pressure) to zero. We can also neglect the volumetric heat source due to moderation, since heat is mainly transferred to the fluid through the fuel rods surface. Here, we have defined $h_{ij}^{\prime }$ as the turbulent enthalpy transfer between sub-channels i,j, $\varphi$ is the dissipation due to viscous stresses and $q_i^{\prime }$ is the average linear power $\left[\frac{\mathrm{kW}}{\mathrm{m}}\right]$ going into control volume of subchannel i from the fuel rods.

2.2. Closure Models

In this section we present the additional equations needed to close the system of the subchannel equations.

• 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}.$$

  (5)

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

  $$Frictio{n}_{ij}=-\frac{1}{2}{g}_{ij}\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 }^{*}}.$$

  (6)

  where $K_{ij}$ is an overall loss coefficient, encompassing lateral concentrated form and friction losses. ${\rho }^{*}$ is the donor cell density.
• Friction factor [22]:

  $$f_w\to \left\{\begin{array}{cc}\frac{1}{64},\hfill & Re<1\hfill \\ \frac{64}{Re},\hfill & 1\le Re<5000\hfill \\ 0.316R{e}^{-0.25},\hfill & 5000\le Re<30000\hfill \\ 0.184R{e}^{-0.20},\hfill & 30000\le Re\hfill \end{array}\right.$$

  (7)
• Turbulent momentum transfer:

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

  (8)

  where $C_T$ is a turbulent modeling parameter.
• Turbulent enthalpy transfer:

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

  (9)
• Turbulent crossflow:

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

  (10)

  where $\beta$ is the turbulent mixing parameter or thermal diffusion coefficient and $\overline{G}$ is the average mass flux of the adjacent subchannels. The $\beta$ term is the tuning parameter for the mixing model. Physically, it is a non-dimensional coefficient that represents the ratio of the lateral mass flux due to mixing to the axial mass flux. In single-phase flow no net 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}$) [26], for enthalpy and momentum, respectively. Simplifying for single-phase, we have made the approximation that the rate of turbulent exchange for energy and momentum are equal:

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

  (11)
• Equation 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)$.h

2.3. Algorithm

For the purposes of implementing a single-phase subchannel code within MOOSE, we developed a hybrid numerical method of solving the subchannel equations presented in the previous section. Hybrid in this context means that the user will have the option of solving portion of the domain at a time by dividing the domain into chunks or blocks as seen in Figure 2. NB is the total number of blocks and NC is the number of axial cells per block. For now, all blocks have the same size (same number of axial cells) but in the future we could allow for blocks of different sizes.

The core of the code is to construct a combined residual function based on the lateral momentum equation, Equation (15). To solve this equation we use a Jacobian Free Newton-Krylov type Method (JFNKM). The workhorse of our code is the non linear equation solvers (SNES) found in the Portable, Extensible Toolkit for Scientific Computation (PETSc) [27]. The SNES class includes methods for solving systems of nonlinear equations of the form: $\mathbf{F}\left(\mathbf{x}\right)=0$. PETSc’s default method for solving the nonlinear equation is Newton’s method. The general form of the n-dimensional Newton’s method becomes:

$${\mathit{x}}_{\mathit{k}+\mathbf{1}}={\mathit{x}}_{\mathit{k}}-\mathcal{J}{\left({\mathit{x}}_{\mathit{k}}\right)}^{-1}\mathit{F}\left({\mathit{x}}_{\mathit{k}}\right),k=0,1,\dots ,$$

(12)

In practice, the Newton iteration (12) is implemented by the following two steps:

$$1. Approximatelysolve\mathcal{J}\left({\mathit{x}}_{\mathit{k}}\right)\Delta {\mathit{x}}_{\mathit{k}}=-\mathit{F}\left({\mathit{x}}_{\mathit{k}}\right)$$

(13)

$$2. Update{\mathit{x}}_{\mathit{k}+\mathbf{1}}={\mathit{x}}_{\mathit{k}}+\Delta {\mathit{x}}_{\mathit{k}}$$

(14)

A JFNK method uses a Krylov method to solve the linear system (13), that appears in each Newton step. Key to a Newton type solver is the forming of the Jacobian matrix $\mathcal{J}$. We usually prefer to approximate the Jacobian. One such approximate comes from “Picard linearization” which writes the nonlinear system as: $\mathit{F}\left(\mathit{x}\right)=\mathit{A}\left(\mathit{x}\right)\mathit{x}-\mathit{b}$, where $\mathit{A}\left(\mathit{x}\right)$ usually contains the lower-derivative parts of the equation. The Picard iterative procedure can be performed by defining $\mathcal{J}\left(\mathit{x}\right)$ to be $\mathit{A}\left(\mathit{x}\right)$. Sometimes this iteration exhibits better global convergence than Newton linearization. In a JFNK method such as the one we use here, we do not calculate the Jacobian, rather we approximate the action of the Jacobian on a vector: $\mathcal{J}\left(\mathit{u}\right)\mathit{\alpha }\approx \frac{\mathit{F}(\mathit{u}+h\ast \mathit{\alpha })-\mathit{F}\left(\mathit{u}\right)}{h}$

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

(15)

The main unknown variable in this non linear residual is the crossflow $w_{ij}$. The combined residual function that we created, calculates the non linear residual $f\left(w_{ij}\right)$ after it updates the other main flow variables, such as mass flow ${\dot{m}}_i$, turbulent crossflow $w_{ij}^{\prime }$, pressure drop $\Delta P_i$ and pressure $P_i$ through Equations (1), (2) and (10), using the current $w_{ij}$ as needed. So every time this function is called by the Newton solver the flow variables get updated. This means that as we solve the lateral momentum equation, at the same time we are solving for all the flow variables. Just as in COBRA [19], in this code as well, the variable $P_i$ is 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. The combined residual vector $f\left(w_{ij,k}\right)$ has a size equal to the number of gaps per level, times the number of axial levels in the block. Each block solution calculates all the unknowns that live in that block at the same time. By defining the number of blocks of the domain equal to one, the combined residual function includes all of the unknown variables in all axial levels. This way all the variables are tightly coupled and solved together at the same time. By defining the number of blocks equal to the number of axial cells, i.e., the size of the block to have one axial level, we retrieve a marching type scheme. By giving a value in between those two extremes, we can have a local tight coupling without paying the price of full scale coupling. Our hybrid scheme is shown in Figure 3.

Once the main flow variables converge in a block, the enthalpy conservation equation, Equation (4) is solved and we retrieve enthalpy $\left(h\right)$ in all the nodes of the block (in the case where no heat is added to the fluid, the enthalpy does not need to be calculated unless we have a non-uniform enthalpy inlet distribution). Using enthalpy, pressure and the equations of state we can calculate the temperature $T_i$ and the fluid properties such as density ${\rho }_i$ and viscosity ${\mu }_i$. After we update the fluid properties, the procedure is repeated until the temperature solution has converged. Once the temperature solution converges we move on to the next block. Once the temperature solution converges in all blocks we check if pressure has converged in all blocks. If not, we repeat the procedure starting again from the first block, until pressure has converged.

3. MOOSE Integration

3.1. MOOSE Solution Method and Structure

MOOSE is a parallel computational framework targeted at the solution of systems of coupled, nonlinear partial differential equations (PDEs), that often arise in the simulation of nuclear systems [14]. In the case of individual physics, a MOOSE application solves for the variables of the system simultaneously using a JFNK scheme. In the multi-physics case, the user can opt to either add physics to the problem formulation and solve for the variables of the system simultaneously, or solve the individual physics separately and couple them via Picard iterations.

The MOOSE framework works on minimizing the residual of a system of PDEs at quadrature points of a finite element mesh using the JFNK method. Forming the Residual begins with writing a weak form of the specific set of PDEs. The weak form is derived from the governing equations by substituting the solution with a finite element expansion and multiplying with test/basis function ${\psi }_i$. Then this weak form is integrated over the whole domain. Using integration by parts and the Gauss theorem produces the boundary conditions and other integrals all of which are called Kernels.

3.2. Subchannel Integration into MOOSE

Developing a subchannel application within the MOOSE computational framework posed a unique challenge. This is because MOOSE is a framework that works on minimizing a residual of partial differential equations at quadrature points of a finite element mesh, which is different to the standard subchannel integration approach. Specifically, we could not use the finite element problem context in MOOSE to solve the subchannel equations but we still had to inherit and make use of the entire finite element mesh, to make coupling to other applications as seamless as it is for other FEM-based solvers in MOOSE. This allows us to keep using the field transfers across different MOOSE applications and it was accomplished using the ExternalProblem class found in MOOSE.

ExternalProblem is an intermediate class that is used as an extension point (inheritance) for externally coupled codes (MOOSE-Wrapped Apps). To achieve the external coupling, this interface overrides the “solve” method along with a method for testing convergence of the last solve. The ExternalProblem object contains the extra interface function necessary for creating and maintaining current solution information on an “external” Mesh object. Specifically, the ExternalProblem creates a final override for the solve () method that triggers three different callbacks: syncSolutions (Direction::TO_EXTERNAL_APP); externalSolve (); syncSolutions (Direction::FROM_EXTERNAL_APP);

Because Subchannel is not based on a FEM-based solver, its field variables are not prime variables but auxiliary variables. Auxiliary variables are used to compute or store intermediate quantities that are not the main variables, the ones being solved for, of the equation system [28]. The prime/main unknowns are often referred to as nonlinear variables. On the other hand, “auxiliary” variables act, with respect to the interaction with other objects, the AuxKernels.

An AuxKernel [29] computes values that are stored in “Auxiliary Variables”. There is an infinite number of ways to generate these values but some common ones are: based on other variables/data in the simulation (i.e., coupled data), read from an external file and interpolated from a Function. The Auxiliary system is extremely flexible on purpose. The data values computed by AuxKernels and stored in Auxiliary Variables (AuxVariables) are often used for visualization (i.e., written to output files) but can also be coupled back into other calculations (including in Kernels) or provided as the input to other systems such as Postprocessors.

In order for MOOSE to understand which variables should be maintained on the ExternalMesh, another overridable method was created that is called during the problem setup to assist the adding the necessary AuxVariables to the simulation.

3.3. MOOSE MultiApp

The MOOSE framework was developed to simplify the creation of fully coupled, nonlinear, multi-physics applications. Here, “fully coupled” means that all of the coupled partial differential equations (PDEs) are solved simultaneously within one Newton-based step. Multiple MOOSE-based applications have been created this way over the years [30]. Each MOOSE-based application is made up of physics “modules” that define the PDEs to be solved. Multiple MOOSE-based applications can be stitched together to create applications that comprise the physics of the whole system [31]. Fully coupled multi-physics is useful for dealing with problems where the physics are strongly interacting, but not all multi-physics problems are (or need to be) fully coupled. Examples include loosely coupled systems with multiple time scales and codes which couple with external software.

3.4. BISON

BISON is a 3D finite element-based nuclear fuel performance code applicable to a variety of fuel forms including light water reactor fuel rods, TRISO particle fuel, and metallic rod and plate fuel. It is a multiphysics fuel analysis tool that solves fully-coupled thermomechanical problems. BISON is also based on the MOOSE framework [32] and implicitly solves coupled thermomechanical equations over the domain of a single fuel rod using the JFNK method [33]. The coupling of Subchannel and BISON, and potentially a newtronics code will provide a multi-physics framework for a full reactor core analysis. The BISON governing equations consist of fully-coupled partial differential equations for energy, species, and momentum conservation. The energy balance is given in terms of the heat conduction Equation (16).

$$\rho C_p\frac{\partial T}{\partial t}+\overrightarrow{\mathbf{\nabla }}·\overrightarrow{q}-E_f\dot{F}=0,$$

(16)

where $T,\rho$ and $C_p$ are the temperature, density and specific heat, respectively, $E_f$ is the energy released in a single fission event, and $\dot{F}$ is the volumetric fission rate. The product gives the volumetric heat source. The heat flux is given in Equation (17)

$$\begin{array}{c}\hfill \overrightarrow{q}=-k\overrightarrow{\mathbf{\nabla }}T,\end{array}$$

(17)

where k denotes the thermal conductivity of the material. $\dot{F}$ can be an input from a separate neutronics calculation (such as RATTLESNAKE or GRIFFIN), or computed based on rod average power and axial profile data.

3.5. Coupling Subchannel with BISON

Coupling between Subchannel and BISON is achieved through the MOOSE MultiApp system. Subchannel is the main app and BISON is the MultiApp. In this coupling Subchannel calculates an average pin surface temperature and the transfer system sends it to a BISON instance as a Dirichlet boundary condition. We have the same number of MultiApps as number of pins we want to simulate. To calculate the pin surface temperature, we have to calculate the Prandtl, Reynolds and Nusselt number, because the heat transfer can be described by the Dittus-Boelter equation, which is:

$$\begin{array}{c}\hfill Nu=0.023R{e}^{0.8}P{r}^{0.4},0.6\le Pr\le 160,\\ \hfill Re>10000.\end{array}$$

(18)

where $Nu$ is the Nusselt number, $Re$ is the Reynolds number and $Pr$ is the Prandtl number. To calculate the Prandtl number, we need to know the following fluid properties:

• The thermal conductivity of the fluid $k\left[\frac{\mathrm{W}}{\mathrm{mK}}\right]$
• The dynamic viscosity of the fluid $\mu \left[\frac{\mathrm{N}\sec }{{\mathrm{m}}^2}\right]$
• The specific heat $C_p\left[\frac{\mathrm{Kj}}{\mathrm{K}\mathrm{g}\mathrm{K}}\right]$

$$\begin{array}{c}\hfill Pr=\frac{\mu C_p}{k}.\end{array}$$

(19)

To calculate the Reynolds number we need to know the following quantities:

• The hydraulic diameter of the subchannel $D_h\left[\mathrm{m}\right]$ which in turn depends on:

  −

  Pitch $P\left[\mathrm{m}\right]$

  −

  Pin Diameter $D\left[\mathrm{m}\right]$

• The dynamic viscosity of the fluid $\mu \left[\frac{\mathrm{N}\sec }{{\mathrm{m}}^2}\right]$
• The fluid density $\rho \left[\frac{\mathrm{Kg}}{{\mathrm{m}}^3}\right]$
• The fluid velocity $V\left[\frac{\mathrm{m}}{\sec }\right]$

$$\begin{array}{c}\hfill Re=\frac{\rho V{D}_h}{\mu }.\end{array}$$

(20)

The convective heat transfer coefficient, $h\left[\frac{\mathrm{KW}}{{\mathrm{m}}^2\mathrm{K}}\right]$, is given directly by the definition of Nusselt number:

$$\begin{array}{c}\hfill Nu=\frac{h{D}_h}{k}\to h=\frac{Nuk}{D_h}\end{array}$$

(21)

Finally, we can calculate the pin surface temperature $T_{pin}\left[\mathrm{K}\right]$, simply using the Newton’s Law of Cooling:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill q^{″}=h(T_{pin}-T_{bulk})\Rightarrow q^{\prime }=\pi Dh(T_{pin}-T_{bulk})\Rightarrow \\ \hfill T_{pin}=\frac{q^{\prime }}{\pi Dh}+T_{bulk}.\end{array}\end{array}$$

(22)

Because each pin is in contact with four subchannels the above calculation is repeated four times and the final result is taken to be the arithmetic mean. The averaged value of $T_{pin}$ is given as a Dirichlet boundary condition to BISON.

Each pin solution requires one BISON instance. Each BISON instance is considered to be one MultiApp. BISON uses a 2D azimuthally symmetric (RZ) model and calculates heat conduction. After BISON solves for temperature it calculates the linear heat rate at the surface using the following equation:

$$\begin{array}{c}\hfill q^{\prime }=\pi D{q}^{″}=-\pi Dk\frac{\mathrm{d}T}{\mathrm{d}x}.\end{array}$$

(23)

The kernel that calculates this quantity is derived from MOOSE DiffusionFluxAux Kernel. The DiffusionFluxAux AuxKernel is used to compute the flux vector for diffusion problems. The flux is computed as $J=-D\frac{\partial C}{\partial X}$, where J is the diffusion flux vector, D is the diffusivity or diffusion coefficient, C is the concentration variable, and X is the coordinate. After the BISON calculation is complete, the linear heat flux is passed on to subchannel and the calculations are repeated via a Picard iteration. The data transfer between the two codes that allows the coupling is shown in Figure 4.

It should be noted that because the field variables of the Subchannel application are auxiliary variables and not primary variables, the calculation of a non-linear residual for Picard iterations is not possible for now. We can nevertheless by-pass this hurdle, by manually defining the number of Picard iterations. In a way we high-jack the Picard iteration interface in MOOSE to achieve the coupling. For the simple demonstration example presented in this paper one data exchange is enough for convergence.

Because the Subchannel mesh is 3-D and exists on dimensions, while the BISON model is 2-D and exists on the x-y plane, the meshes don’t align and special care had to be taken in order to accurately transfer the coupled quantities. This was accomplished by using the MOOSE NearestPointLayeredAverage, LayeredSideAverage and MultiAppUserObjectTransfer2 objects. The first object computes averages of a variable storing partial sums for the specified number of intervals in a direction (x, y, z). Given a list of points this object computes the layered average closest to each one of those points. Similarly the second object computes side averages of a variable. MultiAppUserObjectTransfer2 inherits from MultiAppUserObjectTransfer2 and takes care of transferring variables on possibly different meshes while conserving a user defined property of each variable. This object swaps the y- and z-coordinate, so we can map the subchannel solution on a RZ-mesh of a heating pin.

4. Results

4.1. Verification

In this section we present the verification of the friction model and the verification of the turbulent mixing model. The problems discussed are chosen such that they isolate the specific physics in question and are the same used in the CTF Void drift Validation study [34]. The problem geometry is presented in Figure 5.

4.1.1. Two-Channel Friction Model Verification

In this case, we want to study a problem where the effects of friction are clearly discernible. We know that crossflow in the single-phase case is driven by a lateral pressure gradient and turbulence. By deactivating turbulence in our test model case, crossflow can only be the result of lateral pressure imbalance; which for a model with no form losses (spacer grids), can only be driven by unequal frictional losses. Friction loss depends on the hydraulic diameter, so it makes sense to devise a two-channel problem, with channels that have an unequal flow area.

Channel-2 has double the hydraulic diameter of Channel-1. The length of the model is set to 10 m to allow the flow to fully develop. A sketch of the two-channel problem is shown in Figure 5a. The different frictional pressure drops create a lateral pressure gradient that drives flow from the small channel to the large channel. Moving up the channels, velocity decreases in Channel-1, which decreases the frictional pressure drop. At the same time, velocity grows larger in Channel-2, which increases frictional pressure drop. This continues until the frictional pressure drop is the same in both channels, at which point crossflow between the two channels ceases. At this point, the channels are in mechanical equilibrium. An analytical expression is derived for the flow distribution at the point of mechanical equilibrium.

$${\dot{m}}_{in}={\dot{m}}_2\left(1+{\left(\frac{D_{h,2}}{D_{h,1}}\right)}^{\frac{C_2-1}{C_2+2}}\frac{A_1}{A_2}\right).$$

(24)

where $m_{in}$ is the mass flow at the inlet $m_1,m_2$ are the mass flows of the two channels at equilibrium, $D_{h,2},D_{h,1}$ the hydraulic diameters, $A_1,A_2$ the flow surface areas and $C_2$ is a friction modeling parameter when we define the friction factor: $f=C_{1}R{e}^{C_2}$ as a function of Reynolds number. We can compare the analytical prediction with the code results in Figure 6a. The code converges to the analytical prediction.

4.1.2. Turbulent Mixing Verification

In this case, we want to study a problem where the effects of turbulence are clearly discernible. Turbulence causes both momentum and enthalpy mixing through the terms:

$$\begin{array}{c}\hfill Dra{g}_{ij}=-C_T\sum _j{w}_{ij}^{\prime }\Delta U_{ij},\end{array}$$

(25)

$$\begin{array}{c}\hfill h_{ij}^{\prime }=\sum _j{w}_{ij}^{\prime }\Delta h_{ij}.\end{array}$$

(26)

Because the models for turbulent mixing are gradient-driven based on $\Delta h,\Delta U$, in order to observe the effects of turbulence, it is necessary to make a gradient in either enthalpy or velocity. It is easier, to focus on the energy equation and deactivate the density calculation. The problem geometry consists of two identical channels connected by a gap and is seen in Figure 5b. To test the turbulent mixing model, the temperature of one channel is raised by 10 degrees Celsius. Turbulent enthalpy mixing will transfer heat from the one channel to the other. The solution given by the code is then compared to the analytical solution.

$$\begin{array}{c}\hfill h_1=\frac{(h_{1,in}+h_{2,in})}{2}-\frac{1}{2}(h_{2,in}-h_{1,in})exp(-\frac{2\frac{\mathrm{d}{w}_{12}^{\prime }}{\mathrm{d}z}}{\dot{m}}z)\end{array}$$

(27)

$$\begin{array}{c}\hfill h_2=\frac{(h_{1,in}+h_{2,in})}{2}+\frac{1}{2}(h_{2,in}-h_{1,in})exp(-\frac{2\frac{\mathrm{d}{w}_{12}^{\prime }}{\mathrm{d}z}}{\dot{m}}z)\end{array}$$

(28)

where $h_{1,in},h_{2,in}$ are the enthalpies at the inlet of the two channels and $\frac{\mathrm{d}{w}_{12}^{\prime }}{\mathrm{d}z}$ is the turbulent crossflow between the two channels per unit length. We see that the code prediction and the analytical solution are in good agreement in Figure 6b. The maximum absolute error is calculated to be $Erro{r}_{max}=0.2450969285540623\mathrm{kj}/\mathrm{kg}$ and the maximum relative error $Erro{r}_{relative,max}=0.0002814742382126506$.

4.2. Calibration

In the governing equations, Equations (8) and (10), we have introduced the turbulent modeling parameter $C_T$ and the turbulent diffusion coefficient $\beta$. In this chapter we calibrate these parameters using appropriate experimental data.

4.2.1. Turbulent Diffusion Coefficient $\beta$

The $2\times 3$ facility is an air-water facility that was operated by Kumamoto university. The purpose was to quantify the effects of mixing and void drift [35]. In these experiments, the turbulent mixing rates and the fluctuations of static pressure difference between subchannels were measured. The author derived a way to use the die concentration measurements, in order to calculate the turbulent mixing rates ($w_{ij}^{\prime }$) between subchannels [36].

The experimental results are presented through the dimensionless number $\frac{W_{ij}}{dz\mu }$ versus $R{e}_{ij}$, where $R{e}_{ij}=\frac{\rho u_{ij}{D}_{h,ij}}{\mu }$ and $u_{ij}=\frac{u_i{A}_i+u_j{A}_j}{A_i+A_j}$. Here, u and A are the mean velocity and the cross-sectional areas, the subscription $i,j$ indicate the subchannel index, $\rho ,\mu$ are the density and viscosity of the fluid. Lastly, $D_{h,ij}$ is the hydraulic diameter of an imaginary channel with subchannels: $i,j$. By changing the parameter $\beta$ and deactivating momentum diffusion we see that the experimental results are bounded by values of $\beta =0.004,\beta =0.009$ as seen in Figure 7a,b.

Plotting the $L_2-$ norm of the difference between code prediction and experimental results, as well as the sum of the absolute difference between code prediction and experimental results, versus the parameter $\beta$, we get Figure 7d. A mixing coefficient of $0.006$ leads to the least amount of error. Figure 7c presents the code prediction vs. the experimental results for the value of the optimum parameter $\beta =0.006$. It is important to note, that the mixing coefficient is simply a tuning parameter that depends on the actual geometry of the specific facility being modeled. The Kumamoto facility is not an exact representation of a PWR, nevertheless this study is useful for demonstrating that the code is capable of predicting the correct mixing rate, if $\beta$ is calibrated accordingly [34].

4.2.2. Turbulent Modeling Parameter $C_T$

After calibrating the turbulent diffusion coefficient $\beta$ we turned our attention to the turbulent modeling parameter $C_T$. This is a tuning parameter that informs on how much momentum is transferred/diffused between subchannels, due to turbulence. The CNEN $4\times 4$ test [37] performed at Studsvik laboratory for studying the flow mixing effect between adjacent subchannels was chosen to tune this parameter. This experiment consists in velocity and temperature measurements taken at the outlet of a 16-rod assembly test section. Analysis of the velocity distribution at the exit of the assembly can be used to calibrate the turbulent parameter $C_T$.

For mass-flux = 5425.20 $\left[\frac{\mathrm{kg}}{{\mathrm{m}}^2\mathrm{s}}\right]$ we plot the code predictions along with the experimental results in Figure 8a,b. A value of $C_T=0$ deactivates the turbulent momentum diffusion. For this case, flow redistribution within the assembly will be the result of unequal friction losses, due to the unequal flow areas of the subchannels. Flow tends to get transferred from the smaller channels to the larger ones. It is evident, judging by Figure 8a, that the friction model alone, is not enough, for the code to predict the correct flow redistribution. We can see in Figure 8a that for this value of the modeling parameter, velocity redistribution is overestimated.

This is not unexpected, because turbulent momentum diffusion counteracts the effect of friction. We know that turbulence tends to flatten velocity profiles and transfer momentum towards the wall, while we have already seen that friction transfers momentum towards the center. A value of $C_T=3$ enables the turbulent diffusion model. We can see in Figure 8b that for this value of the modeling parameter, velocity redistribution is underestimated. We can therefore ascertain that the optimum value for $C_T$ will be between 0–3.

Another point we should note is the effect of the spacer grid in the middle of the assembly. The test section was fitted with a grid of low form coefficient (k = 0.3), intended to maintain the alignment of the rods [37]. The effect of the spacer grid is felt through the pressure drop calculation. The form coefficient raises the pressure drop in middle of the assembly which in turn raises the local pressure. Because the pressure drop is a function of velocity squared, pressure drop is higher for the channels with higher velocity, meaning the larger channels vs. the smaller ones. That means that locally flow will be transferred from the larger channel to the smaller one and that is exactly what we see in the figures. We plot the absolute average error in the exit velocity prediction vs. $C_T$ in Figure 8d. We deduce that the optimum value for the modeling parameter is $C_T=2.6$. The code calculation using the optimum value is shown in Figure 8c.

4.3. Validation

4.3.1. Validation (PSBT $5\times 5$)

The PSBT $5\times 5$ benchmark is an international benchmark developed by the Organisation for Economic Co-operation and Development (OECD), the Nuclear Regulatory Commission (NRC), and the Nuclear Power Engineering Center (NUPEC) [38]. In this work, we utilize the steady-state mixing test, detailed in volume III of the PSBT benchmark [39]. The purpose of this test is to validate/tune the mixing models of the participating codes. The participating codes predict the fluid temperature distribution at the exit of the heated section of a rod bundle assembly and compare it with the experimental values provided by the benchmark. The rod bundle geometry description is presented in

Table 1. $5\times 5$ PSBT rod bundle specifications.

Item	Value
Rods array	$5\times 5$
Number of heated rods	25
Heated rod outer diameter (mm)	$9.50$
Heated rod pitch (mm)	$12.60$
Axial heated length (mm)	3658
Flow channel inner width (mm)	$64.9$
Axial power shape	Uniform
Number of mixing vaned (MV) spacers	7
Number of non mixing vaned (NMV) spacers	2
Number of simple spacers (SS)	8
MV spacer location (mm)	457, 914, 1372, 1829, 2286, 2743, 3200
NMV spacer location (mm)	0, 3658
Simple spacer location (mm)	229, 686, 1143, 1600, 2057, 2515, 2972, 3429

(reprinted from [39]).

The rod bundle has a lateral power profile in which the right side of the assembly is under-heated. The radial power profile is shown in Figure 9. The rods on the left side transfer 100% of available rod power to the fluid, while the rods on the right transfer 25%. This causes an uneven temperature distribution at the exit of the assembly. The linear heat rate of the rods is uniform.

Note that the turbulent mixing parameter used for all the cases was: $\beta =0.08$, which differs significantly from the value of $\beta =0.006$ that we estimated to be the optimum value in Section 4.2.1. Looking at the temperature profile in Figure 10b we see that $\beta =0.006$ severely underestimates the enthalpy diffusion in comparison to the experimental results. So clearly $\beta$ had to to be adjusted to a much higher value. That value was adopted from [40] and produces code predictions that are in a much better agreement with the experimental results.

The reason behind the need to adjust $\beta$ to a much higher value, has to do with the geometry of the PSBT facility. Note that there is a preferential mixing direction of the experimental results in the diagonal direction (towards one corner of the assembly), exhibited in Figure 10b, while the code results for both values of $\beta$ are symmetric on axis x as expected. The experimental results exhibit a non symmetric distribution that we cannot capture with a constant value of $\beta$. There is a temperature gradient towards the corner due to an additional mixing effect, which may reduce the exit temperature differences between subchannels and finally lead to an increase of the value of the optimum $\beta$ [40].

The reason for the additional mixing effect is thought to be the special mixing vanes that the NUPEC facility incorporates in its design. Specifically, the temperature gradient appearing in the experimental data was attributed to the thermal mixing in the diagonal direction of the test bundle, which may be caused by the alignment of mixing vanes mounted in the spacer grids [40]. This is the physical reason behind the need to use an increased value of $\beta =0.08$.

This illustrates the fact that modeling parameters such as $\beta$ should ideally be calibrated for specific geometries and in no way can they be applied generally without proper justification. Nevertheless our results show that, provided we choose the optimum parameters adjusted for the specific geometries, we can accurately predict the exit temperatures. Figure 11 presents the cumulative mean absolute error in the exit temperature in comparison with other subchannel codes [39]. We note that for the temperature mixing test, our code performs adequately in comparison to the other codes. Our mean absolute error is calculated to be $Error=2.176$ which places us in 5th place out of the nine codes.

4.3.2. Validation (PNL $2\times 6$)

The PNNL $2\times 6$ benchmark [41] was performed at Pacific Northwest National Laboratory in 1977 for investigating the buoyancy effect in the mixed (combined free and forced) convection regime for specific flow coast-down transients. The objective of the study was to develop an understanding of the thermal hydraulic phenomena at low flows in a rod bundle subjected to non-uniform lateral power distributions. The study was performed using a non-uniformly electrically heated $2\times 6$ rod bundle contained in a flow housing.

Local fluid velocity and temperature measurements were obtained. This benchmark allows us to validate the subchannel code in natural/forced convection conditions not tested until this point. For each of the test conditions of this study, fluid axial velocity and temperature (both local and bulk average at the inlet and outlet) were measured within the heated length of the rod bundle. For the cases presented here, velocity measurements were made along the X axis at Y = 0.0 which is the centerline of the central subchannels. Velocity data were recorded at selected points along this line. Cases 5, 9 and 13 where chosen to test the code. Case 5 was chosen to demonstrate the transient capability of the code, while cases 9, 13 where chosen to demonstrate the code behaviour in a buoyancy driven scenario. The measuring point in case 5 was the center of the central subchannel at an axial location equal to 0.125× the active length of the rod bundle. The nominal run conditions for these cases are shown in

Table 2. Nominal run conditions.

Case Number	Initial Flow [GPM]	Final Flow [GPM]	Transient Time [s]	Power Gradient	${\mathit{Q}}_{\mathit{H}}$ $\left[\frac{\mathbf{kW}}{\mathbf{rod}}\right]$	${\mathit{Q}}_{\mathit{L}}$ $\left[\frac{\mathbf{kW}}{\mathbf{rod}}\right]$	Re Initial/Final
5	3.08	1.08	150	0:0	0.0	0.0	1200/420
9	3.08	3.08	S.S	1:0	0.91	0.0	1290/1290
13	3.08	3.08	S.S	2:1	0.91	0.455	1340/1340

.

The code results vs. the experimental results are shown in Figure 12a–c. The code predicted average values in all three cases are lower than the measured centerline experimental values. That is because the experimental results are point-wise instantaneous velocity measurements and the code results are the surface averaged velocities in each subchannel. The experimental maximum values measured at the subchannel centers are less than the analytically predicted value of $U_{max}=2\times U_{average}$ for laminar flow inside a circular pipe. This happens because even though subchannel codes treat subchannels as equivalent to pipes, we also include a version of turbulence modeling and crossflow which tend to flatten the velocity profiles. In turbulent flows momentum is transferred towards the wall region, hence the maximum velocity in turbulent profile is less than that in the laminar case.

Furthermore we notice that in the case where the power ratio is reduced (case 13 vs. case 9), the velocity profile gets more flat. and the code prediction does indeed follow that trend. This happens because buoyancy effects are more pronounced when we have more extreme gradients in heat rate. In all cases, our results follow the trend of the experimental results. Our results are therefore consistent with what we physically expect to see and what is experimentally measured. Lastly, it should be noted that the initial bump in the code prediction in Figure 12c relates to the time required for the boundary condition information to reach the point of the measurement downstream.

In both validation problems and in all other instances of running the code we have observed rapid convergence of the non-linear iterations.

4.4. 4 × 4 Sub-Assembly Coupling Example

In this section we present an example of Subchannel/BISON coupling, with a geometry typical for an LWR. For this purpose we are going to use the SmearedPelletMeshGenerator class provided by BISON. This class creates a 2D-RZ axisymmetric mesh of a fuel rod. More specifically, it generates a smeared pellet mesh for use in simulating fuel rods with pellet fuel. In a smeared mesh, the fuel stack is treated as a single rectangle of fuel (i.e., dishes and chamfers are not included).

This geometry includes an unmeshed gap between the fuel and the cladding. To model the heat transfer through the gap we use the GapHeatTransferLWR object [42]. The gap conductivity in this case is computed based on the gases in the gap. For the example presented here, conduction through radiation and through contact is neglected.In BISON the heat transfer coefficient in the gap through the gases $h_g$ is modelled using the form proposed by Ross and Stoute (1962) [43].

$$\begin{array}{c}\hfill h_g=\frac{k_g\left(T_g\right)}{d_g+C_r(r_1+r_2)+g_1+g_2},\end{array}$$

(29)

where $k_g$ is the conductivity of the gas in the gap, $d_g$ is the gap width, $C_r$ is a roughness coefficient, $r_1$ and $r_2$ are roughnesses of the surfaces, and $g_1$ and $g_2$ are jump distances, which become important for small gap widths and low gas pressures. The jump distances provide a reduction in gap conductance when the mean free path of the gas molecules is significant in comparison to the gap width, and the continuum approximation is no longer valid. The gas temperature $T_g$ is the average of the two surfaces.

The fuel and the cladding are modeled using the BISON UO2Thermal and ZryThermal objects, respectively [44,45]. The geometry of the demonstration problem is based on typical PWR specifications and the layout of the rod bundle is shown in Figure 13. There are 16 non-uniform heated rods and 25 subchannels. All general specifications can be found in

Table 3. $4\times 4$ PWR sub-assembly specifications.

Type	Value	Units
Reactor	PWR
Layout	$4\times 4$
Fuel	UO2
Cladding	Zyrcalloy-4
Fuel pellet radius	$0.409$	cm
Gap	$0.009$	cm
Clad inner radius	$0.418$	cm
Clad thickness	$0.057$	cm
Clad outer radius	$0.475$	cm
Pin pitch	$1.276$	cm
Pin length	360	cm
Gap bottom	$0.127$	cm
Gap top	$0.127$	cm
Clad thickness top	$0.224$	cm
Clad thickness bottom	$0.224$	cm
Average linear heat rate per pin	$18.5$	kW/m
Bundle power	$1.0656$	MW
Coolant inlet temperature	$560.15$	K
Outlet pressure	$15.5$	MPa
Mass flux at inlet	$3035.27$	kg/sec-m2

adjusted from [46]. We have 36 axial cells in subchannel, 1000 layers in the averaged objects and 100 axial cells in the cladding and fuel pellet in the BISON model.

The user has the option to provide a power profile both in the lateral direction (power per pin) as well as in the axial direction (volumetric heat rate in fuel region of the pin (BISON)/linear heat rate in the surface of the pin (subchannel)). The lateral power distribution of the pins is shown the in

Table 4. Lateral power profile for the $4\times 4$ PWR sub-assembly.

1.0	0.9	0.8	0.7
1.1	1.0	0.9	0.8
1.2	1.1	1.0	0.9
1.3	1.2	1.1	1.0

. The coldest subchannel is 5 and the coldest pin is 4. The hottest subchannel is 21 and the hottest pin is 13. We use BISON to model the two corner pins (4, 13). In a standalone application the power distribution data is provided directly to the subchannel code by defining a linear heat rate per pin. In this coupling scenario, the power distribution data is provided to BISON for the two pins that are simulated by BISON and to subchannel for the rest of the pins.

In addition, the energy deposition from the fuel region of the pin is not uniform in the axial direction. As a simple approximation, the axial power distribution follows the cosine function with the peak energy deposition occurring halfway along the active fuel length. This also has a direct effect on the temperature profiles inside the cladding material. The linear heat rate of the pin surface is given by.

$$\begin{array}{c}\hfill q_{pin}^{\prime }\left(z\right)=q_{pin}^{\overline{\prime }}\frac{\pi }{2}sin\left(\pi \frac{z}{L}\right).\end{array}$$

(30)

where $q_{pin}^{\overline{\prime }}$ is the average linear heat rate per pin and L is the active pin length. The pin power profile for the hot/cold Pin is shown in Figure 14.

We run the subchannel code and plot the subchannel and pin surface temperature fields in Figure 15a,b. We plot cross-sections of the temperature field in Figure 16. We plot the axial temperature profiles for the hot and cold corner (pin and subchannel) in Figure 17a. We plot the center-line axial temperature calculated in BISON for the hot and cold pin in Figure 17b and the radial temperature profile at $z=2.4$ m in Figure 17c. The axial temperature profile is affected by the pin power profile which has a sinusoidal shape. The subchannel code can calculate the flow variables such as bulk fluid temperature and use that to compute pin surface temperature. The fuel performance code BISON solves heat conduction within the pin and calculates the temperature distribution. The purpose of this demonstration is to showcase the capabilities of the subchannel code within a multi-physics context in MOOSE. No additional wrapper was needed to be developed to accomplish this coupling.

5. Conclusions and Future Work

5.1. Summary of This Effort

Gen-IV reactors present a significant challenge in their analysis due to their complexity, advanced performance and new design features. The diversity of the reactor Gen-IV designs, proposed by the industry in recent years, necessitates M&S software that permits flexible multi-physics capabilities. Multi-physics tools such as MOOSE, were developed for the solution of systems of coupled PDEs, that often arise in the simulation of nuclear systems.

Subchannel codes offer an efficient compromise for the simulation of a nuclear reactor core, between CFD and system codes. They use a quasi-3D model formulation and allow for a finer mesh than system codes without the high computational cost of CFD. That is why thermal-hydraulic analysis of a nuclear reactor core is mainly performed using the subchannel type of codes to estimate the various thermal-hydraulic safety margins and the various quantities of interest.

Since most of the nuclear reactors are designed based on the subchannel analysis codes, it is imperative to develop a subchannel code within MOOSE. A MOOSE-based subchannel code is a new addition the fleet of INL developed codes based on the MOOSE framework. The advantage over other subchannel codes is that it fits naturally into the MOOSE family of codes and can be easily coupled to them. This work presents the development of a single-phase, transient subchannel code, within the MOOSE computational framework, including validation and verification and demonstration of multi-physics coupling.

5.2. Future Work

We conclude this presentation by reflecting on possible ways to improve and expand on this effort. We can summarize future work as follows:

5.2.1. Re-Circulation

Subchannel cannot calculate negative mass flows in the axial direction. For this reason we cannot capture recirculation phenomena in cases where we have the formation of eddies due to high forced convection and low flow conditions, such as those seen in [41], or due to blockages in the flow path such as those seen in [47]. Possible accident scenarios that could produce these types of flow effects could be the loss of power in a large part of the core, or pin swelling and pin rupture in the middle of the core. In order to deal with these kinds of scenarios we would have to edit the subchannel equations to allow for back flow, make use of a staggered velocity-pressure grid and assign fluid properties appropriately on the correct node, according to the sign of mass flow.

5.2.2. Two-Phase Flow

Subchannel is a single-phase code. It cannot calculate fluid quality hence it is not applicable to BWR’s or in cases where we have a second phase in part of the coolant. A two phase capability is important in order to predict CHF or DNB. These are parameters of great importance, which constrain the thermal power capability of an LWR. In order to add a two phase capability we would need add equations so that we could translate the mixture variables the subchannel code solves for, into the two separate phase variables (gas/liquid), that describe the two-phase conditions:

• One equation that relates the mass flow of the mixture with the mass flow of the phases and the volume fraction of the phases.
• One equation that relates the flow quality with the mass flow of the mixture and the mass flow of the phases.
• One equation that relates the density of the mixture with the density of the phases.
• One equation that relates the velocity of the mixture with the velocity of the phases.
• Equations of state for the two fluids.

The OECD/NRC PSBT benchmark [38] consists of two separate phases: a void distribution benchmark and a DNB benchmark, that can be used to validate a two-phase subchannel code for based on a PWR. Similarly the OECD/NRC BFBT benchmark [48] consists of two phases as well: Phase I—Void Distribution Benchmark, Phase II—Critical Power Benchmark, that can be used to validate a two-phase subchannel code for based on a BWR.

5.2.3. Hexagonal Lattice

Subchannel code is based on a square lattice geometry. Gen-IV SFRs have the fuel placed in hexagonal assemblies. Efforts are currently been made jointly by both INL and Argonne, to expand the code capabilities in order to model SFR’s. For this purpose the mapping of the has been re-figured to accommodate the hexagonal lattice. The fluid properties of the coolant changed to that of liquid sodium. The effect of the wire wrap needs added, both in the axial pressure drop calculation and the mixing model. Finally heat conduction within the fluid is taken into consideration and included in the energy conservation equation. To confirm the prediction capability of the SFR version of the code, calculations and comparisons with available experimental data of 19-, 37-pin and 127- pin sub-assemblies can be performed [49]. The validation of the code can also make use with Oak Ridge National Laboratory (ORNL) 19 pin tests [50]. The temperature profiles at the end of the heated length for low and high power cases can be compared between experimental results and other codes [51].

5.2.4. Expand Multi-Physics Scope

The subchannel code has been coupled to BISON which is a fuel performance code. In this coupling scheme, subchannel calculates the pin surface temperature and sends that to BISON. BISON calculates heat conduction within the pin and gives back linear heat rate to subchannel. A more complete model would have to include neutronics calculations. This new effort would be similar to what the Virtual Environment for Reactor Applications (VERA) [52,53] is designed to do.

VERA has been developed by the Consortium for Advanced Simulation of Light Water Reactors over the past 12 years to address difficult problems facing commercial light water reactors. At the core of VERA is the VERA core simulator (VERA-CS), composed of the MPACT neutron transport code [54,55] the CTF subchannel thermal-hydraulics (TH) code [56] and the ORIGEN depletion code from the SCALE package [57]. In addition to the core simulator, other codes are available for more specific analyses, such as the BISON fuel performance code [32,58] and the MAMBA coolant chemistry code [59]. In the most general model we would have three main codes coupled together. Neutronics, fuel performance and a thermal-hydraulics code:

• NEUTRONICS. A neutronics code would generate cross sections, calculate the scalar flux and compute the heat deposition in the fuel and coolant/moderator. The data transfer to the fuel performance code would include fission rate density, power density, and fast neutron flux.
• FUEL PERFORMANCE. A fuel performance code would use the flux/power information to calculate burnup in the fuel along with heat conduction and nuclide inventory. The data transfer to the neutronics code would include temperature, to be used for updating the cross sections. The local temperature and density conditions are used to update the macroscopic cross sections in both the pins and the moderator/coolant.
• THERMAL-HYDRAULICS. A subchannel code would calculate subchannel bulk temperature. It would give that information to the fuel performance code and the neutronics code, while it would receive axial heat rate from the fuel performance code and a volumetric heat rate from the neutronics code.

The objectives VERA was to designed to accomplish will be achieved entirely within the MOOSE framework. With GRIFFIN being the neutronics code, BISON being the fuel performance code and Subchannel being the thermal-hydraulics code. This way there would be no other need for external codes/wrappers to deal with the coupling and data transfers, as for example is the case with the Lightweight Integrating Multi-physics Environment (LIME) problem manager [52,60]. LIME [61,62] was used as the initial basis for non-invasive (e.g., black-box) coupling of capabilities in VERA. It was specifically designed for integration of existing components that solve different phenomena of a multi-physics problem. The next logical step would be a full core analysis using the framework described above. The CASL challenge problems [63] or the NuScale core design [64] would be ideal cases to test the codes.