CFD Turbulence Study of PWR Spacer-Grids in a Rod Bundle

Nuclear fuel bundles include spacers essentially for mechanical stability and to influence 
the flow dynamics and heat transfer phenomena along the fuel rods. This work presents the analysis of the turbulence effects of a split-type and swirl-type spacer-grid geometries on single phase in a PWR (pressurized water reactor) rod bundle. Various computational fluid dynamics (CFD) calculations have been performed and the results validated with the experiments of the OECD/NEA-KAERI rod bundle CFD blind benchmark exercise on turbulent mixing in a rod bundle with spacers at the MATiS-H facility. Simulation of turbulent phenomena downstream of the spacer-grid presents high complexity issues; a wide range of length scales are present in the domain increasing the difficulty of defining in detail the transient nature of turbulent flow with ordinary turbulence models. This paper contains a complete description of the procedure to obtain a validated CFD model for the simulation of the spacer-grids. Calculations were performed with the commercial code ANSYS CFX using large eddy simulation (LES) turbulence model and the CFD modeling procedure validated by comparison with measurements to determine their suitability in the prediction of the turbulence phenomena.


Introduction
A PWR requires a design with enough safety margins regarding the critical heat flux (CHF). The value of the CHF is largely altered by the presence of the spacer-grids geometry in the fuel assemblies and is predicted generally by means of empirical departure from nucleate boiling (DNB) correlations.
In the last years, the influence of flow obstacles on the CHF has been studied. Pioro et al. [1] confirm, with various types of flow obstacles, the strong influence on CHF of the distance from an upstream flow obstacle and how the CHF enhancement decays exponentially with the distance from the flow obstacle. Furthermore, several studies that include experimental research [2][3][4] and computational simulations [5][6][7][8][9] have been conducted in order to study the effect of the spacer-grid on the coolant.
In recent years, the use of CFD numerical tools in nuclear engineering area has grown rapidly and some codes based in those numerical techniques are well-established, and stateof-the-art of them is employed in the nuclear engineering area. CFD tools are widely used to provide supplementary information relevant to safety margin, and especially their use is becoming widespread in modeling situations showing a strong multidimensional and turbulent flow behavior. However, it is not very mature in some specific areas, and hence the possibility of validating the complex physical phenomena happening in nuclear structures as spacer-grids is highly appreciated as these results could be of great value for future studies of spacer-grid including heat transfer from the rods and as a basis of spacer-grid simplifications.
Turbulent flow structures in the subchannels of a rod bundle are largely influenced by the design of the spacer-grids and modifying the heat transfer from heat structure to the coolant in operational transients or accidents scenarios. Local characteristics of the turbulent flow in subchannels of a rod bundle are sensitive to numerical and physical models used in CFD analysis of which the tools should be widely validated using a proper set of experimental data.
The simulation of the spacer-grid requires advanced turbulent model. Since a long time ago, many contributions have been focused on the use of LES technique showing its capability reproducing the important features of wallbounded turbulent flows. Since 1963, the Smagorinsky [10] formulation for resolving the subgrid scales has been used and the near-wall predictions have improved [11][12][13][14][15]. Thence, further progress and the availability of faster computers have permitted the use of LES turbulence models in industrial applications as can be shown in several detailed reviews over this year [16][17][18][19][20].
The second international benchmark exercise on the turbulent mixing in a rod bundle is the OECD/NEA-KAERI rod bundle CFD benchmark exercise based on the MATiS-H (Measurement and Analysis of Turbulent Mixing in Subchannels-Horizontal) experiments [21]. These experiments provide data to reveal fundamental characteristics of the turbulent mixing in subchannels and confirm the use of the CFD codes as a tool for later use in modeling spacer-grids characteristics, like pressure drop, quantifying the CHF margin reliably for normal operation condition, and operational transients, and also allowing eventually the use of CFD codes for predicting the DNB under accidental conditions.
The following sections contain an exhaustive description of how the case is configured to obtain an even balance between maintaining a low computational demands and fulfilling the requirements according to the NEA Best Practice Guidelines (BPG) [22] of the CFD simulation. A description of the experimental facility and the conducted experiments are described in Section 2 Benchmark Test Description; the results obtained in those experiments will be used for the validation of the CFD models. A comprehensive description of the model is described in Section 3 CFD Simulation Setup; it accounts for the whole modeling process used to obtain results to validate with the experiments. These models are simulated using a LES turbulence model within reach of a regular workstation without significant loss of accuracy in the predicted results. Results for spacer-grids with split and swirl mixing vanes were presented in the blind benchmark, and results of velocity, RMS velocity, and vorticity at different streamwise and spanwise locations downstream of the vanes of the spacer-grid are used for the validation of the model in Section 4 Results and Discussion as well as an analysis of the resulting flow characteristics in the rod bundle. Finally, Section 5 Conclusions summarizes the conclusions and remarks of the work presented in this paper.

Benchmark Test Description
The MATiS-H, located at the Korea Atomic Energy Research Institute (KAERI), Daejeon, Korea, is illustrated in Figure 1. The test rig consists of a channel of 170 mm side length and 4.670 mm long and a 5 × 5 rod bundle of 25.4 mm of rod diameter and 3.863 mm long operating at Re ∼ 50000. The hydraulic diameter of the channel, , is 24.27 mm. The spacer-grid with mixing vanes, 2.6 times larger than the size of PWR fuel bundles, is located inside the channel for generating lateral turbulent mixing in subchannels.
The movable flow straightener allows expected identical inlet flow conditions upstream of the spacer.  is set at 100 to have a fully developed flow profile at the inlet to the spacer-grid. A measurement section is fixed at a position 10 mm upstream of the end of the rod bundle. The spacer-grid, of dimension , can be moved axially to increase the downstream distance of the measurement section to determine the distance upstream of the spacergrid ( in Figure 1). Different experiments are carried out varying from 0.5 to 10 to cover the decay of the turbulence produced by the spacer-grid.
The experiments have been performed for two different types of mixing vanes: split-type (Figure 2(b)) and swirltype (Figure 2(a)). In both cases, axial and lateral velocities were measured with Laser Droppler Anemometry (LDA) probes for , , and measurements. Turbulent intensities and vorticity in subchannels were then evaluated from the velocity measurements.

CFD Simulation Setup
In this section, the simulation setup developed for investigating spacer-grid effects in MATiS-H is presented. An extensive description of geometry modeling, meshing, and CFD setup of the physical models will be provided.

Computational Domains.
The main purpose of the simulation is to obtain results downstream of the spacer-grid. The model was classified in different domains for better adapting to the experimental scenario. The proposed decomposition of the test model is outlined in Figure 3. (ii) Domain A is used to compute suitable inlet profiles for Domain B as the known value at the inlet of the test rig is the mass flow rate. It is a bare rod bundle of 90 and starts from the movable flow straightener upstream of the spacer-grid (see Figure 1).
(iii) Domain C was modeled in order to determine the influence of the proximity of the measurement section (10 mm). The simulation results corroborate the low influence of the immediate vicinity of this domain on the measurement section. As a result, Domain C was not considered in the final CFD simulations with spacer-grids. Domain C does not have such influence, and all the measurements to validate the model can be done with only one simulation per spacer-grid of 10 downstream of the spacer-grid.

Boundary Conditions.
Assuming the decomposition mentioned in the previous section, Figure 4 shows a schematic representation of the boundary conditions for the CFD simulations. The inlet boundary condition for Domain B is established with the velocity, turbulence kinetic energy, and turbulence energy dissipation data obtained at the end of Domain A. The uniformity of the flow at the inlet of Domain A must be guaranteed by the movable flow straightener (Figure 1).
The assumed fully developed flow at the inlet of Domain B was tested evaluating the turbulence decay in the streamwise direction in the 90 of the bare rod bundle and it was concluded that the values of kinetic turbulence energy are almost stabilized at the final of this domain.
The outlet boundary condition is always established as a constant pressure for both simulations. It has been tested that the pressure drop on the isothermal system does not influence the water properties, as density and viscosity. The Δ value would be neglected and the pressure value obtained located at the inlet of the test rig ( Figure 1) is used as a constant pressure outlet boundary condition. To obtain the exact pressure value at the outlet, it requires the simulations of straighteners, rod bundle supports, flow breaker, and out end, but it is not the focus of the present paper.
The boundary condition at the wall is established as automatic near-wall treatment, which automatically switches from wall functions to a low-Reynolds near-wall formulation as the mesh is refined [23].

Geometry Modeling and Simplifications.
The experiment proposed in this benchmark has a relative huge test section considering the flow characteristics happening there. Consequently, a modeled simulation of the whole experiment becomes heavier and the efforts to reduce computational requirements are appreciated. Hence, a set of measures and simplifications have been performed to achieve a feasible simulation as described in this section.

Grid Strap Simplification.
The specifications of the spacer-grid geometry model show a gap between walls and the end of chamfer strap. This narrow distance can lead to numerical errors; however, the reduction in the cross section cannot be neglected. For such cases, a geometrical simplification in the grid straps has been performed. Straight straps are considered, and full intersection between rods and external wall channel is accomplished as shown in Figure 5.

Periodicity
Pattern. The channel itself, without spacergrid influence, has both horizontal and vertical symmetry. However, the presence of spacer-grids produces a local variation of the flow depending on the vanes direction. It leads to an asymmetric flow downstream of the vanes that complicates the model simplification with symmetries and requires the meshing of the entire model. The spacer-grids used in the experiments of the benchmark have an asymmetric geometrical configuration of the vanes. Nevertheless, during the preliminary steps of the geometric modeling procedure, a repetitive pattern was detected. The configuration of the vanes produces the opposite effect in one half of the domain than in the other half. In Figure 6, it can be noticed the vanes directions for both types of spacer-grids. Figure 6 shows the bending direction of the vanes for the split-type spacer-grid marked with R (Right), L (Left), U (Up), and D (Down). Looking at the vanes colored in red and applying a rotation of 180 ∘ around the center of the channel, it can be appreciated that each vane has the opposite direction regarding the vane in the same new position (e.g., up-down and left-right). It means that things happening in one half of the domain can be considered without modeling the other part of the geometry.  Science and Technology of Nuclear Installations  A periodic boundary condition has been applied to this model. This boundary condition is used when the physical geometry of interest and the expected pattern of the flow solution have a periodically repeating nature. The interface model defines the way the solver models the flow physics across the interface.
Only half of the domain can be modeled applying the required cut (from the center of the channel) in the geometry and implementing boundary interface using "rotational periodicity" [24]. The two sides of the periodic interface can be mapped by a single rotational transformation about an axis. The axis of the rotational transformation must be defined in the streamwise direction with the center in the middle of the channel. To test the correct performance of that simplification, the whole model without periodicity conditions and the new model with half domain with the rotational periodicity were simulated (Figure 7). The results show the velocity profiles for the case with periodicity boundary conditions matching those obtained with the entire model. Therefore, it will save half of the mesh nodes and much more in terms of computing time as the relation between the number of nodes and the run time does not increase linearly.

Grid Generation.
The mesh was created using the ANSYS ICEM CFD 13.0 software and tested its quality in the aspect ratio, angles, skewness, and other quality criteria. The mesh type of Domain A is based on a block mesh as the absence of geometry irregularities allows higher quality mesh. However, Domain B containing the spacer-grid has irregularities that hinder obtaining a suitable mesh as possible.
Aware of the importance of meshing, one of the major parts of the effort devoted to the project has focused on obtaining a mesh as uniform as possible with an optimum size to cover the variations of fluid characteristics.
As demonstrated in [25] using a grid following the Taylor microscales, the spanwise mesh resolution and the cubicity of the cells are a crucial feature of the grid in channel flows. In that contribution, special attention has been taken in each subchannel having the consistent spanwise and streamwise mesh resolution, a uniform distribution, and a beneficial cell topology. Due to the different flow characteristics, the mesh subdomain upstream of the spacer-grid has an axial size twice of the spacer-grid and upstream regions. The cross section mesh distribution is the same in both cases.
In addition, since the complexity of the geometry, a large number of nodes should be considered. An agreement between the number of nodes (meeting the BPG), which affects the accuracy of the prediction, and the computational resources required has to be found.
Domain B mesh was generated with the following procedure: (1) scaling the geometry to reduce the cell size in the axial direction to obtain an exhaustive control of the cell axial size (the coefficient to scale varies for the upstream and spacer-grid areas); (2) generation of the volume mesh with the octree mesh method; (3) creation of the specified prism elements in successive layers away from the wall. A total of 10 layers with an exponential growth rate between them are applied; due to the narrow distance between some spacer-grid walls, the number of layers must be adapted to avoid bad quality of the elements and collisions between layers; (4) performing a "Tetra to Hexa" conversion with the ICEM tool "12 tetra to 1 hexa" to save cells and provide spatial uniformity; (5) resizing the mesh to its original axial size.
As a result of this procedure, two meshes were obtained with the characteristics shown in Table 1.
The result of this mesh procedure and its influence on the flow will be shown. Figure 8(a) establishes the intersection planes for the mesh analysis, and Figure 8(b) shows the mesh distribution in the cross plane.
Regarding the cross section mesh, an ordinary subchannel has been selected and the zoom detail in the near wall  presented ( Figure 9). The wall detail shows the fine mesh near the wall and the complex transition between the layers and the core region. The height of the first prism layer is defined to reach the target + and a minimum number of prism layers are established. The height ratio (ratio in which the height of the next layer is specified with exponential law) is dynamically adapted in order to assort the geometry irregularities.
The mesh of the two longitudinal planes shown in Figure 10 is as follows: one intersects rods and the other intersects spacer-grid wall.
The mesh distribution near the wall must deal to reach the criteria of NEA BPGs in a good agreement between the + and the Courant number. The range of + for both spacergrids is from 0.06 to 2.98 with average values of 1.13 for splittype and 1.08 for swirl-type. Figure 11 shows the local effects on the Courant number of the flow near the wall. Far away from the wall near 1 values are obtained.
The vane configuration has a direct influence in the meshing procedure, particularly, in the transversal and longitudinal mesh sizes. In the current simulation, the maximum cell size has been considered equivalent for both spacergrid types. It is required to highlight the influence of the discretization in order to know the assumed error in benefit of the computational resources and know the possible weakness points of the used model. In case of the swirl vane, gradual changes of the flow occur, but for split vanes higher gradients are present. Consequently, for the same spatial discretization, the vane area shows more local numerical irregularities ( Figure 12(b)) for the split-type rather than for the swirl-type (Figure 12(a)).
A detailed zoom of the mesh near the vanes for split and swirl spacer is shown in Figure 13.

Turbulence
Modeling. The turbulence model selected for solving Domain A is the RANS shear stress transport (SST) model [23] because of the simplicity of the flow in this area. The same turbulence model has been used for the simulation to initialize values in Domain B for the LES simulation. The Smagorinsky subgrid-scale model [23] was used because of its solid mathematical and physical background instead of the WALE model, although theoretically it is more appropriate for wall bounded flows. To account for the near-wall effect, the turbulent viscosity is damped by means of a length minimum and viscosity damping functions.  scheme (2nd or higher order) has been chosen when using the SST turbulence model. For the LES simulation, we have used the central differencing scheme and second order backward Euler transient scheme. The target maximum number of coefficient loops per each time step has been selected to optimize the computational time. The selected value of 4 maximum number of coefficient loops showed optimal computational time without loss of precision. The convergence criteria for the RMS residual were set to a target value of 1 ⋅ 10 −6 for both turbulence models. For LES turbulence model, as an unsteady simulation, a total time must be provided and it is needed to specify enough time to provide a valid statistic. An optimal time has been determined based on temporal evolution of some monitored variables in a set of points. The statistical begins at a sufficient time that it is not influenced by the initial values. The summary of simulations for SST and LES turbulence models is shown in Table 2.

Results and Discussion
The CFD results presented were run in a 2 × 6-core Intel Xeon E5645 at 2.40 GHz workstation using ANSYS CFX Academic Research, Release 13.0. For Domain B simulation, Table 3 shows a summary of times and related information for each spacer-grid type. A fixed time step was selected to permit a resulting RMS Courant number near 1. A total time large enough to obtain statistics steadiness have been determined. The next sections will show a set of results to study the flow pattern modifications due to its pass through the spacergrid elements. The validation of the CFD model through comparison with the experimental data is performed at strategic locations in the domain considering both mean and profile values (Section 4.1). As additional information, some figures are obtained from the CFD results in order to see the flow behavior produced by the spacer-grids as a

Comparison between CFD Simulation and Experimental
Results: Model Validation. The following pages show the comparison between experimental data and the CFD results. The comparison includes profiles for time-averaged velocity and time-averaged RMS values of the fluctuating components for all three velocity components at two different distances from the wall of wall = 3.71 mm and core = 68.44 mm. In order to observe the decay of the turbulence, three different axial positions 1 , 4 , and 10 from the origin of the system of coordinates as marked in Figure 14 are analyzed. Furthermore, an ordinary subchannel has been chosen in order to analyze the local effects as the circulation and vorticity in the individual subchannel.
Generally, the split spacer-grid produces more turbulence than the swirl-type spacer-grid and its effects go to 10 downstream. This situation differs depending on the distance from the wall being in the center of the channel more influenced by the neighboring vanes.
The comparison of the CFD results with the measured data at the two different distances from the wall for mean and RMS velocity values in the three axial positions for the three components , V, and is shown in Figures 15 and 16 for both swirl-type and split-type spacer-grid spacers. The horizontal axis represents the position from the center of the channel to the wall and is normalized with the rod-to-rod pitch distance ( = 33.12 mm). Mean and RMS velocities in the vertical axis are normalized with the streamwise bulk velocity ( bulk = 1.5 m/s). Sets of charts in Figures 15 and 16    The results obtained in the simulations show a good capability to capture the turbulence phenomena and the procedure of production of turbulence and dissipation. Velocity profiles for all the components as those shown in Figures 15(a) to 15(c) and 16(a) to 16(c) clearly show a good agreement for peaks and valleys, while in Figures 15(g) to 15(i) and 16(g) to 16(i) small differences are presented; nevertheless, the trend remains fitted to the experimental data. It is attributed to a need of a local grid refinement in the mesh for the split spacer-grid where high velocity gradients appear. The assumption of a small error was accepted in benefit of the computational resources required as well as the influence on the time step in order to keep the Courant number under the required value. These mentioned figures also show how the turbulence model used is capturing with the same agreement the velocity profiles near the wall and far away from the wall. The dissipation turbulence downstream of the spacer-grid is noticed from distance 1 to 4 and 10 where the profiles are almost flattened.
RMS fluctuating velocities show general overpredicted values especially near the focus of turbulence (Figures 15(d), 15(j), 16(d), and 16(j)). The prediction improves as the fluctuation decreases downstream of the spacer-grid ( Figures  15(e), 15(f), 15(k), and 15(l)). The turning points determining the relative maximum and minimum values seem to fit the experimental data. As mentioned before, the grid size plays an important role, but there is a big influence of the subgrid-scale model and the wall damping functions requiring further investigation. More experimental results are required to validate the turbulence behavior between the rods in a subchannel.
Values of circulation expressed as ∯ , where is the streamwise component of vorticity, are compared with experimental data for the designated subchannel (see Figure 17). As expected, the subchannel circulation for the split-type spacer-grid is around twice bigger than the swirltype spacer-grid. Simulation results have the same trend in the split-type, but the swirl-type has a discrepancy at 1 .

Study of the Flow Structure Produced by the Spacer-
Grids. Once the model has been validated, some information about the flow behavior produced by the spacer-grids will be studied. Following figures extracted from the CFD results are illustrated with reflection applied in the postprocessing stage for a better visualization. Figure 18 shows the effect of split and swirl spacer-grids on the flow. Streamlines representing the flow leaving the spacer around only one rod are illustrated. It is appreciated as for the split-type Figure 18(b) the flow influences directly the other subchannels, while for swirl-type Figure 18(a) the effect remains in the same subchannel.
In Figure 19, the turbulence dissipation process in three different measuring planes is appreciated. Values of timeaveraged velocity, time-averaged RMS values of the fluctuating velocity in the streamwise component are presented.
In Figure 19, it is highlighted that regardless of the type of spacer-grid, the flow distribution at 10 is quite similar, but the swirl has a more symmetric distribution of velocity and RMS fluctuation values since 4 downstream from the spacer-grid. The split-type shows some undesired minimum     peaks at the corners. For the mean velocities, it can be observed that the split type influences the flow far after the vanes. At 10 , the presence of high gradients at two opposite corners can be noticed, while the swirl maintains almost uniformity from 4 . The RMS are quite stable in both cases and seem more similar. In case of the axial vorticity evolution, there are higher values in the subchannels of the split-type in comparison with the swirl-type.
In order to illustrate the local effect that the spacergrids have in each subchannel, the evolution of the mean streamwise vorticity along the downstream of the spacer is shown in Figure 20 for the designated subchannel (see Figure 14). The turbulence produced by the vanes of the split-type spacer-grid affects fairly the near-wall region, and it could have a positive effect on the local heat extraction from the wall. In case of the swirl-type, the vorticity gradient is more uniform increasing the capability for homogenized temperature distribution in the core of the subchannel.

Conclusions
A CFD model for a rod bundle with a spacer-grid has been performed. The design procedure of modeling, meshing, physical model applied, and considered assumptions to model the related experiment has been described in the previous chapter, and the results and discussion were presented for two types of spacer-grid, split, and swirl mixing vanes.
It has been demonstrated that the described procedure to simulate PWR spacer-grids in rod bundles is valid to obtain realistic fluid mechanic predictions. The LES turbulence model implemented on a general purpose commercial code as the one used in this paper was satisfactory predicting the turbulent results. The simplifications added in the modeling stage reduce the simulation time substantially and allow a LES simulation with the performance of a regular workstation.
The turbulence production and dissipation evolution obtained in the simulation fit well the experimental data.
Velocity and RMS fluctuating evolve with a similar trend as the experimental data for two distinct spacer-grids. The circulation data in the subchannel adequately conforms the evolution along the flow direction verified with the experiments. Therefore, an analysis focusing on the mixing temperature can be done as the turbulence phenomena are being modeled properly. Further work should focus on the influence of the turbulence enhancement due to spacer-grids in the temperature profiles in the subchannels and the resulting CHF. Furthermore, the validated CFD model will be useful as a basis to create correlations to consider the spacer-grid turbulence generation without its detailed model to use in subchannel codes.