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

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.


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 trapconverter 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 trapconverter 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 trapconverter 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 thermalhydraulic 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][11][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].

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 U 3 Si 2 -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 × 10 14 n/cm −2 s. 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.  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 Section 2.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.

Description to
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 Figure 3. 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].

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 where is the momentum source term in the momentum equations along the axis ( , , or ), is the density, is the viscosity, V is the velocity, and and 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 where is the permeability and 2 is the inertial loss coefficient. The two coefficients in (2) are the diagonal elements of and 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 resistancevelocity 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: Follower fuel assembly: 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: where Δ is the medium thickness along the axis. As Figure 1(b) shows, the flow channels in the platetype 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 The energy conservation equation for the solid is where f is the fluid temperature, s is the solid temperature, ℎ fs is the volumetric heat transfer coefficient between the fluid and solid, fs 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 h s and h f 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 h , ℎ fs , and fs have to be specified before the calculation. In this study, s h was specified by fitting the axial energy density distribution based on previous model results [2]. Then, s h 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 10 4 . Therefore, ℎ fs was specified by where the diameter and fs were specified based on the specific geometry of the standard and follower fuel assemblies.

Boundary Conditions and Model
Options. Table 1 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.
The initial conditions were based on the specific thermalhydraulic 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. Figure 4 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.

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. Figure 5 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 m 3 /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.

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.  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 − turbulence model, a pressure outlet, and an inlet mass flow rate of 300 kg/s. The normalized flow distribution factor was expressed as where is the mass flow rate in each assembly. The predictions shown in Figure 7 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.
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.  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.

Validation of the Thermal Nonequilibrium
The peak fuel temperatures in Table 2 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. Figure 9(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.

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. Figure 10 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 11 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 12 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.

Temperature Distribution.
The temperature distribution in the core at an inlet mass flow rate of 400 kg/s is shown in Figure 13 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.
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 Section 3.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 Figure 14 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. Figures 14 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, ℎ fs , between the fluid and the solid in the energy conservation equation was specified by the empirical Dittus-Boelter equation [19]. Moreover, ℎ fs 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. Figure 16 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.
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 Section 3.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,   Temperature  390  385  380  375  370  365  360  355  350  345  340  335  330  325  320 315 310 Figure 13: Temperature distribution in the core at an inlet mass flow rate of 400 kg/s. so the pressure drop does not increase much with increasing velocity. The curves in Figure 17 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.

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.

Additional Points
Highlight. (i) The flow and temperature fields in a High Flux Reactor core are predicted using a CFD model. (ii) The porous media approach can accurately simulate the HFRR core. (iii) The porous media energy equation uses the thermal nonequilibrium model.

Conflicts of Interest
The authors declare that they have no conflicts of interest.