# Prediction of Flow and Temperature Distributions in a High Flux Research Reactor Using the Porous Media Approach

Authors: Shanfang Huang; Daxin Gong; Chao Li; Xiaoyu Guo; Guanbo Wang; Huaqiang Yin; Kan Wang
Journal: Science and Technology of Nuclear Installations
 (2017)

---

## Abstract

High thermal neutron fluxes are needed in some research reactors and for irradiation tests of materials. A High Flux Research Reactor (HFRR) with an inverse flux trap-converter target structure is being developed by the Reactor Engineering Analysis Lab (REAL) at Tsinghua University. This paper studies the safety of the HFRR core by full core flow and temperature calculations using the porous media approach. The thermal nonequilibrium model is used in the porous media energy equation to calculate coolant and fuel assembly temperatures separately. The calculation results show that the coolant temperature keeps increasing along the flow direction, while the fuel temperature increases first and decreases afterwards. As long as the inlet coolant mass flow rate is greater than 450 kg/s, the peak cladding temperatures in the fuel assemblies are lower than the local saturation temperatures and no boiling exists. The flow distribution in the core is homogeneous with a small flow rate variation less than 5% for different assemblies. A large recirculation zone is observed in the outlet region. Moreover, the porous media model is compared with the exact model and found to be much more efficient than a detailed simulation of all the core components.


---

## Body

## 1. Introduction

There are over 240 research reactors in operation in the world [1]. These research reactors have been constructed to generate neutrons for different purposes, such as irradiation testing of materials and the production of isotopes [2]. The irradiation effect is an important indicator for new materials. Current mainstream research reactors can generate high thermal neutron fluxes for irradiation testing of materials. However, these reactors cannot easily produce high thermal neutron fluxes that also have high energy levels. An inverse flux trap-converter target structure can convert high thermal neutron fluxes into high energy fusion neutrons (about 14 MeV) which are needed in some irradiation tests. For example, the China Mianyang Research Reactor (CMRR) has an inverse flux trap-converter target structure so it can also produce high energy neutrons. These research reactors are characterized by large leakages due to their compact reactor design [3]. The Reactor Engineering Analysis Lab (REAL) at Tsinghua University has designed a research reactor with an inverse flux trap-converter target structure modeled after JRR-3M [4, 5] and CARR (China Advanced Research Reactor) [6, 7] to further investigate these research reactors.To generate high thermal neutron fluxes, the reactor core must be compact with a high core energy density [2]. The coolant channels in a plate-type fuel assembly are quite narrow, so a high coolant flow rate is needed to effectively remove the heat. With these conditions, the fuel cladding surface experiences high heat fluxes which combined with the high coolant flow rate result in a challenge for the thermal-hydraulic design of the HFRR.For safety concerns, the fuel temperature should be lower than 400°C (673 K) and the cladding surface temperature must not exceed 220°C (493 K). No nucleate boiling is allowed in the assemblies and flow instabilities should be avoided [8]. Gong et al. [9] already developed an exact model for these criteria to numerically investigate the heat transfer characteristics in the HFRR fuel assemblies at inlet flow rates of 4.5 m/s to 7.5 m/s. Their results showed that the peak temperatures in both the fuel and the cladding are lower than the design temperatures. The cladding temperature was always lower than the local saturation temperature, so no bubble generation occurred in the fuel assemblies. However, these calculations assumed that the coolant velocities were equal in each of the flow channels even though the flow rates are known to differ. When the differences are large, some fuel assemblies will not be effectively cooled. To ensure safe operation of the reactor, the flow distribution and cooling conditions in each fuel assembly must be accurately known. Thus, full core flow and temperature analyses are needed.Existing CFD models have been well developed and tested to predict the temperature and flow distributions in reactor cores [10–12]. The model in Gong et al. [9] simulated in detail the temperatures in the fuel, the cladding, and the coolant in the assembly channel. However, modeling of the flow and temperature fields throughout the full core with an exact model would require enormous computer resources and processing times. Therefore, a simplified approach is needed to model the entire core while avoiding these limitations.Porous media models are approximate methods that can be used to simplify core simulations. The porous media model adds a source term into the momentum equation to simulate the solid resistance to the fluid flow in the calculation region [13]. The porous media model has been successfully used to simulate several types of reactors to investigate various practical problems, such as the thermal-hydraulic analysis of CARR [14], the modeling of the Missouri S&T Reactor’s natural convection [15], the simulations of the fluid flow and heat transfer in the Advanced CANDU Reactor [16], and the CFD analyses of the MYRRHA primary cooling system [17].The purpose of this research is to develop a full core CFD model using the porous media model for safety analyses of the HFRR. The flow and temperature distributions are predicted throughout the entire core for various inlet flow rates based on neutronics analyses of the HFRR design [2].


## 2. Model Description

### 2.1. HFRR Overview

As Gong et al. [9] described, HFRR is a pool-type research reactor which is moderated and cooled by light water and reflected by beryllium and heavy water, with a maximum thermal power of 20 MW and low-enriched uranium (LEU) plate-type fuel elements. Figure 1 shows the fuel element configuration and the schematic of the follower fuel assembly. The core contains 18 standard fuel elements in yellow, 6 follower fuel elements in grey, and 4 irradiation tubes in blue. The fuel meat is U3Si2-Al and the fuel cladding is Al (6061). Gong et al. [9] showed more details about the structure of the standard and follower fuel elements in their paper. The core center has an inverse flux trap made of Be which has a peak thermal neutron flux of 4.04 × 1014 n/cm−2s. Figure 2 shows the 3D schematic of the reactor core as well as the sketch of the HFRR layout. Besides the fuel elements, the core also contains a heavy water tank installed around the core and a lower plenum. During normal operation, the coolant flows down the core by forced convection. The inlet coolant temperature is 35°C and the working pressure is 0.152 MPa.Figure 1
HFRR core configuration and schematic of follower fuel assembly [9]: (a) fuel element configuration; (b) schematic of the follower fuel assembly.
(a)
(b)Figure 2
HFRR reactor layout.


### 2.2. Description to Mathematical Models

#### 2.2.1. Geometric Modeling

GAMBIT was used to build the geometric model of the HFRR core and generate the grid, with FLUENT then being used for the calculations. The mesh was simplified by not including the upstream drive mechanisms and the lower plenum. The heavy water tank was also ignored in the simulations because the heavy water tank has a separate circulation system which is not directly involved in the cooling of the fuel assemblies. The coolant flow in the inverse flux trap was also ignored with the core center modeled as a regular solid because the fuel assemblies are mainly cooled by the coolant flow through the narrow rectangular channels in the fuel assemblies. The interfaces between the fuel assemblies were assumed to be adiabatic in the mesh described in Section2.2.4. Therefore, the coolant channels between fuel assemblies were also ignored. These simplifications significantly reduced the simulation complexity without introducing large errors so that the heat transfer in the reactor core was reasonably predicted.The 18 standard fuel assemblies have exactly the same geometries, so they will have the same flow resistance characteristics. The 6 follower fuel assemblies were also assumed to have the same flow characteristics. The main purpose of this study was to show whether the flow was homogeneous throughout the core instead of the exact flow field inside each assembly. Therefore, each fuel assembly was treated as a virtual channel with the same flow resistance characteristics as the real assembly. The virtual assemblies were then modeled by a porous media model with the same boundary conditions as the real assembly.The virtual fuel assembly was modeled by the geometric model shown in Figure3. The upper blue region simulates the 0.35 m long core entrance section. The lower blue region simulates the 1 m long exit region where the coolant flows out of the active zone before entering the lower plenum. Each standard fuel assembly and follower fuel assembly had a rectangular cross section. The active zone in the standard fuel assemblies is marked in yellow with the nonfuel part in grey. The fuel region in the two follower elements completely inserted into the core is in deep red while that region in the four partly inserted follower elements is purple. The HFRR parameters were obtained from Gong et al. [9].Figure 3
Geometry of HFRR fuel assemblies: (a) overview of the geometry model; (b) detailed structure of the virtual fuel assemblies.


#### 2.2.2. Mathematical Model

As described in the FLUENT 14 user manual [13], the porous media model is an approximation method. In this model, the computational domain containing both fluid and solid regions is treated as a fluid region with a momentum source term added to the standard momentum equation to represent the flow losses due to the porous media. This momentum source term is expressed as(1)Si=-∑j=13Dijμvj+∑j=13Cij12ρvvj,where Si is the momentum source term in the momentum equations along the I axis (X, Y, or Z), ρ is the density, μ is the viscosity, v is the velocity, and Dij and Cij are the coefficient matrices. The momentum source term consists of a viscous loss term (the first term) and an inertial loss term (the second term). For simple homogeneous porous media, (1) can be simplified to(2)Si=-μαvi+C212ρvvi,where α is the permeability and C2 is the inertial loss coefficient. The two coefficients in (2) are the diagonal elements of D and C when they are simplified to diagonal matrices.Whether the porous media model can reasonably simulate the actual flow distribution largely depends on the loss coefficients. The loss coefficients are often measured in hydraulic flow experiments to get the fitting resistance-velocity relationships. However, there have not been any flow tests in the HFRR yet, so the fitting pressure drop-inlet coolant velocity relationship was determined based on the hydraulic experiments and calculations for CARR [14].Standard fuel assembly:(3)Δp=-183.58822v+3098.73222v2.Follower fuel assembly:(4)Δp=8029.13485v+8930.8323v2.The relationship between the pressure drop and the momentum source term was used to calculate the viscous loss coefficients and the inertial loss coefficients for the fuel assemblies:(5)Δp=-SiΔn,where Δn is the medium thickness along the I axis.As Figure1(b) shows, the flow channels in the plate-type fuel assemblies are really narrow, which leads to little cross flow in the radial directions. Therefore, the resistance coefficients in the other two directions were specified to be 1000 times the axial flow resistance coefficients to suppress the radial flow.The temperature distribution in the HFRR core was then found by solving the porous media energy equations for the thermal nonequilibrium model equations. The porous media model assumes that a solid region overlaps the fluid region. The fluid and solid regions are then connected by a surface heat transfer coefficient. The energy conservation equation for the fluid is(6)∂∂tγρfEf+∇·v→ρfEf+p=∇·γkf∇Tf-∑ihiJi+τ--·v→+Sfh+hfsAfsTs-Tf.The energy conservation equation for the solid is(7)∂∂t1-γρsEs=∇·1-γks∇Ts+Ssh+hfsAfsTf-Ts,where Tf is the fluid temperature, Ts is the solid temperature, hfs is the volumetric heat transfer coefficient between the fluid and solid, Afs is the interfacial area concentration, that is, the ratio of the actual fluid-solid interfacial area in the fuel assembly to the volume of the porous media in each cell, and Ssh and Sfh are the enthalpy source terms in the solid and fluid regions. The present study had no source term in the fluid. The last terms in (6) and (7) are the energy source terms in the thermal nonequilibrium model, which represent the energy transfer between the fluid and the solid.S
s
h, hfs, and Afs have to be specified before the calculation. In this study, Ssh was specified by fitting the axial energy density distribution based on previous model results [2]. Then, Ssh was input into the solid regions of the assemblies by a User Defined Function (UDF). Sudo et al. [4] and Ma et al. [18] demonstrated through experiments that the Dittus-Boelter turbulence equation [19] is applied to flow up or down narrow rectangular channels, as long as the Reynolds number is higher than 104. Therefore, hfs was specified by(8)Nu=0.023Re0.8Pr0.4,where the diameter and Afs were specified based on the specific geometry of the standard and follower fuel assemblies.


#### 2.2.3. Boundary Conditions and Model Options

Table1 summarizes the boundary conditions and the model options used in the HFRR simulations. In the HFRR, the solid material occupied a large part of the assemblies, so the porosity was quite small and the turbulence needed to be suppressed in the whole porous media region by enabling Laminar Zone Option in FLUENT. The inlet condition for the computational domain was the mass flow inlet condition with different mass flow rates. The outlet used the pressure outlet condition with the outlet gauge pressure set as 0.Table 1
Boundary conditions and model options.
Boundary conditions and model options
Enabled options
Turbulence model
Standardk-ɛ model
Laminar Zone Option
Enabled
Inlet boundary condition
Mass flow inlet
Outlet boundary condition
Pressure outlet
Energy source term
UDF
Discretization scheme
1st-order upwind
Pressure-velocity coupling
SimpleThe initial conditions were based on the specific thermal-hydraulic parameters for the HFRR presented in Gong et al. [9]. The solid material in the porous media region contained the fuel and the cladding which could not be distinguished after the homogenization in the porous media model. Therefore, the solid material properties were set as the volume average of the fuel and cladding properties. In addition, the model also considered the changes in the coolant density and viscosity with temperature. The HFRR neutronics calculations using RMC [2] were fit with a polynomial to describe the axial power distribution in each fuel assembly that was specified in a UDF.Figure4 shows the average power and the peak factor for each fuel assembly for the HFRR criticality conditions [9]. The assemblies with the highest average power for criticality are indicated by the solid line (standard fuel assembly) and dashed line (follower fuel assembly) in Figure 4. Numbers 7 and 12 standard fuel assemblies have the highest average powers of 1.33 MW. Numbers 21 and 22 follower fuel assemblies located in the active zone have the highest average follower fuel assembly power of 0.993 MW, as shown in Figure 4(b). The temperatures in these two follower fuel assemblies are expected to be high.Figure 4
Fuel assembly number and average power and peak factor of each fuel assembly in the HFRR [9]: (a) fuel assembly number; (b) average power and peak factors.
(a)
(b)


#### 2.2.4. Mesh Generation

A structured hexahedral mesh was used in this study. A grid independence study was conducted to determine the number of elements for grid independent results with several meshes with varying refinements used to obtain the optimal number of cells for the simulation. Figure5 shows the coolant flow rates in channel 1 used to judge the convergence for the various meshes. The coolant flow rate in channel 1 is around 0.0165 m3/s for the meshes with more than 2,977,478 elements. Despite the differences in the number of cells, the differences in the flow rate were quite small. This is because the porous media model makes the flow uniform, that is, almost independent of the number of cells. The mesh with 2,977,478 cells was used to balance the calculation time and accuracy.Figure 5
Volumetric flow rate in channel 1 for various numbers of elements.


### 2.3. Model Validation

The porous media model parameter settings must be reasonable for the simulation to be accurate. The empirical formulas for the resistance coefficients came from the experimental data of CARR [14] with the thermal nonequilibrium model used to solve the energy equation for the temperature distribution. The feasibilities of these options were verified first.


#### 2.3.1. Validation of the Resistance Coefficients

The resistance coefficients were validated based on the flow distribution calculations in CARR [14], whose experimental data is regarded as good reference data. As Tian et al. [6] described, the CARR core contains 17 standard fuel assemblies and 4 follower fuel assemblies. The numerical model of CARR shown in Figure 6 was similar to the model used in the present study. Only the fuel region together with its upstream and downstream sections was considered. The mesh was generated in the same way as for the entire core simulation of the flow distribution in Liu et al. [14]. The simulation used the standard k-ε turbulence model, a pressure outlet, and an inlet mass flow rate of 300 kg/s. The normalized flow distribution factor η was expressed as(9)ηi=Qi1/n∑j=1nQj,where Qi is the mass flow rate in each assembly.Figure 6
CARR simulation model: (a) assembly configuration of CARR; (b) outline of the geometry model.
(a)
(b)The predictions shown in Figure7 agree well with the experimental data with relative errors between −4.26% and 3.4%. The results indicate that the porous media model is reasonable and can be used for simulating the flow distribution.Figure 7
Flow distribution factor in each CARR fuel assembly.The errors between the simulation and the measurements come from the differences between the model and the actual core. For example, the model does not include some physical structures such as the control rods and their guide tubes. These physical structures will affect the flow field, induce mixing, and reduce the temperature stratification in the core. Moreover, the simplified decay tank at the core outlet will also affect the flow distribution symmetry.


#### 2.3.2. Validation of the Thermal Nonequilibrium Model

The porous media model energy equation can use the thermal equilibrium model or the thermal nonequilibrium model. The thermal equilibrium model assumes that the fluid and solid temperatures are the same. However, for these simulations of channels with high power generation rates in the solid, the thermal equilibrium model is not sufficiently accurate, especially when the temperature distribution of the solid needs to be accurately predicted. Since the solid is not treated separately in the thermal equilibrium model, there is no conduction equation for the solid that can be solved to calculate the temperature distribution.The thermal nonequilibrium model was validated in the present study by comparing the predicted core temperature distributions from the thermal nonequilibrium model, the thermal equilibrium model, and the exact model (Gong et al. [9]). The exact model does not apply the porous media approach. It assumes that the coolant flow rate inside the standard fuel assembly is about 5.5 m/s, based on the calculated mass flow rate in Section 3.1. The boundary conditions and model options of the exact model were specified according to Gong et al. [9]. The three models did not consider the power differences between the fuel assemblies with the energy source terms in the standard and follower fuel assemblies specified as their average power generation rates. The predicted temperatures are shown in Figures 8 and 9 and Table 2. Figure 8 shows the temperature contours of a section at the edge of the heated region and the whole calculation region obtained by the thermal equilibrium model. Figure 9 shows the temperature distributions on a vertical plane and in the whole core region presented by the thermal nonequilibrium model. The peak fuel temperatures are listed in Table 2, with the result of the exact model used as the reference.Table 2
Peak fuel temperature.
Model
Peak fuel temperature (K)
Relative error/%
Exact model
345
/
Thermal nonequilibrium model
344
−0.289
Thermal equilibrium model
375
+8.695Figure 8
Temperature contours predicted by the thermal equilibrium model: (a) a section at the edge of the heat region and (b) the whole calculation region.
(a)
(b)Figure 9
Temperature contours predicted by the thermal nonequilibrium model: (a) whole calculation region, (b) vertical cross sections for coolant (left) and fuel (right) temperatures, respectively, and (c) fuel transverse cross section with the peak temperature.
(a)
(b)
(c)The peak fuel temperatures in Table2 predicted by the thermal nonequilibrium model agree well with those by the exact model, with the thermal equilibrium model giving larger differences. The thermal equilibrium model assumes that the fluid and solid in a computational cell are in thermal equilibrium; hence, they have the same temperatures. As shown in Figure 8, the temperatures gradually increase along the flow direction, even with the cosine distribution of the input power. The fluid and solid temperatures are quite different in the thermal nonequilibrium model. As shown in Figure 9, the fuel temperature is higher in the middle and lower at both ends along the axial direction, with the fluid temperature gradually increasing along the flow direction. Thus, the thermal nonequilibrium model provides more reasonable results for this study.Figure9(c) shows the temperatures in the fuel cross section with the peak fuel temperature. Since there are no power differences between the fuel assemblies in these models, the peak temperatures in each standard fuel assembly are the same. With the overlapping solid mesh in the porous media model, each fuel assembly is essentially independent with what can be considered to be adiabatic interfaces since the temperatures are the same. Therefore, the temperature differences between the standard and follower fuel assemblies are obvious, as shown in Figure 9. There are also coolant channels between the fuel assemblies that absorb some heat by convection. However, the heat exchange in these channels was ignored because the convection in these channels was significantly lower than the conduction and convection within the assemblies, especially for these plate-type fuel assemblies. The gaps between the assemblies were also ignored in the simulations. If the interfaces were not adiabatic, the conduction would affect the temperature distribution. Therefore, the adiabatic option is reasonable.


## 2.1. HFRR Overview

As Gong et al. [9] described, HFRR is a pool-type research reactor which is moderated and cooled by light water and reflected by beryllium and heavy water, with a maximum thermal power of 20 MW and low-enriched uranium (LEU) plate-type fuel elements. Figure 1 shows the fuel element configuration and the schematic of the follower fuel assembly. The core contains 18 standard fuel elements in yellow, 6 follower fuel elements in grey, and 4 irradiation tubes in blue. The fuel meat is U3Si2-Al and the fuel cladding is Al (6061). Gong et al. [9] showed more details about the structure of the standard and follower fuel elements in their paper. The core center has an inverse flux trap made of Be which has a peak thermal neutron flux of 4.04 × 1014 n/cm−2s. Figure 2 shows the 3D schematic of the reactor core as well as the sketch of the HFRR layout. Besides the fuel elements, the core also contains a heavy water tank installed around the core and a lower plenum. During normal operation, the coolant flows down the core by forced convection. The inlet coolant temperature is 35°C and the working pressure is 0.152 MPa.Figure 1
HFRR core configuration and schematic of follower fuel assembly [9]: (a) fuel element configuration; (b) schematic of the follower fuel assembly.
(a)
(b)Figure 2
HFRR reactor layout.


## 2.2. Description to Mathematical Models

### 2.2.1. Geometric Modeling

GAMBIT was used to build the geometric model of the HFRR core and generate the grid, with FLUENT then being used for the calculations. The mesh was simplified by not including the upstream drive mechanisms and the lower plenum. The heavy water tank was also ignored in the simulations because the heavy water tank has a separate circulation system which is not directly involved in the cooling of the fuel assemblies. The coolant flow in the inverse flux trap was also ignored with the core center modeled as a regular solid because the fuel assemblies are mainly cooled by the coolant flow through the narrow rectangular channels in the fuel assemblies. The interfaces between the fuel assemblies were assumed to be adiabatic in the mesh described in Section2.2.4. Therefore, the coolant channels between fuel assemblies were also ignored. These simplifications significantly reduced the simulation complexity without introducing large errors so that the heat transfer in the reactor core was reasonably predicted.The 18 standard fuel assemblies have exactly the same geometries, so they will have the same flow resistance characteristics. The 6 follower fuel assemblies were also assumed to have the same flow characteristics. The main purpose of this study was to show whether the flow was homogeneous throughout the core instead of the exact flow field inside each assembly. Therefore, each fuel assembly was treated as a virtual channel with the same flow resistance characteristics as the real assembly. The virtual assemblies were then modeled by a porous media model with the same boundary conditions as the real assembly.The virtual fuel assembly was modeled by the geometric model shown in Figure3. The upper blue region simulates the 0.35 m long core entrance section. The lower blue region simulates the 1 m long exit region where the coolant flows out of the active zone before entering the lower plenum. Each standard fuel assembly and follower fuel assembly had a rectangular cross section. The active zone in the standard fuel assemblies is marked in yellow with the nonfuel part in grey. The fuel region in the two follower elements completely inserted into the core is in deep red while that region in the four partly inserted follower elements is purple. The HFRR parameters were obtained from Gong et al. [9].Figure 3
Geometry of HFRR fuel assemblies: (a) overview of the geometry model; (b) detailed structure of the virtual fuel assemblies.


### 2.2.2. Mathematical Model

As described in the FLUENT 14 user manual [13], the porous media model is an approximation method. In this model, the computational domain containing both fluid and solid regions is treated as a fluid region with a momentum source term added to the standard momentum equation to represent the flow losses due to the porous media. This momentum source term is expressed as(1)Si=-∑j=13Dijμvj+∑j=13Cij12ρvvj,where Si is the momentum source term in the momentum equations along the I axis (X, Y, or Z), ρ is the density, μ is the viscosity, v is the velocity, and Dij and Cij are the coefficient matrices. The momentum source term consists of a viscous loss term (the first term) and an inertial loss term (the second term). For simple homogeneous porous media, (1) can be simplified to(2)Si=-μαvi+C212ρvvi,where α is the permeability and C2 is the inertial loss coefficient. The two coefficients in (2) are the diagonal elements of D and C when they are simplified to diagonal matrices.Whether the porous media model can reasonably simulate the actual flow distribution largely depends on the loss coefficients. The loss coefficients are often measured in hydraulic flow experiments to get the fitting resistance-velocity relationships. However, there have not been any flow tests in the HFRR yet, so the fitting pressure drop-inlet coolant velocity relationship was determined based on the hydraulic experiments and calculations for CARR [14].Standard fuel assembly:(3)Δp=-183.58822v+3098.73222v2.Follower fuel assembly:(4)Δp=8029.13485v+8930.8323v2.The relationship between the pressure drop and the momentum source term was used to calculate the viscous loss coefficients and the inertial loss coefficients for the fuel assemblies:(5)Δp=-SiΔn,where Δn is the medium thickness along the I axis.As Figure1(b) shows, the flow channels in the plate-type fuel assemblies are really narrow, which leads to little cross flow in the radial directions. Therefore, the resistance coefficients in the other two directions were specified to be 1000 times the axial flow resistance coefficients to suppress the radial flow.The temperature distribution in the HFRR core was then found by solving the porous media energy equations for the thermal nonequilibrium model equations. The porous media model assumes that a solid region overlaps the fluid region. The fluid and solid regions are then connected by a surface heat transfer coefficient. The energy conservation equation for the fluid is(6)∂∂tγρfEf+∇·v→ρfEf+p=∇·γkf∇Tf-∑ihiJi+τ--·v→+Sfh+hfsAfsTs-Tf.The energy conservation equation for the solid is(7)∂∂t1-γρsEs=∇·1-γks∇Ts+Ssh+hfsAfsTf-Ts,where Tf is the fluid temperature, Ts is the solid temperature, hfs is the volumetric heat transfer coefficient between the fluid and solid, Afs is the interfacial area concentration, that is, the ratio of the actual fluid-solid interfacial area in the fuel assembly to the volume of the porous media in each cell, and Ssh and Sfh are the enthalpy source terms in the solid and fluid regions. The present study had no source term in the fluid. The last terms in (6) and (7) are the energy source terms in the thermal nonequilibrium model, which represent the energy transfer between the fluid and the solid.S
s
h, hfs, and Afs have to be specified before the calculation. In this study, Ssh was specified by fitting the axial energy density distribution based on previous model results [2]. Then, Ssh was input into the solid regions of the assemblies by a User Defined Function (UDF). Sudo et al. [4] and Ma et al. [18] demonstrated through experiments that the Dittus-Boelter turbulence equation [19] is applied to flow up or down narrow rectangular channels, as long as the Reynolds number is higher than 104. Therefore, hfs was specified by(8)Nu=0.023Re0.8Pr0.4,where the diameter and Afs were specified based on the specific geometry of the standard and follower fuel assemblies.


### 2.2.3. Boundary Conditions and Model Options

Table1 summarizes the boundary conditions and the model options used in the HFRR simulations. In the HFRR, the solid material occupied a large part of the assemblies, so the porosity was quite small and the turbulence needed to be suppressed in the whole porous media region by enabling Laminar Zone Option in FLUENT. The inlet condition for the computational domain was the mass flow inlet condition with different mass flow rates. The outlet used the pressure outlet condition with the outlet gauge pressure set as 0.Table 1
Boundary conditions and model options.
Boundary conditions and model options
Enabled options
Turbulence model
Standardk-ɛ model
Laminar Zone Option
Enabled
Inlet boundary condition
Mass flow inlet
Outlet boundary condition
Pressure outlet
Energy source term
UDF
Discretization scheme
1st-order upwind
Pressure-velocity coupling
SimpleThe initial conditions were based on the specific thermal-hydraulic parameters for the HFRR presented in Gong et al. [9]. The solid material in the porous media region contained the fuel and the cladding which could not be distinguished after the homogenization in the porous media model. Therefore, the solid material properties were set as the volume average of the fuel and cladding properties. In addition, the model also considered the changes in the coolant density and viscosity with temperature. The HFRR neutronics calculations using RMC [2] were fit with a polynomial to describe the axial power distribution in each fuel assembly that was specified in a UDF.Figure4 shows the average power and the peak factor for each fuel assembly for the HFRR criticality conditions [9]. The assemblies with the highest average power for criticality are indicated by the solid line (standard fuel assembly) and dashed line (follower fuel assembly) in Figure 4. Numbers 7 and 12 standard fuel assemblies have the highest average powers of 1.33 MW. Numbers 21 and 22 follower fuel assemblies located in the active zone have the highest average follower fuel assembly power of 0.993 MW, as shown in Figure 4(b). The temperatures in these two follower fuel assemblies are expected to be high.Figure 4
Fuel assembly number and average power and peak factor of each fuel assembly in the HFRR [9]: (a) fuel assembly number; (b) average power and peak factors.
(a)
(b)


### 2.2.4. Mesh Generation

A structured hexahedral mesh was used in this study. A grid independence study was conducted to determine the number of elements for grid independent results with several meshes with varying refinements used to obtain the optimal number of cells for the simulation. Figure5 shows the coolant flow rates in channel 1 used to judge the convergence for the various meshes. The coolant flow rate in channel 1 is around 0.0165 m3/s for the meshes with more than 2,977,478 elements. Despite the differences in the number of cells, the differences in the flow rate were quite small. This is because the porous media model makes the flow uniform, that is, almost independent of the number of cells. The mesh with 2,977,478 cells was used to balance the calculation time and accuracy.Figure 5
Volumetric flow rate in channel 1 for various numbers of elements.


## 2.2.1. Geometric Modeling

GAMBIT was used to build the geometric model of the HFRR core and generate the grid, with FLUENT then being used for the calculations. The mesh was simplified by not including the upstream drive mechanisms and the lower plenum. The heavy water tank was also ignored in the simulations because the heavy water tank has a separate circulation system which is not directly involved in the cooling of the fuel assemblies. The coolant flow in the inverse flux trap was also ignored with the core center modeled as a regular solid because the fuel assemblies are mainly cooled by the coolant flow through the narrow rectangular channels in the fuel assemblies. The interfaces between the fuel assemblies were assumed to be adiabatic in the mesh described in Section2.2.4. Therefore, the coolant channels between fuel assemblies were also ignored. These simplifications significantly reduced the simulation complexity without introducing large errors so that the heat transfer in the reactor core was reasonably predicted.The 18 standard fuel assemblies have exactly the same geometries, so they will have the same flow resistance characteristics. The 6 follower fuel assemblies were also assumed to have the same flow characteristics. The main purpose of this study was to show whether the flow was homogeneous throughout the core instead of the exact flow field inside each assembly. Therefore, each fuel assembly was treated as a virtual channel with the same flow resistance characteristics as the real assembly. The virtual assemblies were then modeled by a porous media model with the same boundary conditions as the real assembly.The virtual fuel assembly was modeled by the geometric model shown in Figure3. The upper blue region simulates the 0.35 m long core entrance section. The lower blue region simulates the 1 m long exit region where the coolant flows out of the active zone before entering the lower plenum. Each standard fuel assembly and follower fuel assembly had a rectangular cross section. The active zone in the standard fuel assemblies is marked in yellow with the nonfuel part in grey. The fuel region in the two follower elements completely inserted into the core is in deep red while that region in the four partly inserted follower elements is purple. The HFRR parameters were obtained from Gong et al. [9].Figure 3
Geometry of HFRR fuel assemblies: (a) overview of the geometry model; (b) detailed structure of the virtual fuel assemblies.


## 2.2.2. Mathematical Model

As described in the FLUENT 14 user manual [13], the porous media model is an approximation method. In this model, the computational domain containing both fluid and solid regions is treated as a fluid region with a momentum source term added to the standard momentum equation to represent the flow losses due to the porous media. This momentum source term is expressed as(1)Si=-∑j=13Dijμvj+∑j=13Cij12ρvvj,where Si is the momentum source term in the momentum equations along the I axis (X, Y, or Z), ρ is the density, μ is the viscosity, v is the velocity, and Dij and Cij are the coefficient matrices. The momentum source term consists of a viscous loss term (the first term) and an inertial loss term (the second term). For simple homogeneous porous media, (1) can be simplified to(2)Si=-μαvi+C212ρvvi,where α is the permeability and C2 is the inertial loss coefficient. The two coefficients in (2) are the diagonal elements of D and C when they are simplified to diagonal matrices.Whether the porous media model can reasonably simulate the actual flow distribution largely depends on the loss coefficients. The loss coefficients are often measured in hydraulic flow experiments to get the fitting resistance-velocity relationships. However, there have not been any flow tests in the HFRR yet, so the fitting pressure drop-inlet coolant velocity relationship was determined based on the hydraulic experiments and calculations for CARR [14].Standard fuel assembly:(3)Δp=-183.58822v+3098.73222v2.Follower fuel assembly:(4)Δp=8029.13485v+8930.8323v2.The relationship between the pressure drop and the momentum source term was used to calculate the viscous loss coefficients and the inertial loss coefficients for the fuel assemblies:(5)Δp=-SiΔn,where Δn is the medium thickness along the I axis.As Figure1(b) shows, the flow channels in the plate-type fuel assemblies are really narrow, which leads to little cross flow in the radial directions. Therefore, the resistance coefficients in the other two directions were specified to be 1000 times the axial flow resistance coefficients to suppress the radial flow.The temperature distribution in the HFRR core was then found by solving the porous media energy equations for the thermal nonequilibrium model equations. The porous media model assumes that a solid region overlaps the fluid region. The fluid and solid regions are then connected by a surface heat transfer coefficient. The energy conservation equation for the fluid is(6)∂∂tγρfEf+∇·v→ρfEf+p=∇·γkf∇Tf-∑ihiJi+τ--·v→+Sfh+hfsAfsTs-Tf.The energy conservation equation for the solid is(7)∂∂t1-γρsEs=∇·1-γks∇Ts+Ssh+hfsAfsTf-Ts,where Tf is the fluid temperature, Ts is the solid temperature, hfs is the volumetric heat transfer coefficient between the fluid and solid, Afs is the interfacial area concentration, that is, the ratio of the actual fluid-solid interfacial area in the fuel assembly to the volume of the porous media in each cell, and Ssh and Sfh are the enthalpy source terms in the solid and fluid regions. The present study had no source term in the fluid. The last terms in (6) and (7) are the energy source terms in the thermal nonequilibrium model, which represent the energy transfer between the fluid and the solid.S
s
h, hfs, and Afs have to be specified before the calculation. In this study, Ssh was specified by fitting the axial energy density distribution based on previous model results [2]. Then, Ssh was input into the solid regions of the assemblies by a User Defined Function (UDF). Sudo et al. [4] and Ma et al. [18] demonstrated through experiments that the Dittus-Boelter turbulence equation [19] is applied to flow up or down narrow rectangular channels, as long as the Reynolds number is higher than 104. Therefore, hfs was specified by(8)Nu=0.023Re0.8Pr0.4,where the diameter and Afs were specified based on the specific geometry of the standard and follower fuel assemblies.


## 2.2.3. Boundary Conditions and Model Options

Table1 summarizes the boundary conditions and the model options used in the HFRR simulations. In the HFRR, the solid material occupied a large part of the assemblies, so the porosity was quite small and the turbulence needed to be suppressed in the whole porous media region by enabling Laminar Zone Option in FLUENT. The inlet condition for the computational domain was the mass flow inlet condition with different mass flow rates. The outlet used the pressure outlet condition with the outlet gauge pressure set as 0.Table 1
Boundary conditions and model options.
Boundary conditions and model options
Enabled options
Turbulence model
Standardk-ɛ model
Laminar Zone Option
Enabled
Inlet boundary condition
Mass flow inlet
Outlet boundary condition
Pressure outlet
Energy source term
UDF
Discretization scheme
1st-order upwind
Pressure-velocity coupling
SimpleThe initial conditions were based on the specific thermal-hydraulic parameters for the HFRR presented in Gong et al. [9]. The solid material in the porous media region contained the fuel and the cladding which could not be distinguished after the homogenization in the porous media model. Therefore, the solid material properties were set as the volume average of the fuel and cladding properties. In addition, the model also considered the changes in the coolant density and viscosity with temperature. The HFRR neutronics calculations using RMC [2] were fit with a polynomial to describe the axial power distribution in each fuel assembly that was specified in a UDF.Figure4 shows the average power and the peak factor for each fuel assembly for the HFRR criticality conditions [9]. The assemblies with the highest average power for criticality are indicated by the solid line (standard fuel assembly) and dashed line (follower fuel assembly) in Figure 4. Numbers 7 and 12 standard fuel assemblies have the highest average powers of 1.33 MW. Numbers 21 and 22 follower fuel assemblies located in the active zone have the highest average follower fuel assembly power of 0.993 MW, as shown in Figure 4(b). The temperatures in these two follower fuel assemblies are expected to be high.Figure 4
Fuel assembly number and average power and peak factor of each fuel assembly in the HFRR [9]: (a) fuel assembly number; (b) average power and peak factors.
(a)
(b)


## 2.2.4. Mesh Generation

A structured hexahedral mesh was used in this study. A grid independence study was conducted to determine the number of elements for grid independent results with several meshes with varying refinements used to obtain the optimal number of cells for the simulation. Figure5 shows the coolant flow rates in channel 1 used to judge the convergence for the various meshes. The coolant flow rate in channel 1 is around 0.0165 m3/s for the meshes with more than 2,977,478 elements. Despite the differences in the number of cells, the differences in the flow rate were quite small. This is because the porous media model makes the flow uniform, that is, almost independent of the number of cells. The mesh with 2,977,478 cells was used to balance the calculation time and accuracy.Figure 5
Volumetric flow rate in channel 1 for various numbers of elements.


## 2.3. Model Validation

The porous media model parameter settings must be reasonable for the simulation to be accurate. The empirical formulas for the resistance coefficients came from the experimental data of CARR [14] with the thermal nonequilibrium model used to solve the energy equation for the temperature distribution. The feasibilities of these options were verified first.


### 2.3.1. Validation of the Resistance Coefficients

The resistance coefficients were validated based on the flow distribution calculations in CARR [14], whose experimental data is regarded as good reference data. As Tian et al. [6] described, the CARR core contains 17 standard fuel assemblies and 4 follower fuel assemblies. The numerical model of CARR shown in Figure 6 was similar to the model used in the present study. Only the fuel region together with its upstream and downstream sections was considered. The mesh was generated in the same way as for the entire core simulation of the flow distribution in Liu et al. [14]. The simulation used the standard k-ε turbulence model, a pressure outlet, and an inlet mass flow rate of 300 kg/s. The normalized flow distribution factor η was expressed as(9)ηi=Qi1/n∑j=1nQj,where Qi is the mass flow rate in each assembly.Figure 6
CARR simulation model: (a) assembly configuration of CARR; (b) outline of the geometry model.
(a)
(b)The predictions shown in Figure7 agree well with the experimental data with relative errors between −4.26% and 3.4%. The results indicate that the porous media model is reasonable and can be used for simulating the flow distribution.Figure 7
Flow distribution factor in each CARR fuel assembly.The errors between the simulation and the measurements come from the differences between the model and the actual core. For example, the model does not include some physical structures such as the control rods and their guide tubes. These physical structures will affect the flow field, induce mixing, and reduce the temperature stratification in the core. Moreover, the simplified decay tank at the core outlet will also affect the flow distribution symmetry.


### 2.3.2. Validation of the Thermal Nonequilibrium Model

The porous media model energy equation can use the thermal equilibrium model or the thermal nonequilibrium model. The thermal equilibrium model assumes that the fluid and solid temperatures are the same. However, for these simulations of channels with high power generation rates in the solid, the thermal equilibrium model is not sufficiently accurate, especially when the temperature distribution of the solid needs to be accurately predicted. Since the solid is not treated separately in the thermal equilibrium model, there is no conduction equation for the solid that can be solved to calculate the temperature distribution.The thermal nonequilibrium model was validated in the present study by comparing the predicted core temperature distributions from the thermal nonequilibrium model, the thermal equilibrium model, and the exact model (Gong et al. [9]). The exact model does not apply the porous media approach. It assumes that the coolant flow rate inside the standard fuel assembly is about 5.5 m/s, based on the calculated mass flow rate in Section 3.1. The boundary conditions and model options of the exact model were specified according to Gong et al. [9]. The three models did not consider the power differences between the fuel assemblies with the energy source terms in the standard and follower fuel assemblies specified as their average power generation rates. The predicted temperatures are shown in Figures 8 and 9 and Table 2. Figure 8 shows the temperature contours of a section at the edge of the heated region and the whole calculation region obtained by the thermal equilibrium model. Figure 9 shows the temperature distributions on a vertical plane and in the whole core region presented by the thermal nonequilibrium model. The peak fuel temperatures are listed in Table 2, with the result of the exact model used as the reference.Table 2
Peak fuel temperature.
Model
Peak fuel temperature (K)
Relative error/%
Exact model
345
/
Thermal nonequilibrium model
344
−0.289
Thermal equilibrium model
375
+8.695Figure 8
Temperature contours predicted by the thermal equilibrium model: (a) a section at the edge of the heat region and (b) the whole calculation region.
(a)
(b)Figure 9
Temperature contours predicted by the thermal nonequilibrium model: (a) whole calculation region, (b) vertical cross sections for coolant (left) and fuel (right) temperatures, respectively, and (c) fuel transverse cross section with the peak temperature.
(a)
(b)
(c)The peak fuel temperatures in Table2 predicted by the thermal nonequilibrium model agree well with those by the exact model, with the thermal equilibrium model giving larger differences. The thermal equilibrium model assumes that the fluid and solid in a computational cell are in thermal equilibrium; hence, they have the same temperatures. As shown in Figure 8, the temperatures gradually increase along the flow direction, even with the cosine distribution of the input power. The fluid and solid temperatures are quite different in the thermal nonequilibrium model. As shown in Figure 9, the fuel temperature is higher in the middle and lower at both ends along the axial direction, with the fluid temperature gradually increasing along the flow direction. Thus, the thermal nonequilibrium model provides more reasonable results for this study.Figure9(c) shows the temperatures in the fuel cross section with the peak fuel temperature. Since there are no power differences between the fuel assemblies in these models, the peak temperatures in each standard fuel assembly are the same. With the overlapping solid mesh in the porous media model, each fuel assembly is essentially independent with what can be considered to be adiabatic interfaces since the temperatures are the same. Therefore, the temperature differences between the standard and follower fuel assemblies are obvious, as shown in Figure 9. There are also coolant channels between the fuel assemblies that absorb some heat by convection. However, the heat exchange in these channels was ignored because the convection in these channels was significantly lower than the conduction and convection within the assemblies, especially for these plate-type fuel assemblies. The gaps between the assemblies were also ignored in the simulations. If the interfaces were not adiabatic, the conduction would affect the temperature distribution. Therefore, the adiabatic option is reasonable.


## 2.3.1. Validation of the Resistance Coefficients

The resistance coefficients were validated based on the flow distribution calculations in CARR [14], whose experimental data is regarded as good reference data. As Tian et al. [6] described, the CARR core contains 17 standard fuel assemblies and 4 follower fuel assemblies. The numerical model of CARR shown in Figure 6 was similar to the model used in the present study. Only the fuel region together with its upstream and downstream sections was considered. The mesh was generated in the same way as for the entire core simulation of the flow distribution in Liu et al. [14]. The simulation used the standard k-ε turbulence model, a pressure outlet, and an inlet mass flow rate of 300 kg/s. The normalized flow distribution factor η was expressed as(9)ηi=Qi1/n∑j=1nQj,where Qi is the mass flow rate in each assembly.Figure 6
CARR simulation model: (a) assembly configuration of CARR; (b) outline of the geometry model.
(a)
(b)The predictions shown in Figure7 agree well with the experimental data with relative errors between −4.26% and 3.4%. The results indicate that the porous media model is reasonable and can be used for simulating the flow distribution.Figure 7
Flow distribution factor in each CARR fuel assembly.The errors between the simulation and the measurements come from the differences between the model and the actual core. For example, the model does not include some physical structures such as the control rods and their guide tubes. These physical structures will affect the flow field, induce mixing, and reduce the temperature stratification in the core. Moreover, the simplified decay tank at the core outlet will also affect the flow distribution symmetry.


## 2.3.2. Validation of the Thermal Nonequilibrium Model

The porous media model energy equation can use the thermal equilibrium model or the thermal nonequilibrium model. The thermal equilibrium model assumes that the fluid and solid temperatures are the same. However, for these simulations of channels with high power generation rates in the solid, the thermal equilibrium model is not sufficiently accurate, especially when the temperature distribution of the solid needs to be accurately predicted. Since the solid is not treated separately in the thermal equilibrium model, there is no conduction equation for the solid that can be solved to calculate the temperature distribution.The thermal nonequilibrium model was validated in the present study by comparing the predicted core temperature distributions from the thermal nonequilibrium model, the thermal equilibrium model, and the exact model (Gong et al. [9]). The exact model does not apply the porous media approach. It assumes that the coolant flow rate inside the standard fuel assembly is about 5.5 m/s, based on the calculated mass flow rate in Section 3.1. The boundary conditions and model options of the exact model were specified according to Gong et al. [9]. The three models did not consider the power differences between the fuel assemblies with the energy source terms in the standard and follower fuel assemblies specified as their average power generation rates. The predicted temperatures are shown in Figures 8 and 9 and Table 2. Figure 8 shows the temperature contours of a section at the edge of the heated region and the whole calculation region obtained by the thermal equilibrium model. Figure 9 shows the temperature distributions on a vertical plane and in the whole core region presented by the thermal nonequilibrium model. The peak fuel temperatures are listed in Table 2, with the result of the exact model used as the reference.Table 2
Peak fuel temperature.
Model
Peak fuel temperature (K)
Relative error/%
Exact model
345
/
Thermal nonequilibrium model
344
−0.289
Thermal equilibrium model
375
+8.695Figure 8
Temperature contours predicted by the thermal equilibrium model: (a) a section at the edge of the heat region and (b) the whole calculation region.
(a)
(b)Figure 9
Temperature contours predicted by the thermal nonequilibrium model: (a) whole calculation region, (b) vertical cross sections for coolant (left) and fuel (right) temperatures, respectively, and (c) fuel transverse cross section with the peak temperature.
(a)
(b)
(c)The peak fuel temperatures in Table2 predicted by the thermal nonequilibrium model agree well with those by the exact model, with the thermal equilibrium model giving larger differences. The thermal equilibrium model assumes that the fluid and solid in a computational cell are in thermal equilibrium; hence, they have the same temperatures. As shown in Figure 8, the temperatures gradually increase along the flow direction, even with the cosine distribution of the input power. The fluid and solid temperatures are quite different in the thermal nonequilibrium model. As shown in Figure 9, the fuel temperature is higher in the middle and lower at both ends along the axial direction, with the fluid temperature gradually increasing along the flow direction. Thus, the thermal nonequilibrium model provides more reasonable results for this study.Figure9(c) shows the temperatures in the fuel cross section with the peak fuel temperature. Since there are no power differences between the fuel assemblies in these models, the peak temperatures in each standard fuel assembly are the same. With the overlapping solid mesh in the porous media model, each fuel assembly is essentially independent with what can be considered to be adiabatic interfaces since the temperatures are the same. Therefore, the temperature differences between the standard and follower fuel assemblies are obvious, as shown in Figure 9. There are also coolant channels between the fuel assemblies that absorb some heat by convection. However, the heat exchange in these channels was ignored because the convection in these channels was significantly lower than the conduction and convection within the assemblies, especially for these plate-type fuel assemblies. The gaps between the assemblies were also ignored in the simulations. If the interfaces were not adiabatic, the conduction would affect the temperature distribution. Therefore, the adiabatic option is reasonable.


## 3. Results and Discussion

The porous media model with the nonequilibrium model was used to calculate the flow and temperature distributions in the full HFRR core. Simulations are carried out for five different inlet mass flow rates: 300 kg/s, 350 kg/s, 400 kg/s, 450 kg/s, and 500 kg/s. The predictions of the present model are compared with those of the exact model. This work provides insight into the thermal-hydraulic design of the HFRR and can be used for the HFRR project application.


### 3.1. Flow Distribution

Figure10 shows the flow distribution inside each assembly at an inlet mass flow rate of 300 kg/s in accordance with the order in Figure 4(a). No coolant flows through the core center due to the inverse flux trap. The normalized flow distribution factors of the standard and follower fuel assemblies are marked in Figure 10. Standard fuel assemblies 8 and 11 have the highest coolant mass flow rates with normalized flow distribution factors of 1.0267. The normalized flow distribution factors for the standard fuel assemblies near the follower fuel assemblies are around 0.985, almost 4% lower than the peak flow rate. The flow rates in the standard fuel assemblies vary by less than 3% from the average. Since the flow resistance coefficients in the follower fuel assemblies are much larger than those in the standard fuel assemblies, the coolant flow rates in the follower fuel assemblies are much lower with smaller variations. Therefore, the predictions indicate that the flow rates inside the fuel assemblies are quite uniform, which means that the inverse flux trap structure has little impact on the flow distribution.Figure 10
Flow distribution inside each assembly for an inlet mass flow rate of 300 kg/s.Figure11 shows the variations of the flow distribution factor in each standard fuel assembly for various inlet conditions. The flow field variations are similar but gradually increase with increasing inlet coolant mass flow rate. At the highest inlet mass flow rate of 450 kg/s, the largest difference between the flow rate in each assembly and the average flow rate is less than 5%. Thus, the results show that the flow distribution in the standard fuel assemblies is almost uniform and there will be sufficient coolant flow within each assembly.Figure 11
Variation of the flow distribution factor in each standard fuel assembly for various inlet flow rates.Figure12 shows the flow fields at the core entrance and outlet at an inlet mass flow rate of 300 kg/s. Since the flow resistance coefficients in the follower fuel assemblies are much larger than those in the standard fuel assemblies, more coolant enters the standard fuel assemblies with velocities 2 or 3 times those in the follower fuel assemblies. A large recirculation zone develops at the fuel assembly outlet with many vortices at the bottom of the core due to the inverse flux trap that blocks the flow in the middle. The coolant flowing out the assemblies is then fully mixed and the velocity is nearly homogeneous there. Thus, the inverse flux trap promotes mixing of the fluid streams exiting the core.Figure 12
Flow fields at the core entrance and outlet for an inlet mass flow rate of 300 kg/s.


### 3.2. Temperature Distribution

The temperature distribution in the core at an inlet mass flow rate of 400 kg/s is shown in Figure13 with only the solid temperature presented in the active core region and the coolant temperature presented in the upper and lower regions. Follower fuel assemblies 21 and 22 have the highest temperature of 395 K, with the highest coolant temperature rise through those assemblies, 12 K. The highest temperatures in the standard fuel assemblies occur at numbers 7 and 12 standard fuel assemblies. Their temperature is about 375 K and the temperature rise of the coolant through these assemblies is about 12 K. As shown in Section 2.2.3, numbers 21 and 22 follower fuel assemblies have the highest average power levels among the follower fuel assemblies. At the same time, the coolant flow rates through these two assemblies are relatively low, as mentioned in Section 3.1. Therefore, numbers 21 and 22 follower fuel assemblies are hottest due to the reduced cooling.Figure 13
Temperature distribution in the core at an inlet mass flow rate of 400 kg/s.Further studies were conducted to ensure the safety of the hottest standard and follower fuel assemblies. Fuel assemblies 7 and 22 are symmetric, and therefore only one was analyzed. The coolant inlet velocities to follower fuel assemblies 12 and 21 were based on the distribution results in Section3.1. Then, the exact model was used to calculate more detailed temperature distributions in these two assemblies.The standard assembly temperatures predicted by the porous media model and the exact model are shown in Figure14 with the follower fuel assembly temperatures shown in Figure 15. The two temperature profiles predicted by the porous media in each figure represent the average solid temperature and the average fluid temperature variations along the axis obtained by the nonequilibrium model, while the corresponding curves of the exact model represent the fuel temperature and the coolant temperature in the channel.Figure 14
Calculated temperatures in the number 12 standard fuel assembly.Figure 15
Calculated temperatures in number 21 follower fuel assembly.Figures14 and 15 show that the temperatures predicted by the porous media model agree well with those predicted by the exact model. The temperature variations predicted by the porous media model are somewhat flatter than those predicted by the exact model because the exact model exactly simulates the axial power distribution while the porous media model is an approximate method. In the porous media model, the heat transfer coefficient, hfs, between the fluid and the solid in the energy conservation equation was specified by the empirical Dittus-Boelter equation [19]. Moreover, hfs and the solid thermal conductivity were assumed to be constant along the whole channel. However, at the inlet and the outlet of the exact model, the flow field in the narrow rectangular channel is not fully developed compared to that modeled by the porous media method. This difference leads to less turbulence and decreases the heat transfer, so the power profiles show different shape at the inlet and outlet region. Nevertheless, these approximations led to small differences between the two models with a maximum difference of less than 1.4%, which indicates that the porous media model can reasonably predict the temperature distribution with the thermal nonequilibrium model. In addition, for simulations of the full reactor core, the porous media model has the additional advantage of significantly reducing the computational cost.Figure16 shows the variations of the peak temperatures in fuel assemblies 12 and 21 predicted by the exact model for various inlet conditions. As the inlet coolant mass flow rate increases from 300 kg/s to 500 kg/s, the peak temperature in number 12 standard fuel assembly gradually decreases from 350 K to 339 K, with all the temperatures lower than the local saturation temperature. Therefore, the predictions indicate that all the temperatures in the standard fuel assemblies in the reactor core are always low enough and the flow will be stable for the given inlet conditions. However, the peak fuel temperature in number 21 follower fuel assembly decreases from 421 K to 389 K, which is much higher than the temperatures in the standard fuel assembly which may lead to nucleate boiling. More accurate calculations are needed to confirm whether nucleate boiling occurs in number 21 follower fuel assembly.Figure 16
Variations of the peak temperatures in number 12 and number 21 fuel assemblies for various inlet conditions.The exact model was then used to analyze the safety of number 21 follower fuel assembly. The inlet coolant velocity in number 21 follower fuel assembly was obtained from the predictions in Section3.1. Figure 17 shows the variations of the peak cladding temperature and the local saturation temperature for various inlet conditions. As the inlet coolant flow rate increases, the peak cladding temperature gradually decreases. However, the increased coolant velocity also increases the pressure drop which reduces the local saturation temperature downstream, but the decreases are quite small with the saturation temperature always around 385 K. The flow resistances in the channels are relatively small, so the pressure drop does not increase much with increasing velocity.Figure 17
Peak cladding temperature and the local saturation temperature in number 21 follower fuel assembly for various inlet conditions.The curves in Figure17 show that as long as the inlet mass flow rate is greater than 450 kg/s the peak cladding temperature of the hottest follower fuel assembly will be lower than the local saturation temperature. The safety margin then increases as the inlet mass flow rate increases. Therefore, the inlet coolant flow rate should exceed 450 kg/s to ensure the safety of all the fuel assemblies in the HFRR.


## 3.1. Flow Distribution

Figure10 shows the flow distribution inside each assembly at an inlet mass flow rate of 300 kg/s in accordance with the order in Figure 4(a). No coolant flows through the core center due to the inverse flux trap. The normalized flow distribution factors of the standard and follower fuel assemblies are marked in Figure 10. Standard fuel assemblies 8 and 11 have the highest coolant mass flow rates with normalized flow distribution factors of 1.0267. The normalized flow distribution factors for the standard fuel assemblies near the follower fuel assemblies are around 0.985, almost 4% lower than the peak flow rate. The flow rates in the standard fuel assemblies vary by less than 3% from the average. Since the flow resistance coefficients in the follower fuel assemblies are much larger than those in the standard fuel assemblies, the coolant flow rates in the follower fuel assemblies are much lower with smaller variations. Therefore, the predictions indicate that the flow rates inside the fuel assemblies are quite uniform, which means that the inverse flux trap structure has little impact on the flow distribution.Figure 10
Flow distribution inside each assembly for an inlet mass flow rate of 300 kg/s.Figure11 shows the variations of the flow distribution factor in each standard fuel assembly for various inlet conditions. The flow field variations are similar but gradually increase with increasing inlet coolant mass flow rate. At the highest inlet mass flow rate of 450 kg/s, the largest difference between the flow rate in each assembly and the average flow rate is less than 5%. Thus, the results show that the flow distribution in the standard fuel assemblies is almost uniform and there will be sufficient coolant flow within each assembly.Figure 11
Variation of the flow distribution factor in each standard fuel assembly for various inlet flow rates.Figure12 shows the flow fields at the core entrance and outlet at an inlet mass flow rate of 300 kg/s. Since the flow resistance coefficients in the follower fuel assemblies are much larger than those in the standard fuel assemblies, more coolant enters the standard fuel assemblies with velocities 2 or 3 times those in the follower fuel assemblies. A large recirculation zone develops at the fuel assembly outlet with many vortices at the bottom of the core due to the inverse flux trap that blocks the flow in the middle. The coolant flowing out the assemblies is then fully mixed and the velocity is nearly homogeneous there. Thus, the inverse flux trap promotes mixing of the fluid streams exiting the core.Figure 12
Flow fields at the core entrance and outlet for an inlet mass flow rate of 300 kg/s.


## 3.2. Temperature Distribution

The temperature distribution in the core at an inlet mass flow rate of 400 kg/s is shown in Figure13 with only the solid temperature presented in the active core region and the coolant temperature presented in the upper and lower regions. Follower fuel assemblies 21 and 22 have the highest temperature of 395 K, with the highest coolant temperature rise through those assemblies, 12 K. The highest temperatures in the standard fuel assemblies occur at numbers 7 and 12 standard fuel assemblies. Their temperature is about 375 K and the temperature rise of the coolant through these assemblies is about 12 K. As shown in Section 2.2.3, numbers 21 and 22 follower fuel assemblies have the highest average power levels among the follower fuel assemblies. At the same time, the coolant flow rates through these two assemblies are relatively low, as mentioned in Section 3.1. Therefore, numbers 21 and 22 follower fuel assemblies are hottest due to the reduced cooling.Figure 13
Temperature distribution in the core at an inlet mass flow rate of 400 kg/s.Further studies were conducted to ensure the safety of the hottest standard and follower fuel assemblies. Fuel assemblies 7 and 22 are symmetric, and therefore only one was analyzed. The coolant inlet velocities to follower fuel assemblies 12 and 21 were based on the distribution results in Section3.1. Then, the exact model was used to calculate more detailed temperature distributions in these two assemblies.The standard assembly temperatures predicted by the porous media model and the exact model are shown in Figure14 with the follower fuel assembly temperatures shown in Figure 15. The two temperature profiles predicted by the porous media in each figure represent the average solid temperature and the average fluid temperature variations along the axis obtained by the nonequilibrium model, while the corresponding curves of the exact model represent the fuel temperature and the coolant temperature in the channel.Figure 14
Calculated temperatures in the number 12 standard fuel assembly.Figure 15
Calculated temperatures in number 21 follower fuel assembly.Figures14 and 15 show that the temperatures predicted by the porous media model agree well with those predicted by the exact model. The temperature variations predicted by the porous media model are somewhat flatter than those predicted by the exact model because the exact model exactly simulates the axial power distribution while the porous media model is an approximate method. In the porous media model, the heat transfer coefficient, hfs, between the fluid and the solid in the energy conservation equation was specified by the empirical Dittus-Boelter equation [19]. Moreover, hfs and the solid thermal conductivity were assumed to be constant along the whole channel. However, at the inlet and the outlet of the exact model, the flow field in the narrow rectangular channel is not fully developed compared to that modeled by the porous media method. This difference leads to less turbulence and decreases the heat transfer, so the power profiles show different shape at the inlet and outlet region. Nevertheless, these approximations led to small differences between the two models with a maximum difference of less than 1.4%, which indicates that the porous media model can reasonably predict the temperature distribution with the thermal nonequilibrium model. In addition, for simulations of the full reactor core, the porous media model has the additional advantage of significantly reducing the computational cost.Figure16 shows the variations of the peak temperatures in fuel assemblies 12 and 21 predicted by the exact model for various inlet conditions. As the inlet coolant mass flow rate increases from 300 kg/s to 500 kg/s, the peak temperature in number 12 standard fuel assembly gradually decreases from 350 K to 339 K, with all the temperatures lower than the local saturation temperature. Therefore, the predictions indicate that all the temperatures in the standard fuel assemblies in the reactor core are always low enough and the flow will be stable for the given inlet conditions. However, the peak fuel temperature in number 21 follower fuel assembly decreases from 421 K to 389 K, which is much higher than the temperatures in the standard fuel assembly which may lead to nucleate boiling. More accurate calculations are needed to confirm whether nucleate boiling occurs in number 21 follower fuel assembly.Figure 16
Variations of the peak temperatures in number 12 and number 21 fuel assemblies for various inlet conditions.The exact model was then used to analyze the safety of number 21 follower fuel assembly. The inlet coolant velocity in number 21 follower fuel assembly was obtained from the predictions in Section3.1. Figure 17 shows the variations of the peak cladding temperature and the local saturation temperature for various inlet conditions. As the inlet coolant flow rate increases, the peak cladding temperature gradually decreases. However, the increased coolant velocity also increases the pressure drop which reduces the local saturation temperature downstream, but the decreases are quite small with the saturation temperature always around 385 K. The flow resistances in the channels are relatively small, so the pressure drop does not increase much with increasing velocity.Figure 17
Peak cladding temperature and the local saturation temperature in number 21 follower fuel assembly for various inlet conditions.The curves in Figure17 show that as long as the inlet mass flow rate is greater than 450 kg/s the peak cladding temperature of the hottest follower fuel assembly will be lower than the local saturation temperature. The safety margin then increases as the inlet mass flow rate increases. Therefore, the inlet coolant flow rate should exceed 450 kg/s to ensure the safety of all the fuel assemblies in the HFRR.


## 4. Conclusions

A safety analysis of the full HFRR core was conducted using a CFD model. The porous media model was used to simulate the fuel region in the core with the thermal nonequilibrium model used for the energy equation The model predicted the flow and temperature distributions in the entire core.The results show that the flow distribution predicted by the porous media model compares well with experimental data. The thermal nonequilibrium model then provides reasonable predictions of the very different fluid and solid temperature distributions at substantially lower cost than a full Navier-Stokes model.The flow distribution is almost uniform across the entire core with slightly lower coolant flow rates in the follower fuel assemblies than in the standard fuel assemblies. The results show that all the flow channels have sufficient coolant flow to prevent nucleate boiling. Thus, the safety analysis result satisfies the thermal-hydraulic design criteria. The calculations also show that the inverse flux trap increases the mixing and the heat transfer between the flows from all the assembly channels.The predicted temperature distributions show that the inlet coolant flow rate must be greater than 450 kg/s to ensure that boiling will not occur in any of the HFRR fuel assemblies with the peak cladding temperature of the hottest fuel assembly always lower than the local saturation temperature with no nucleate boiling.These calculations are helpful for the thermal-hydraulic design and nuclear safety analysis of the HFRR.



---

*Source: 7152730-2017-12-11.xml*