FIRST STEPS TO COUPLED HYDRAULIC AND MECHANICAL CALCULATIONS WITHIN A PARAMETER STUDY TO DEFINE POSSIBLE CORE DESIGNS FOR THE CONVERSION OF FRM II

The Forschungs-Neutronenquelle Heinz Meier-Leibnitz (FRM II) is actively participating in the worldwide efforts on developing high-density uranium fuels in order to reduce the enrichment of fuels used in high flux research reactors. This work is part of a parameter study to define possible compatible FRM II core designs for conversion. As a first step, a code-to-code verification is performed and experimental data is used for validation. The Gambill experiment was performed in the early 1960’s in support of the HFIR program and provides results regarding the heat transfer coefficient and friction factors of water flowing through an electrically heated thin rectangular channel. A comparison is made between the Gambill Test and the results simulated by Ansys CFX and STAR-CCM+.


INTRODUCTION
The Forschungs-Neutronenquelle Heinz Maier-Leibnitz (FRM II) is located in Garching and it is Germany's most powerful neutron source. Operating with a thermal power of 20 MW and an undisturbed thermal neutron flux of 8 × 10 14 1 scm 2 , FRM II is the neutron source with the highest flux-to-power ratio worldwide. In order to reduce the use of highly-enriched uranium (HEU, 235 U/U > 20 wt %) in the civil cycle, FRM II takes an active part in the international efforts to develop a high density uranium fuel. With these new fuels, the enrichment used in high performance research reactors can be significantly reduced. To enable efficient cooling of such a high density fuel, special fuel geometries are needed which demands for a broad assessment based on the combination of neutronics, thermohydraulics and mechanics. This work is part of a parameter study that aims to identify possible FRM II core designs suitable for conversion. The focus is set on creating a coupled code system model accounting for neutronics, thermohydraulics and mechanics. The numerical simulations can be used in order to study the influence of various parameters that could be changed for conversion. Several Computational Fluid Dynamics (CFD) and mechanical codes are available as potential candidates to perform a code-to-code verification based on experimental results. FRM II is a high-flux research reactor and as such, it has significantly higher power densities than power reactors. Hence, validation work that has been done for power reactors cannot be directly used for research reactors. In order to perform the validation, experimental data obtained in thermohydraulic conditions close to a high-flux research reactor are required. However, availability of applicable experimental data for comparison is limited. One of relevant experiments that could be used for validation is an experiment performed by W. R. Gambill and R. D. Bundy in support of the High Flux Isotope Reactor (HFIR) program in the early 1960's [1]. It involves measurements of friction factors as well as heat transfer coefficients for forced convection of water flowing through an electrically heated rectangular channel. This experiment is for brevity referred throughout this paper as Gambill Test or Gambill Experiment. The Involute Working Group (IWG) is an alliance between Argonne National Laboratory (ANL), Oak Ridge National Laboratory (ORNL), Institute Laue-Langevin (ILL) and Technical University of Munich (TUM). IWG is proceeding towards the validation of CFD solvers and techniques for high performance research reactors safety analysis. As a first step, Gambill Test in analyzed using Ansys CFX and STAR-CCM+ commercial solvers. A comparison of the results is performed between the two codes and the experimental data.

GAMBILL EXPERIMENT
The Gambill experiment was performed in the early 1960s by W. R. Gambill and R. D. Bundy, as part of the High Flux Isotope Reactor (HFIR) program [2]. The shape of the cooling channel between the fuel plates of HFIR was of importance during the design, since it could influence the friction factors, heat transfer coefficients, and consequently the thermal behavior of the core. A study performed by Levy et al. [3] in 1959 suggested that the heat transfer coefficient for a forced convection in a rectangular channel is 15-30% smaller than that of a round tube and 32-57% lower than Sieder-Tate correlation would predict (see chapter 3). This outcome was surprising and of main concern to Gambill and Bundy, since the resulting heat transfer coefficients were not as high as expected.
Since HFIR would operate with high heat flux, these uncertainties were chosen to be dealt with by performing an experiment representative of HFIR conditions and geometry. One of the key goals of the experiment was to measure the heat transfer coefficients and to examine what correlation would be the best suited for the HFIR thermohydraulic safety analysis. The experimental setup consisted of water circulating at high velocity in an electrically heated thin rectangular channel. The light water coolant was pumped from a 1150 liters tank vertically into the heated test section. Twelve thermocouples were used, each positioned approximately at 5 cm (2 inch) axial intervals on alternately opposing sides of the test section. The coolant pressure was measured by gauges calibrated during steady-state operation. The inlet temperature was regulated by a steam-heated coil, and the test section was thermally insulated.
One of expected problems was an uneven distribution of the heat flux in the peripheral area of the channel. In order to prevent the occurrence of a peak heat flux on the section edges, it was decided that the thickness of the aluminum walls should not be maintained constant. An analysis made for a rectangular channel revealed the necessary thickness of the curved wall required to maintain a uniform heat flux distribution in the periphery. The final dimensions of the cross-section were as shown in Fig. 1.
Many experimental runs were performed for several test conditions, but only the results of run 8, test 7 are explained and analyzed in detail in [1]. Therefore this run is taken as reference for this work. In this test, the saturation temperature was taken from the steam tables [4] to be consistent with the measured terminal static pressure. The bulk temperature t b was assumed to be increasing linearly along the channel length from 69 • C at the inlet to 154 • C at the outlet. A rough calculation showed negligible errors for this approximation. The thermocouples at each axial position measured the temperature at the outer wall of the channel. The wall temperatures were derived by first calculating the temperature drop through the wall ∆t w with: where φ is the heat flux in the coolant, k is the thermal conductivity of the aluminum and α stands for the thickness of the aluminum channel. The wall temperature increased from 111 • C near the inlet to 209 • C near the outlet. The wall roughness was 1.6 µm. The heat transfer coefficient was calculated with the following standard formula: where x is the distance from the heated inlet. The resulting curve was integrated in order to obtain an average h. A plot summarizing run 8, test 7 is shown in Fig. 10.

EMPIRICAL CORRELATIONS
In order to describe the heat transfer process and different fluid flow conditions in fluid dynamics, dimensionless quantities such as Reynolds Re, Prandlt P r and Nusselt number N u are used.
In most broad terms, Re is defined as a ratio between inertia and viscous forces and it is a measure of dynamic similarity between fluid flows. Re is calculated as: where v is the fluid velocity, d the hydraulic diameter and ν the kinematic viscosity. Physics of Reactors Transition to a Scalable Nuclear Future P r is used to describe the thermal properties in a fluid by the ratio of momentum diffusivity ν to thermal diffusivity a: For heat transfer, P r is used for estimating the ratio of the thicknesses of the velocity and thermal boundary layer. P r values for FRM II and HFIR are 4.56 [5] and 2.16 [6], respectively.
N u is often referred to as a dimensionless heat transfer coefficient and it is defined as the ratio of convective to conductive heat transferred in the direction normal to the boundary: where h represents the heat transfer coefficient, d is the hydraulic diameter and k the thermal conductivity of the fluid.
Heat transfer coefficients are not straight-forward to measure. Hence, they are often deduced by applying empirical correlations between the Nusselt number, Reynolds number and Prandtl number. For forced convection: Empirical relations for the Nusselt number have been developed for different geometries and different flow conditions such as the Sieder-Tate [7] and Hausen [8] correlations. Gambill and Bundy compared their experimental Nusselt number to the ones obtained with these correlations.

Sieder-Tate Correlation
Sieder-Tate correlation [7] represents an implicit function for turbulent flows which accounts for the viscosity gradient between the wall and bulk. In contrast to the previous studies performed on this field, Sieder and Tate recognized that a large viscosity gradient is found not only on the axial direction, but on the radial as well. Therefore, the calculated Nusselt number using this method is more adequate than N u from previous correlations because of the incorporation of the viscosity correction factor µ b /µ w , where µ b and µ w are referred to as the viscosity at bulk liquid temperature and at wall temperature. It is valid for Re ≥ 10000 and 0.7 < P r < 16700. The calculated Nusselt number is expressed as: The subscript bm stands for the average bulk.

Hausen Correlation
Hausen modified the Sieder-Tate correlation to be used for both transition and fully developed turbulent flows [8]. The Nusselt number can be expressed as: where d is the hydraulic diameter and L h the heated length.
All property values are considered at t b , except µ w . Eq. 8 provides the average heat-transfer coefficient between the starting point of the heated section and the distance from the inlet.

MODELING OF TURBULENCE AND NEAR-WALL REGION WITH CFD
As mentioned before, the flow in the Gambill experiment is turbulent, with Re ≈ 2.1 × 10 5 . The near-wall region can be divided into three layers: the viscous sublayer, located close to the wall where viscosity plays a dominant role, the logarithmic layer where the turbulence is superior and the buffer layer, located between them, where turbulence and viscosity have the same order of magnitude of influence [9]. One important parameter when it comes to mesh description and quality is the dimensionless wall distance y+ defined as: where y is the absolute distance from the wall, u τ is the friction velocity and ν is the kinematic viscosity. y+ can be interpreted as a local Reynolds number since its value gives information regarding the importance of the viscous and turbulent processes. The viscous sublayer lies in a region within y+ < 5 also called the linear region, the buffer layer within 5 < y+ < 30 and the logarithmic layer at y+ > 30 [9].
As a practical approach to save computational resources, the so-called wall functions are often applied near the wall. The wall functions are empirical formulations that establish reasonable conditions near the wall in order to resolve the viscous layer without the need of fine meshes.
In certain flow conditions, the presence of rough walls can affect the flow considerably as they tend to increase the turbulence near the wall, and therefore increase the heat transfer coefficients. However, in highly turbulent flows, roughness does not play a notably important role.

CFD MODEL
In this section, the CFD turbulence models used in Ansys CFX and STAR-CCM+ are described and the differences between them are pointed out.

Ansys CFX
A short overview of the turbulence models used in Ansys CFX as well as wall functions is given in this subsection. Furthermore, the CFD model built in Ansys CFX is explained in detail.
As mentioned in section 4, wall functions are empirical formulas used to model the near wall region without resolving the viscous sublayer. Automatic and scalable are the two common wall functions implemented in Ansys CFX [10].
Scalable wall functions can be applied to arbitrarily fine meshes and allow the user to perform a consistent mesh refinement independent of Re. The idea of this implementation is to put a lower limit for y+, in order to ensure that all mesh points are not on the viscous sublayer, reducing the chance of numerical errors.
Automatic wall functions have the ability of switching between low-Re formulation to wall function treatment by using a blending function to tune the behavior of the viscous sublayer and logarithmic region. More information regarding wall functions in Ansys CFX can be found in [10].
Regarding roughness, Ansys CFX typically uses the sand-grain roughness to model rough walls. The sand-grain roughness averages the technical roughness, and it is incorporated into the corresponding equations, instead of the absolute roughness [10]. The logarithmic boundary layer velocity profile is shifted upon rough wall usage. Scalable wall functions are recommended for rough walls in Ansys CFX according to [10].
Three different turbulence models were used in the Ansys CFX simulations for this paper: k-, SST and Reynolds Stress model (RS).
k-is a 2-equation model, where k is the turbulent kinetic energy and determines the length scale of the turbulence. This model can struggle to model the region next to a heated wall, unless a fine mesh with y+ < 0.2 is provided. It is a widely used model since it offers a reasonable compromise between accuracy and robustness. For y+>0.2, Ansys CFX uses scalable wall functions.
The Wilcox (standard) k-ω is a 2-equation model very similar to k-except that it determines the time scale of turbulence ω instead of . The strong aspect of this model is the ability to model the near wall region. The convergence behavior, the difficulty to model the free stream flow and high sensitivity to initial conditions are the disadvantages of k-ω [10].
SST model was developed to overcome some of the difficulties of k-ω. It combines the advantages of k-and k-ω models. It uses k-ω approach in order to model the near-wall region and k-to model the free stream. The transition between the two models is accomplished by using blending functions [10]. Automatic wall functions are used with SST in Ansys CFX for meshes with y+>2.
In RS model, all the components of the Reynolds stress tensor are calculated, requiring slightly more computational effort than k-and SST. Many RS models are available in Ansys CFX, however, SSG Reynolds Stress model with scalable wall functions was used here. SSG has the advantage to describe more accurately the distribution of the Reynolds stresses and it is the recommended RS model for most flows by Ansys CFX [10].
Based on the information available in Gambill and Bundy report [1], a CFD model of the experimental section has been generated. The geometry for the model was built in SOLIDWORKS. An unheated extension of 10 cm was attached to the inlet allowing the flow to develop before the water reaches the heated section. The meshes that were used in Ansys CFX were generated with ICEM-CFD, a meshing tool provided by ANSYS suite. Two meshes were built with different mesh densities by varying the y+ values (see Table. 1 for details). The difference between mesh 1a with 30<y+<130 and mesh 2a with 0.2<y+<2 stands on the cell density near the liquid-solid interface. For instance, mesh 1a has fewer cells in the liquid-solid interface, while mesh 2a has notably more cells. The meshes are similar in the middle of the channel and simulations were run with both meshes.
The boundary conditions for run 8 of test 7 included an inlet velocity v 1 = 15.7m/s, inlet temperature t 1 = 69 • C, outlet pressure p 2 = 27.9bar and wall roughness ρ = 1.6µm. Since no specific information was mentioned regarding the heat applied to the aluminum channel in the Gambill experiment report, a value for the heat source was roughly estimated by using the data in Fig. 10 with: where ∆T is the difference between outlet and inlet temperature, here 82 • C, c p and m are the specific heat capacity at constant pressure and mass of water, respectively. The wall temperature t w and the heat flux φ x were calculated at a point in the solid-liquid interface. The bulk temperature was determined as a weighted average of the temperature with the mass flow in a cross section inside the tube. Both t w and t b were assessed at the same axial positions. Adiabatic conditions were used for the outer walls of the channel, as well as for the inlet and outlet faces of the aluminum channel.
The aluminum is heated by an electric current passing through the channel. Therefore, electrical heating in the simulations was applied as a Joule heating, whereas the volumetric heat depends on aluminum electrical resistivity. Since the electrical resistivity increases with temperature, the heat generated in the metal increases as well. Therefore, the applied heat was iteratively increased until the output heat matches the value from Fig. 10, and consequently the water temperature at outlet matched the one from the Gambill report. A heating voltage of approximately 9.7 V was applied to the aluminum channel and was kept constant in all the simulations. Therefore, the applied heat altered from one simulation to the other depending on which turbulence model was used. The resistivity data was taken from [11].
The IAPWS-IF97 material database was used for water properties. The simulations were run for smooth and rough channel walls to observe the roughness effect on the heat transfer parameters. The roughness could only be applied to mesh 1a, since the first cell distance from the wall (y) has to be at least twice as large as the roughness, which was not the case for mesh 2a.

STAR-CCM+
This subsection contains a brief description of the CFD models developed at Argonne National Laboratory in a parallel study under IWG activities. Many more details on these models and accompanying results can be found in the open report [12].
While the general CFD modeling techniques apply equally to Ansys CFX and STAR-CCM+, the specific implementations of turbulence models vary between these solvers. The general description of these models is available in the software documentation [13]. In this work out of available turbulence models in STAR-CCM+ the following three were used: • Realizable K-Epsilon two layer turbulence model with two layer all y+ wall treatment • SST (Menter) K-Omega with all y+ wall treatment • Reynolds Stress Turbulence with all y+ wall treatment All y+ wall treatment is a hybrid treatment that blends the low y+ wall treatment for fine meshes, and the high +y wall treatment for coarse meshes [13]. Low y+ wall treatment resolves the viscous sublayer and few assumptions are needed to predict the flow near the wall boundary. It means that the governing transport equations are solved all the way to the wall. In high y+ treatment, the viscous sublayer is not resolved directly. As a substitute, wall functions that assume flow profiles are used to obtain boundary conditions for the continuum equations. In a similar way as in Ansys CFX implementation, these wall functions can contain wall roughness definitions. As a requirement for that model to work, the centroid of the first near wall cell must be located further away from the wall than the roughness scale.
The geometry of the specimen and the artificial inlet extension were built in the STAR-CCM+ internal CAD modeler. The extension was about 11.5 cm long what is equivalent to about 43 hydraulic diameters of the channel. Usually 40 hydraulic diameters are more than enough for a turbulent flow to fully develop. In the experiment, the actual flow profile at the inlet to the specimen was unknown but fully developed profile was assumed to be most appropriate. For both CFD tools used, Ansys CFX and STAR-CCM+, also the thermal boundary layer was not developed in the same way as in the experiment. In the simulations adiabatic conditions were used on the inlet face of the aluminum specimen. However, in the experiment solid metal components were in contact with the specimen. These components were conducting heat away from the heated specimen and thus the calculated temperature in the specimen near the inlet is expected to be different than in the experiment. However, this should not affect the validity of the models in general.
As with the model built in Ansys CFX, two different meshes were used in STAR-CCM+ denoted here as mesh 1b and mesh 2b. This distinction refers to the resolution of the mesh near the solidfluid interface. The mesh with high y+ above 30 contained in total about 3.3 million of cells while the model with fine mesh, with y+ in the range between 2.0 and 5.0, contained about 5.1 million of elements. Linear velocity profile is assumed for meshes with y+<5. In the mesh regime with y+ between 5 and 30, neither linear nor logarithmic velocity profiles are accurate and the numerical errors may be significant. For that reason, such meshes are in general not recommended.    2 shows a cross section through the model with coarse mesh (1b). Fig. 3 shows details of fine mesh (2b) and coarse mesh (1b) near the solid-fluid interface. It can be noted that the meshes are almost identical in the middle of the channel and the difference is only visible in the resolution near the interface wall.
The properties of water were defined as temperature dependent polynomials of different orders based on the data obtained from NIST Chemistry WebBook [14] for water at pressure of 30 bars.
Heat generation within the specimen can be simulated in STAR-CCM+ with the use of several methods, including: total heat source, volumetric heat source and Joule heating. It was decided that the volumetric heat source is the most appropriate method for the purpose of this study. The volumetric heat source can be defined to be dependent on resistivity of the aluminum which varies with the temperature. The temperature dependent resistivity of aluminum was defined in a similar fashion as in the Ansys CFX model based on [11]. Adiabatic boundary conditions were assumed on the external boundaries of the aluminum specimen. The water temperature at the inlet was set to be the same as in the experiment at 68.5 • C. The power input was initially adjusted so the bulk water temperature at the outlet of the base model closely matches the experimental value of 154.5 • C. The total power input was subsequently kept constant in all simulations as opposed to the constant voltage in Ansys CFX. As a result, the outlet bulk temperature was constant in all simulations as well.

RESULTS
In this section the results from the Gambill experiment are summarized and compared with the CFD simulations.

Gambill Experiment Results
In the Gambill test, the heat transfer parameters were compared to the Sieder-Tate and Hausen correlations. Regarding Sieder-Tate correlation, it was decided that Eq. 7 with a coefficient of 0.024 represents a reasonable design correlation for HFIR program.
However, slightly better result was provided by the Hausen equation. Around 240 local heat transfer coefficients were derived from the measurements. These values are known to scatter significantly since their measurements are relied on point measurements. The local h by the Hausen correlation were calculated by using Eq. 8. At the end, a coefficient of 0.105 in Eq. 8 was considered a reasonable design correlation for HFIR steady-state analysis.

CFD Simulation Results
The wall temperatures, pressure drop and comparison to empirical correlations are presented and explained in this subsection.
The temperature profiles of the wall resulting from the Ansys CFX simulations are compared to the experimental results. The wall temperatures obtained from all three turbulence models for both meshes ( Fig. 4(a)) match well the experimental data. Slight differences are seen in t w in the region close to inlet, especially in the results calculated by RS model in mesh 2a. This could be related to the non-constant heat deposition in the aluminum channel. Nevertheless, all turbulence models calculate well matching t w at the end of the channel.
In Fig. 4(b), t w for smooth and rough walls are compared. Slightly smaller wall temperatures are expected as a result of channel roughness. For this reason, smooth walls can be assumed to be more conservative. The maximum t w difference between smooth and rough walls is approximately 10 • C for all three turbulence models.
The pressure drop ∆p through the channel is shown in Fig. 5 for the performed Ansys CFX simulations. The outlet pressure has the same value for all the simulations since this was put as an initial boundary condition. Mesh 1a and mesh 2a provide substantially the same pressure at the inlet with a maximum deviation of 0.4 bar. However, ∆p is larger by 1bar for all turbulence models used in  conjunction with the rough wall model. The friction factors are larger and therefore, larger is the pressure drop.
The results of the Ansys CFX and STAR-CCM+ simulations compared to the experimental data and Sieder-Tate and Hausen correlations are depicted in Fig. 6-8, representing the used meshes and roughness, respectively. In the simulations performed with STAR-CCM+ for fine mesh (2b) and k-omega model, problems with convergence were encountered and for that reason these results are not included in this paper. When simulations are performed with mesh 2a, the results of the Hausen equation (Eq. 8) from both Ansys CFX and STAR-CCM+ are very similar. Due to a fine mesh in the liquid-solid interface, good results were achieved using k-and RS with less than 10%      The results for mesh 1a with wall roughness are presented in Fig. 8. The N u according to Hausen and Sieder-Tate equations exhibits slightly higher values in the Ansys CFX simulations than in the results obtained by STAR-CCM+, in which the N u values are closer to the empirical correlations. The observed differences between Ansys CFX and STAR-CCM+ may come as a consequence of different approaches in modeling the roughness as well as from slight differences in the total power input as a consequence of keeping voltage constant in the case of Ansys CFX or keeping the power constant in STAR-CCM+. In order to be able to compare the models where roughness is included with proper empirical correlations, Gnielinski equation [15] will be used for future calculations.
Nevertheless, the roughness is supposed to play a minor role for forced convection in such thin channels, like present in the Gambill Test.
The heat transfer coefficients were calculated for all simulation conditions by using Eq. 2 and the results are plotted in Fig. 9. Despite the fact that the heat transfer coefficients as determined in the Gambill report are assumed to be not that accurate considering their point calculation, the CFD simulations obtain very similar result with less than 10% deviation with the experiment in the end of the channel. However, further investigation is required to gain a deeper understanding of the behavior of the heat transfer coefficients in the Gambill Test as well as in the CFD simulations.

CONCLUSIONS AND OUTLOOK
One test condition of the Gambill Test was simulated by using Ansys CFX and STAR-CCM+. Two meshes were generated with different mesh densities at the liquid-solid interface. Additionally, a comparison was made between the simulation run with the two different meshes as well as with different turbulence models, such as k-, SST and RS model. The heat transfer coefficients were compared to the Gambill experimental results and a good agreement was observed, especially for the end of the channel. Simulations with constant heat deposition will be run in the future with Ansys CFX. Additionally, an investigation of the heat transfer process will take place with Abaqus CFD. After performing this code validation for Ansys CFX and STAR-CCM+, the Involute Working Group (IWG) continues its way of validating CFD calculations for high-performance research reactors.