Benchmark Experiments for Steady-State Natural Convection in Fuel Rod Bundles

The purpose of this work is to provide experimental benchmark validation data of steady-state, single-phase, natural convection flow through a vertical nuclear fuel rod bundle for three-dimensional computational fluid dynamics (CFD) simulations. The experimental facility, instrumentation, data acquisition system, BCs, and system response quantities (SRQs) are described in detail. All inflow, BCs, and SRQs are available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection.

The files attached to this article are also available on the Utah State University Library's Digital Commons online database corresponding to Ref. [1]. The geometry, BCs, and SRQs are available under the same file names presented in this work. Three other benchmark cases (one forced, unheated case, and two mixed convection cases with the same heat flux levels presented here) are also available for download to be used as preliminary simulations for CFD modelers. We anticipate that these cases are less challenging for CFD compared to natural convection.

The BCs and SRQs are table formatted as comma-separated values (*.csv) text files. The geometry needed to perform the simulations is included in three, widely used formats: Parasolid (*.x_t), STEP (*.stp), and Binary Stereolithography (*.stl). The files for the geometry may be downloaded in the file Geometry.zip, which is available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection. The BCs and SRQs for the two heat flux cases considered, 400 W/m² and 700 W/m² (referred to as Natural400 and Natural700, respectively), in this study may be downloaded in the Natural400.zip and Natural700.zip files, respectively, which are available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection. In each of these files, the name of the variable is presented at the column header with the units in brackets, and the uncertainty for that variable is indicated by a lowercase u preceding the variable name and the standard deviation by a lowercase s. For example, the column header “W (m/s)” is the W-velocity in units of meter per second and the column header “uW (m/s)” is the uncertainty in the W-velocity in units of meters per second, and the column “sW (m/s)” is the standard deviation of the w-velocity component in units of meters per second.

With the exponential increase in the computing power over the last 30 years, research groups rely heavily on the results of simulation for design and decision making processes. In order to ensure that the solutions provided by these simulations are sufficiently accurate, it is necessary that the models be validated to determine their accuracy and uncertainty. The U.S. Department of Defense (DoD) released a recommended practices guide for the Verification, Validation and Accreditation (VV&A) of Models and Simulations (M&S) [2] for more informed judgment and assessment. In 1998, the American Institute for Aeronautics and Astronautics (AIAA) published a standard to be used for the verification and validation of CFD models [3]. The American Society of Mechanical Engineers (ASME) also published a standard for Verification and Validation (V&V) of CFD and heat transfer in 2009 called V&V-20 which has been extended from the AIAA standard to include heat transfer problems [4]. Verification and validation is a critical part of the model development and application in simulation use.

Using these standards, one can ensure that the simulation results are accurate. According to Oberkampf and Roy [5], “Validation is the process of assessing the physical accuracy of a mathematical model based on comparisons between computational results and experimental data.” In contrast to discovery experiments where the experimentalist is attempting to discover or understand more about a certain physical phenomenon, validation experiments simply look to measure the physical nature of a system for the purpose of providing a complete description of the system which can be used to test the validity of a model [6].

Validation of complete, large-scale systems is generally not feasible. Methodology for a hierarchical approach to validation has been outlined in Refs. [7–10]. Figure 1 shows the validation hierarchy for complete systems which are broken down into smaller subsystems for validation purposes. According to this hierarchy, the present work is considered a benchmark-level case due to the highly coupled fluid flow and heat transfer occurring in the flow domain.

When performing a validation experiment, careful documentation of all BCs, initial conditions (ICs) and SRQs, and their associated uncertainties is vital. Examples of boundary conditions may be locations of physical walls/geometries or the temperature distribution of a surface. Initial conditions are similar to boundary conditions but exist at the start of an experiment and change over time, e.g., the temperature response of a body after starting an electric heater. System response quantities are the outputs of the system that are recorded for comparison to the output of a simulation, such as a velocity profile in the fluid domain or wall shear field at a boundary. System response quantities are provided to the analyst after performing the simulations in order to ensure that the data are used for validation rather than the model calibration [5].

Careful experimental design is needed in order to have sufficient insight to plan for SRQs to be measured. Modelers and experimentalists should be involved through all aspects of the design and execution of a validation experiment for optimum success. Knowledge of the relative difficulty of specific SRQs will aid both parties in exercising the limits of the CFD model especially in cases of complex physics. Figure 2 shows the spectrum of SRQs and their relative difficulty to both measure and predict. Using this spectrum in the planning stages of a validation study allows for more rigorous testing of the models applied.

where U_num is the numerical uncertainty of the model, U_input is the simulation uncertainty due to errors in the boundary conditions, and U_D is the experimental SRQ uncertainty. The validation comparison error is satisfactory if the validation comparison error E is sufficiently smaller than the validation uncertainty, U_val, for the intended application of the model. This work will provide U_D for the SRQs used to validate the CFD models as well as uncertainties of the BCs applied to the CFD models in order to determine U_input.

where t_α∕2,v is the quantile of the t distribution for v = n − 1 degrees-of-freedom, and s is the sample standard deviation. Using this validation metric allows for inclusion of multiple validation experiments in determining the comparison error, decreasing the interval by $1/n$ as more experimental data are included.

where the heat transfer q to or from the object is equal to the product of h, the surface area A_s, and the temperature difference between the surface T_s and the surrounding fluid T_∞. The Rayleigh number, Ra = gβ(T_s − T_∞)L³/(να), is useful for determining transition to turbulence which, for a vertical heated plate, is generally assumed to occur at Ra ≈ 10⁹. Here, g is the acceleration due to gravity, β is the coefficient of thermal expansion, L is the distance along the flat plate, ν is the kinematic viscosity, and α is the thermal diffusivity.

where k_f is the thermal conductivity of the fluid. For natural convection, these correlations are a function of Rayleigh number and Prandtl number (Pr = ν/α). The convection coefficient can be determined by using these empirical correlations. In many cases when the flow conditions cannot easily be determined due to complex geometry, the only practical approach is to determine the Nusselt number experimentally [11].

Few validation studies exist for buoyantly driven flows. De Vahl Davis [12], and DeVahl Davis and Jones [13] were among the first to perform benchmark studies for natural convection in an enclosure, although these benchmarks were using CFD simulations on uniform meshes of up to only 81 × 81 cells due to limited computational resources at the time and were compared to analytical solutions for the laminar flow regime.

The first full validation experiment of this type was performed in 1998 by Leong et al. [14] for natural convection in a cubical cavity with opposing faces at different temperatures. The authors report Nusselt numbers accurate to 1% for Rayleigh numbers in the laminar flow regime. The experimental results were compared to CFD for Ra = 4 × 10⁴ and found to be accurate to within 0.3%. Later, Leong, et al. [15] also performed benchmark experiments using the same apparatus for a wider range of Rayleigh numbers up to 10⁸ in order to determine Nusselt numbers for the cavity at different tilt angles.

Other variations of validation experiments for natural convection in cavities were performed later. Mamun et al. [16] provided an extension of Leong's previous work by using a “double inclined” cubical cavity for Rayleigh numbers ranging from 10³ to 3 × 10⁸. Ampofo and Karayiannis [17] performed validation experiments for turbulent natural convection in a square cavity (Ra > 10⁹).

Betts and Bokhari [18] have provided a detailed validation study for turbulent natural convection in an enclosed tall cavity. This study not only provides integral quantities such as the Nusselt number and heat transfer coefficient but also temperature and velocity profiles through the vertical direction of the flow. The experiment presented is very extensive and provides substantial detail for use in CFD validation. Tian and Karayiannis [19] also performed a turbulent experiment for a square cavity providing both velocity and temperature profiles within the cavity for use in validation. These two studies have provided high-quality validation data for natural convection flows. However, they are only applicable to recirculating cavity flow. Natural convection open channel flow can be a much more difficult problem to simulate due to tight coupling of fluid flow with fluid properties where slight variations in boundary conditions can dramatically affect the solution.

Natural convection flows differ from forced convection flows by their BCs. Forced convection is driven by an inflow BC or applied pressure gradient. For natural convection, however, a temperature or heat input BC induces flow due to the effects of buoyancy. Thus, application of an inflow BC to natural convection over-prescribes the solution, allowing the inflow measurement in natural convection problems to be used as an SRQ.

The rotatable buoyancy tunnel (RoBuT) is a wind tunnel facility at Utah State University designed for performing CFD validation benchmark experiments with buoyancy either aiding or opposing the flow, depending on the tunnel orientation. The wind tunnel is shown in Fig. 3 with the major components labeled. The facility test section is a square cross section of 0.305 m × 0.305 m and is 2 m in length. The coordinate system lies at the center leading edge of the opaque wall as shown with the positive z-axis being aligned with the streamwise direction.

The fluid enters through the bottom of the inlet contraction and flows upward. Flow entering the wind tunnel first passes through an array of 1/4 in, Schedule 40 PVC tubes with $∅3.2 mm$ holes spaced 25.4 mm apart through which olive oil vapor is introduced into the flow for particle image velocimetry (PIV) measurements. It then passes through a single row, aluminum fin/copper tube heat exchanger (which was unused for these benchmark experiments but still present in the flow) followed by a honeycomb flow straightener and finally two high porosity screens. The fluid then accelerates uniformly through the inlet contraction, which had an area ratio of 6.25:1 and was 0.914 m in length, before entering the test section. The flow conditioning system is too complex to model, and so it is instead described by the pressure drop that it generates.

The test section is constructed of three transparent polycarbonate walls with several antireflective windows for optical access during fluid measurement and an opaque aluminum wall for structural rigidity. The construction of the test section and fuel rod bundle will be discussed in Sec. 3.

After leaving the test section, the fluid passes through the outlet transformation which is constructed of four sections that are bolted together. The transformation converts the 0.305 m × 0.305 m test section cross section to a circular cross section of $∅ 102 mm$ over a length of 0.686 m. The inlet contraction and outlet transformation are constructed of fiberglass with a highly polished gel coating. A flexible baffle between the outlet transformation and blower was removed in order for the fluid to freely exit the test section for this study. The baffle was left in place for the mixed convection cases presented in Ref. [1].

Two Laskin nozzle seeding canisters were used to introduce seed to the fluid flow for PIV measurements. The particle diameter of the olive oil vapor produced by the seeders is a function of line pressure and nozzle hole diameter [20]. A TSI, Inc. (Shoreview, MN) aerodynamic particle size spectrometer was used to measure the physical particle diameter and was found to have a mean of 1.7 μm. Air was pumped through the body of the seeding canisters in order to dilute the seed exiting the PVC array before being entrained. Holes of $∅3.2 mm$ were used to evenly disperse the seed into the flow. The volume fraction of the seed was found to be ∼10⁻⁹ and should have no impact on the thermodynamic properties of the air.

The test section has a 0.305 m × 0.305 m (12 in × 12 in) internal cross section and is 2 m in length. Three of the four walls are constructed of transparent Lexan™ polycarbonate and are 12.7 mm thick. The fourth wall is constructed of a 12.7 mm thick Aluminum 6061-T6 plate for structural rigidity and has been painted flat black and coated with a rhodamine fluorescent die for filtering reflections during PIV data acquisition. Each wall is divided into four interlocking sections for ease of manufacturing, assembly, and maintenance. The top wall panels (opposite the aluminum wall) are easily removed for internal cleaning as well as for PIV calibration throughout the duration of the experiment. Rectangular optical windows are placed at the inlet of the test section for PIV measurements to reduce the amount of refraction inside the polycarbonate walls. Circular optical windows are placed at the midpoint in the z-direction between each grid spacer for SRQ PIV measurements in the flow-wise direction.

The coordinate system shown in Fig. 3 is located at the center of the leading edge of the opaque, aluminum wall. The left and right walls are positioned at x = 152.4 mm and x = −152.4 mm, respectively, while the top wall is positioned at y = 304.8 mm.

The fuel rod assembly model consisted of four 1.58 m long 6061-T6 aluminum cylinders of $∅31.75 mm$⁠, each divided into four interlocking sections, referred to herein as rod sections. Contained in each rod section is a 20 W split-sheath cartridge heater from Dalton Electric (Part No. W3H126) of $∅9.53 mm$ and 323.85 mm in length with a 314 mm heated length. The heaters were powered by two BK Precision Model 9174 programmable power supplies. A constant power condition was applied to the heaters in order to achieve a desired rod surface heat flux.

Each rod section was constructed of two cylinder halves machined internally to house the cartridge heaters and thermocouples (TCs). The TCs were embedded from the inside near the outer surface of the rods with thermally conductive epoxy (Dow-Corning 3-6751) and routed outside the rods and test section via an aluminum-wire conduit. The two cylinder halves were bolted together using six #8-32 stainless steel screws that were torqued to 2.3 N·m. The counter bore holes were filled with 6061-T6 aluminum rod using Duralco 4540 machinable epoxy.

Oxidation of aluminum causes wide variation in radiative properties [21,22]. The outer surface of the rod sections was plated with electroless nickel to generate low, constant emissivity for reduction of radiation heat transfer. The total hemispherical emissivity of the plating was measured by Optotherm Thermal Imaging and reported as ϵ_s = 0.10 ± 0.01. Any unplated surfaces of the fuel rod assembly were polished. Figures 4 and 5 show the surfaces of the model that were plated and polished. Ceramic spacers were placed between the rod sections and wire conduit to eliminate conductive heat transfer.

A hemispherical aluminum cap was placed on the upstream end of each fuel rod to limit flow separation which results in long time-scale turbulent fluctuations. Similarly, the trailing end of each fuel rod was treated with 128 mm long aluminum cones to maintain attached flow. Fiberglass tape was used to thermally isolate the cones from the heated rods. Aluminum bars were attached to the trailing end of the rods for suspending the fuel rods in the test section.

Steel 4130 grid spacers, designed similar to corrugated box separators, were used to maintain fuel rod spacing in the assembly. Aluminum borders were machined with slots to maintain the grid spacer shape and were embedded in the polycarbonate walls such that the interior walls of the test section were smooth as shown in Fig. 6. The rods were centered in the grid spacers with fiberglass-tipped, #4-40 set screws. Swirl elements were designed as a simple 45–45–90 deg triangle bent at 45 deg from vertical toward the rods as shown in Fig. 7. The swirl elements were placed in the first four grid spacers in the streamwise direction while the final grid spacer did not have swirl elements in order to accommodate the aluminum bars used to suspend the fuel rods in the test section. Two swirl elements are bent toward each rod resulting in a 17.9% blockage ratio. While considerably less effective at mixing than commercially produced mixing vanes, this shape still generates the same physics that are present in an actual fuel bundle.

The pitch-to-diameter ratio of the assembly is 1.52. Both the rod diameter and pitch-to-diameter ratio are larger than a prototypical PWR assembly which was necessary to allow for instrumentation and assembly purposes but still allows the model to capture the physics associated with dry fuel storage.

Various TCs were placed in the facility for controlling and monitoring room conditions. A proportional-integral-derivative (PID) controller used the heating and air conditioning system to maintain the room temperature at 20 ± 1.0° through the duration of the study. The ambient temperature and relative humidity were measured using an Omega HX93A relative humidity and temperature sensor with measurement uncertainties of ±2.5% for humidity and ±0.6 °C for temperature. Ambient pressure was recorded using an SB-100 barometric pressure sensor from Apogee Instruments with a measurement uncertainty of ±1.5%. An NI USB-9215 A 4-channel, ±10 V, 16-bit, Analog Input DAQ was used to measure the voltage from each of these sensors. Samples at 1 Hz were acquired then averaged and recorded once per minute. The ambient temperature, atmospheric pressure, and relative humidity are sufficient to determine air properties. It should be noted that the atmospheric pressure is much lower than one standard atmosphere due to the elevation where the data were acquired (1456 m above sea level). The air properties are included in the AmbientConditions.csv files, which are available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection, for the respective cases.

Fifteen TCs (in a 5 × 3 array in the x- and y-directions, respectively) were placed immediately downstream of the honeycomb flow straightener at the inlet of the wind tunnel (on the z = −0.941 m plane) for measuring the incoming fluid temperature. The temperature distribution at the wind tunnel inlet and test section inlet was assumed to be the same when calculating mass flow rate through the wind tunnel. The inflow temperatures are found in the InflowTemperatures.csv files, which are available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection.

It is not feasible to model the complex geometries found in the flow conditioning portion of the wind tunnel. Therefore, the pressure drop across this portion of the wind tunnel was measured using a 1-Torr Baratron-differential-pressure sensor connected to a MKS207D Signal Conditioner. The voltage output by the signal conditioner was read using an NI-9205 16-channel, ±10 V, 16-bit, Analog Input module. A total of 3000 samples were acquired at 5 Hz and averaged for each test case. The pressure drop for each case is available in the FlowConditionerPressureDrop.csv files, which are available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection.

Preliminary CFD simulations were performed during the design phase of the experiments and found no wall temperature variation across the lateral direction of each wall though there were slight temperature gradients in the streamwise direction. Thus, 21 TCs were placed within 2.5 mm of the inner surface of each wall at three evenly spaced lateral positions and seven streamwise positions, near the inlet, slightly downstream of each grid spacer and near the outlet.

As mentioned previously, TCs were embedded inside each rod section. Preliminary CFD simulations also indicated variation in surface temperature due to the presence of the swirl elements near the rods. It was determined that four TCs would be placed circumferentially at three streamwise positions in each rod section for a total of 48 TCs per fuel rod. Machining of the internals of each rod section allowed the TCs to be placed at 2.5 mm of the rod surface. TCs were aligned with the small gap between adjacent fuel rods and at 90 deg spacing at a given streamwise position.

All four heaters in the rod sections at a given streamwise position were wired in parallel and powered from a single output from the power supply. Measurement of the unpowered heater resistance indicated a variation of 178.1 ± 0.98 Ω, and the current was assumed to be evenly distributed to each heater. The current supplied to the heaters was read from the programmable power supplies which had a readout uncertainty of ±0.1% and bias uncertainty of 1 mA. The voltage across the four heaters at each streamwise position was measured using a voltage divider circuit and an M-series NI PCI-6221 DAQ card. High accuracy, ultra-stable film resistors from Caddock of 10 MΩ and 1 MΩ (Part No. USF240-10.0 M-0.1%-2PPM and USF240-1.00 M-0.1%-2PPM, respectively) made up the voltage divider circuit each with an accuracy of ±0.1%.

Sixteen TCs (in a 4 × 4 array) were placed at the outlet of the test section (z = 2 m) in an evenly spaced square grid pattern. TCs were suspended by four $∅ 1.58 mm$ steel rods and routed outside the test section through the aluminum wall. The TCs were evenly spaced 61 mm apart in both the x- and y-directions, and the array was centered in the test section outlet.

The TCs placed inside each rod section, inside the polycarbonate walls, and at the test section outlet were welded in-house using an Argon-shielded welder. Fiber-glass shielded Type-K TC wire from Omega with special limits of error was used. Each of these TCs was calibrated using an Isotech FASTCAL-M calibrator with an accuracy of 0.3 °C across a range of 30–180 °C. An average calibration curve was applied to these TCs because they were made from the same spool of wire. Outside this range, an uncertainty of 1.1 °C from the manufacturer was used. All other TCs in the facility (including the inlet TCs described previously) were uncalibrated, and the manufacturer's uncertainty of 1.1 °C was used. The cold-junction compensation of the NI-9213 modules had an accuracy of 0.8 °C applied to the uncalibrated TCs.

All fluid velocity measurements were made using a stereoscopic PIV system that provides nonintrusive velocity field measurements. Two Imager sCMOS cameras from LaVision [23] were employed with a dual cavity Quantel Evergreen (frequency-doubled 532 nm Nd:YAG, 100 mJ/pulse, 25 Hz per cavity) and a programmable timing unit. A focusing optics and a cylindrical lens were placed on the front of the laser to convert the laser beam to a sheet with adjustable thickness. The Imager sCMOS cameras had a 16-bit CMOS sensor that was 2560 × 2160 pixels with a pixel size of 6.5 μm with acquisition speeds of up to 50 frames per second.

The cameras were fitted with Nikon AF Nikkor 28 mm f/2.8 D fixed focal length lenses for inlet velocity measurements. Scheimpflug mounts were used to align the focal plane of camera and lens such that the entire field of view is in focus. These add approximately 12 mm to the lens focal length. The configuration used for acquisition of the inlet velocity profile is shown in Fig. 8. The lenses were changed to Tamron Telephoto AF 180 mm f/3.5 Di Macro fixed focal length lenses, while still using the Scheimpflug mounts, for all SRQ data locations. The cameras were placed at approximately 15 deg half angle vertically at each SRQ location as shown in Fig. 9.

Optical rails were used for camera adjustment in the flow-wise direction. Three linear traverses were used for positioning the cameras in the xy-plane and for adjusting the laser position in the x-direction. The repeatability of the traverses was 0.005 mm. After positioning the cameras and laser, davis version 8.3.0 software [24] by LaVision controlled image acquisition, processing, and vector calculation including uncertainty. All PIV measurements made for this study are stereoscopic resulting in three-components of velocity in the laser plane with uncertainty.

PIV calibration for inflow measurements was performed by mounting a two-plane calibration plate inside the wind tunnel at the inlet plane (z = 0.067 m). The 11.875 in², two-plane calibration plate consisting of holes of 0.125 in diameter is evenly spaced every 0.625 in, with a distance between planes of 0.125 in. The plate was mounted directly to the test section, and the laser was then aligned with the plate. The calibration images were dewarped using a two-dimensional third-order polynomial fit. The laser sheet for the inflow measurements was approximately 4 mm thick.

A much smaller, single plane calibration plate was used for the SRQ measurement calibration. A two-dimensional dot array Max Levy Autograph calibration plate was traversed through the laser sheet thickness for multiple plane calibration. The calibration plate consisted of $∅0.25 mm$ dots spaced 0.5 mm apart and was 50 mm × 50 mm in size with a replication error of less than 0.001 mm. The calibration plate was mounted on a Velmex traverse, and images of the target were acquired at −1 mm, 0 mm, and 1 mm within the laser sheet, where 0 mm corresponds to the center of the sheet, and the sheet was 2 mm thick. Every other dot on the calibration target was used for calibration to limit processing time needed for the dot locating algorithm. The same third-order polynomial model was used for the SRQ calibration as the inflow calibration. Due to the geometric constraints, it was not possible to calibrate inside the test section. Therefore, a sheet of clear polycarbonate was mounted to the exterior of the test section inline with the test section wall during calibration to account for refraction issues caused by wide-angle placement of the cameras.

Prior to acquisition, the lens aperture, laser intensity, and seeding density were adjusted to produce optimal particle image quality and particle density. The time between images, dt, was then adjusted such that the displacement was between 4 and 16 pixels to achieve optimal displacement and limit the number of particles leaving the laser plane between images. Particle image diameter and density were determined based on the method presented in Ref. [25].

Images were processed using multipass, window-deformation in DaVis. A geometric mask was applied to the images to eliminate the influence of walls/rods on the correlation after which the average of the dataset was subtracted from the images to remove background information. Two passes of the cross-correlation were completed with 64 × 64 square interrogations windows followed by four passes of 32 × 32 round interrogations windows with flow-based window deformation and symmetric displacement. Each pass was completed with 75% window overlap. Vector postprocessing was then performed to remove vectors with a peak ratio less than 2. Vectors were “strongly removed and iteratively replaced” by a two-pass median filter to correct spurious vectors if the value was larger than one standard deviation of its neighbors. Vector field statistics were then calculated, with uncertainty, for the dataset including mean and standard deviation vector fields, Reynolds stresses, turbulence kinetic energy, and turbulent shear stress.

where B_r is the root-sum-square of all bias sources for the instrument. Uncertainties of the data provided are reported at 95% confidence.

where s_x is the standard deviation of the velocity vectors at a single location for the dataset, and N is the number of vectors at each location. This method assumes that the samples are statistically independent and follow a normal distribution.

All fluid properties used in this study were calculated using the fluid temperature, atmospheric pressure, and relative humidity measurements described in Sec. 4. The fluid specific heat was estimated using the third-order polynomial fit presented in Vol. 6 of Ref. [29]. The fluid thermal conductivity was estimated using a linear interpolation of tabulated data presented in Vol. 3 of Ref. [29]. The remaining properties were calculated using standard psychrometric relationships and are presented in Appendix C of Ref. [1].

The BCs and SRQs presented for these benchmark experiments are presented in Table 1. All of the SRQs reported in this study are directly measured quantities, with exception of the mass flow rate, which is a double integral on the difficulty spectrum in Fig. 2.

Data were acquired at steady-state for two values of heat flux; 400 W/m² and 700 W/m². Each SRQ velocity dataset consisted of 1000 images acquired at 0.1 Hz in order to have statistically independent samples. The inflow velocity fields consisted of 1000 images acquired at 2 Hz. The laser sheet was 2 mm thick for the SRQ velocity measurements and 4 mm thick for the inflow measurements. The laser sheet for inflow measurement was aligned with the streamwise velocity component traveling through the thickness of the laser sheet. This required a larger thickness than usual in order to have sufficient particle displacement between camera frames for low measurement noise. The 2-mm-thickness of the laser sheet for SRQ measurements was required for sufficient particle displacement in the x-direction for the PIV system, while the width of the laser sheet was aligned with the streamwise component as shown in Figs. 3 and 9.

Prior to acquiring the final sets of data at each location, preliminary data were acquired at 1 Hz and processed to check for sample independence. An autocorrelation was applied to the fluctuating streamwise velocity component $w′$ and integrated from zero lag to the time of the zero-crossing after [28]. For 1000 images, the effective number of samples was 100, indicating the need to acquire images 10 times slower. Sample data were then acquired at 0.1 Hz, and an autocorrelation was applied in the same manner. The resulting autocorrelation indicated statistical independence.

The test section inlet velocity data were acquired at the z = 67 mm plane. The laser sheet was oriented parallel to the xy-plane such that the streamwise velocity component was aligned through the laser thickness. Inlet velocity fields and centerline profiles are shown in Figs. 10–13. The inflow velocity field was integrated to determine volume flow rate through the wind tunnel. Inflow mean velocity field measurements are found in InletVelocityField.csv, which is available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection, and inflow Reynolds stresses and turbulence kinetic energy fields are found in InletReynoldsStressField.csv, which is available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection, with turbulence kinetic energy indicated as “TKE” and Reynolds stresses as “RS” followed by the two directions of the stress, i.e., “RSXY” is the Reynolds shear stress in the xy-direction or $u′v′¯$⁠.

Using the volume flow rate and the measured pressure drop across the inlet flow conditioning, the loss coefficient across the flow conditioning system is determined and shown in Table 2. The loss coefficient was calculated by

where the fluid density ρ and the mass flow rate $m˙=ρQ$ were calculated using standard psychrometric relationships for humid air. The density was determined using the mean inflow temperature of the 5 × 3 array of inlet TCs and the ambient pressure and relative humidity measured in the facility. The volume flow rate Q through the inlet was calculated by integrating the w-velocity component of the test section inflow measurement. The calculated mass flow rate through the test section is found in MassFlowRate.csv, which is available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection. Inflow temperature measurements of the inlet TC array are found in InflowTemperatures.csv, which is available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection.

It should be noted that the inflow velocity profile is not a BC but an SRQ for natural convection. The flow is driven by the effects of buoyancy due to the temperature gradient generated by the heated rods and not by an imposed velocity at the inlet as it would for forced convection. The integral of the inlet velocity field, however, was used to determine the mass flow rate through the wind tunnel in order to determine the loss coefficient across the flow conditioner. A velocity profile is a more rigorous validation metric than the double integral of velocity (see Fig. 2) used to determine the loss coefficient. The user of these data should be aware that using the pressure drop is a more accurate representation of the flow conditioner, and the loss coefficient should be used with caution.

Velocity measurements were made on the x = 0 m and x = −0.06 m planes at z = 0.476 m, 0.825 m, 1.174 m, and 1.524 m. The locations of the velocity profiles presented for validation in this study are shown in Figs. 7 and 14.

Each fuel rod contains four cartridge heaters. The spacing between these heaters needed for routing of the wiring outside the test section resulted in a piecewise heating boundary condition on the inner surface of the fuel rods. Table 3 shows the heating locations for each rod section in the z-direction and the power applied to all four rod sections at the respective z-position for both heating conditions. The power listed is the power supplied to all four cartridge heaters at the given z-position and is assumed to be divided between them evenly.

The fluid temperature entering the wind tunnel was measured using the inlet TC array described in Sec. 4 which were positioned on the z = −0.941 m plane. The pressure drop present upstream of the contraction inlet due to the presence of the flow conditioner and seeding array is presented in Table 2.

The heat lost through the polycarbonate walls for both cases was found to be negligible due to its low thermal conductivity. The loss through the aluminum wall, however, was estimated to be 0.72% of the total power input to the system for the 400 W/m² case and 1.22% of the total power input for the 700 W/m² case.

The test section walls are designated as follows and are shown in Fig. 6: The opaque, aluminum wall of the test section (on the y = 0 mm plane) is referred to as the “plate,” the wall opposite the plate (on the y = 304.8 mm plane) is the “top wall,” the wall that lies on the x = +152.4 mm plane is referred to as the “left wall,” and the wall at the x = −152.4 mm plane is referred to as the “right wall.”

The rods are designated by numbers 1–4 starting at the top right corner going counter-clockwise when looking from the test section inlet in the positive z-direction as shown in Fig. 6. The Rayleigh numbers for the 400 W/m² and 700 W/m² cases were 4.5 × 10⁹ and 5.5 × 10⁹, respectively, and the Reynolds numbers were 3450 and 4600, respectively. The Reynolds number was calculated using the test section inlet hydraulic diameter, while the Rayleigh number was based on the rod length.

Time-averaged inflow data for both cases are shown in Figs. 10–13. The velocity contour plots (in Figs. 10 and 12) show the $w¯$ velocity component across the inlet of the test section. The arrow indicates the direction of increasing contour levels, and the dashed line (at x = 0 m) represents the location of the velocity profiles in Figs. 11 and 13 which show all three components of velocity. Turbulence intensity is also plotted with the velocity components along the dashed line and was calculated by

where the turbulence kinetic energy $k=(1/2)(u′u′¯+v′v′¯+w′w′¯)$⁠, and the bulk velocity $w¯bulk$ is the mean inlet velocity, and its value for each case is listed in the figure caption. Turbulence intensity peaks near the walls as the relative magnitude of the streamwise velocity component decrease, causing a greater impact of the $u′u′¯$ and $v′v′¯$ components on turbulence intensity. Note that I in the center of the tunnel for both cases is near zero, and the value reported is likely representative of the dynamic range of the PIV system. Since I is inversely proportional to $w¯bulk$⁠, the I value of the faster moving case appears smaller.

Average turbulence quantities at the test section inlet are presented in Table 4. Here, average turbulent shear stress $τ¯$ was calculated in the davis software [24] by

While some gradients were calculated from the PIV data, through-plane gradients were not. Therefore, the turbulence dissipation rate was calculated using the standard approximation $ε≅Cμ3/4k3/2/ℓ$⁠, where C_μ = 0.09 and the mixing length ℓ = 0.07D_h. Likewise, specific turbulence dissipation rate was determined by ω = ϵ/C_μ/k.

The nondimensionalized rod temperature in the streamwise direction are shown in Figs. 15 and 16 for the Natural400 and Natural700 cases. High thermal conductivity of the aluminum rods resulted in negligible circumferential variation in rod temperature. The rod and wall temperature measurements are found in the RodTemperatures.csv and WallTemperatures.csv files, respectively, which are available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection. As expected, the temperature increases in the streamwise direction as the mixed mean fluid temperature increases. A sharp increase in the temperature between points 3 and 4 (at ∼0.6 m) occurs due to the presence of the grid spacer.

The centerline temperatures for each wall are shown in Figs. 17 and 18. The heat loss from the test section walls due to slightly elevated temperatures was discussed in Sec. 6.

Fluid heating is concentrated in the center of the test section directly down-stream of the rod bundle. The 4 × 4 array of TCs measured the outlet temperature profile, and contour plots are shown in Figs. 19 and 20. Outlet temperatures are found in the OutletTemperatures.csv file, which is available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection.

Twelve velocity fields were measured for each heating case at the SRQ locations shown in Figs. 7 and 14. The velocity profiles at each of these locations were extracted from the fields and are plotted in Figs. 21 and 22 along with the six Reynolds stress components. The mean velocity from each profile w_bulk has been subtracted from the streamwise velocity component for comparison with the other components, and its value is shown in each y-axis label. Each figure contains four plots aligned vertically for comparison with the flow-wise direction, and the z-position is indicated on each plot area. For the sake of brevity, only the Natural700 case will be presented here. The Natural400 case presented similar SRQ trends as the Natural700 case, differing only by magnitude, and may be found in Ref. [1] or in the attached data files.

For the x = 0 m plane (lying between the fuel rods) in Fig. 21, the effect of the swirl elements is strongly present in the $u¯$ velocity component, where y = 0.154 m corresponds to the centerline of the fuel rod bundle. The increasing magnitude of w_bulk in the flow-wise direction is present due to the acceleration of the flow due to buoyancy. The Reynolds stresses in the z = 0.476 m and z = 0.825 m positions are negligible, however, after entering the transition flow regime, the Reynolds stresses increase rapidly as shown in the z = 1.174 m and z = 1.524 m positions of Fig. 21. This phenomenon is in a good agreement with the local Rayleigh number, which indicates that transition to turbulence occurs at approximately half the length of the rod bundle (between the z = 0.825 m and z = 1.174 m SRQ locations).

The strong through-plane swirl effect of the mixing vanes is visible at all z-locations of the x = 0 m plane (Fig. 21), while the in-plane horizontal velocity component, $v¯$⁠, is nearly zero at all z-positions. The normal Reynolds stresses at this x-position are also the only significant contributors to the Reynolds stress tensor with the Reynolds shear stresses being nearly zero at all z-positions.

The fluid flow on the x = −0.060 m plane (just outside the rod bundle) differs greatly from that of the x = 0 m plane and is shown in Fig. 22. The fluid acceleration due to the growth of the boundary layer increases less quickly due to the much larger mass of bulk fluid accelerating outside the fuel rod bundle. This results in lower velocity magnitudes overall as evidenced by the magnitude of w_bulk when compared to the x = 0 m plane. The Reynolds stresses outside the rod bundle also have significantly decreased magnitudes. At the lowest z-location, a slight decrease in the streamwise velocity immediately downstream of the swirl elements is present. The Reynolds stresses in the z = 0.476 m and z = 0.825 m are nearly zero lying in the laminar flow regime. At z = 1.174 m, the Reynolds stresses become nonzero near the center of the rod bundle (the left side of the plot) as the bulk fluid velocity also increases most dramatically in this region. Upon reaching the z = 1.524 m, turbulence diffuses outward into the flow (toward the right of the x-axis in the figure) causing nonzero Reynolds stresses and increased bulk velocity.

The velocity and Reynolds stress SRQ measurements are presented in the attached files as follows. The location of each profile is called out in the filename beginning with “Velocity_” or “ReynoldsStress_” followed by the x- and z-coordinates of the profile in millimeters. Thus, the file “ReynoldsStress_x-60,z1174.csv,” which is available under the “Supplemental Materials” tab for this paper on the ASME Digital Collection, contains the profile of Reynolds stresses along the y-direction for all y at x = −60 mm and z = 1174 mm.

This study has presented experiments that were performed for steady-state natural convection in nuclear fuel rod bundles under two constant heat flux surface conditions (400 W/m² and 700 W/m²). The geometry, BCs, and SRQs are available for download and use for CFD model validation. Stereoscopic PIV was used for unintrusive velocity measurements. The BCs and SRQs considered are listed in Table 1. Uncertainties for all measured and derived quantities were also calculated and are included in each of the attached files. Tight coupling of fluid properties and flow characteristics, along with a flow driven by thermal BCs rather than inflow, make natural convection difficult to simulate. These data provide a benchmark case that may be used to validate CFD models.

The U.S. Department of Energy Nuclear Engineering University Program provided the funding for these experiments under Grant No. 00128493.

A =

cross-sectional area

A_s =

surface area

B_r =

bias error of a measured variable r

C_μ =

turbulence model constant

D =

experimental data

D_h =

hydraulic diameter

E =

validation comparison error

$Ẽ$ =

estimated model error

g =

acceleration due to gravity

h =

heat transfer coefficient

I =

turbulence intensity

k =

turbulence kinetic energy

K =

pressure loss coefficient

k_f =

fluid thermal conductivity

ℓ =

turbulence model mixing length

L =

characteristic length

$m˙$ =

mass flow rate

n =

number of experiments performed

N =

number of vectors

q =

heat transfer

Q =

volume flow rate

$q″$ =

heat flux

s =

standard deviation

s_x =

standard deviation

S =

simulation solution

S_r =

random error on a measured variable r

t_α∕2,v =

quantile for t distribution with v degrees-of-freedom

T_s =

surface temperature

T_∞ =

free-stream temperature

$u¯$ =

time-averaged velocity in the x-direction

U_D =

experimental data uncertainty

U_input =

simulation uncertainty due to input uncertainty

U_num =

numerical model uncertainty

U_r =

uncertainty of a measured variable r

U_val =

validation uncertainty

$Ux¯$ =

uncertainty of single mean velocity vector

$u′u′¯$ =

Reynolds normal stress in the x-direction

$u′v′¯$ =

Reynolds shear stress in the xy-direction

$v¯$ =

time-averaged velocity in the y-direction

$v′v′¯$ =

Reynolds normal stress in the y-direction

$w¯$ =

time-averaged velocity in the z-direction

$w¯bulk$ =

average inlet velocity

$w′w′¯$ =

Reynolds normal stress in the z-direction

x =

transverse direction

y =

aluminum wall-normal direction

$y¯e$ =

estimated mean based on n experiments

y_m =

model system response quantity

z =

streamwise direction

α =

thermal diffusivity

β =

coefficient of thermal expansion

ΔP =

measured pressure difference

ε =

turbulence dissipation rate

ε_s =

surface emissivity

θ =

dimensionless Temperature

ν =

kinematic viscosity

ρ =

fluid density

$τ¯$ =

turbulent shear stress

ω =

specific turbulence dissipation rate