Investigating the Effects of Interface Topology on Flow Development in Rod Bundle Supported by Spacer Grid with Split Mixing Vane Using STAR-CCM+ ()
1. Introduction
The use of Computational Fluid Dynamics (CFD) for nuclear applications is promoted internationally by the Organization for Economic Cooperation and Development (OECD) through its Nuclear Energy Agency [1]. The CFD approach for analysis and design is enhanced by the application of best practice guidelines for simulations and which does not invalidate a final but essential step of verifying the CFD simulation results through a few experiments than before. Currently, the application of CFD codes in predicting single-phase flows in rod bundles still faces challenges due to the difficulties in accurately predicting turbulent structures such as secondary flows, vortex shedding and flow pulsations that contribute to inter-sub-channel mixing [1] [2] [3]. Current computational methods can not be fully trusted in complex applications and still require significant development of new models as well as thorough validation of the existing models [4]. Among the contributing factors to these failures is the inability to develop an optimal mesh that appropriately captures the turbulent structures at the relevant points in the flow domain. To obtain any fine details during CFD simulations, the high-density mesh of appreciable quality is required over the extent of the computational domain and especially at the location of interest. This requirement in turn poses high computational resource demand which in many cases is limited.
To overcome this shortcoming, shorter segments of the full domain are modeled with imposed periodic interface topology of fully developed type at inlet and outlet boundaries. The boundaries for a periodic interface are separated in space but can be mapped from one to the other through some constant rotation or translation, and it represents a cyclic repeat of information across the boundaries, so that fluxes that cross one boundary are transformed and applied to the other boundary [5]. The approach effectively approximates a repeating geometry and makes it possible to model an infinitely long pipe by just using a short segment of the pipe [5], and therefore, reduces the computational resource demand for CFD simulations.
In the present study, two interface boundary conditions were imposed to study their effect on the development of flow downstream of a spacer grid with and without mixing vanes. A uniform mass flow inlet condition without the imposition of a cyclic periodic boundary condition at the inlet-outlet boundaries was simulated and the results were compared to data obtained when a Cyclic Periodic interface Boundary Condition (CPBC) or topology was applied at both inlet and outlet boundaries of the domain with a fully developed interface type. Figures 1-3 show a generic design of the split-type spacer grid, the split mixing vane arrangement and the detailed dimension of the mixing vanes respectively that were used to support the rod bundles for the present work.
Figure 3. Detailed dimensions of split mixing vane [6].
Spacer grids with various types of attached vanes to promote turbulence mixing are designed as a support for the fuel rods. When used to provide mechanical support, spacer grids reduce flow-induced vibration [7], and also absorb any vibration impact during the process of loading, shipping and handling of the fuel assemblies. The split-type mixing vanes serve as flow deflectors attached to the spacer grid to improve lateral convection, and hence, enhance mixing between and within sub-channels depending on their inclination in the flow domain. The attached vanes enhance the turbulent processes which lead to suppression of local Critical Heat Flux (CHF) [8] [9]. A commercial nuclear fuel assembly with a mixing vane can enhance CHF by up to 10% [8]. The geometry of the spacer grids also helps maintain gaps between fixed rods which serve the purpose of uniform distribution of hydro-dynamic properties of the flow. Despite the adverse increase of pressure drop in the rod bundle due to constriction of the coolant flow area by spacer grids, form drag and skin friction [10] resulting from the additional wall of the spacer grid, most nuclear fuel developments include the design of spacer grids with or without mixing vanes since their thermal hydraulic impact is too enormous to ignore. There are many efforts being made to develop an optimum spacer grid that can mix the flow effectively in sub-channel geometry [11]. Industrial designs of spacer grids are therefore assessed to satisfy desired optimum performance of sustained turbulence mixing whiles restraining pressure drops as much as possible. To achieve these enhancements, the mixing vane designs are optimized prior to their manufacture and application in the experimental investigation by the use of modern computational fluid dynamic codes. The application of computational fluid dynamic codes for design optimization makes it a less expensive practice than performing a number of experiments until the desired effects of the mixing device on the flow are realized.
The present work aims at investigating the extent of predictability of the mean and turbulent flow properties by the computational fluid dynamic code STAR-CCM+ for the flow of water through a 5 × 5 rod bundle supported by a spacer grid with and without a mixing vane. This was achieved by comparing simulation results with experimental data obtained from the MATIS-H test facility at the Korea Atomic Energy Research Institute (KAERI). The study is relevant in assessing the suitability of the models employed to accurately predict the physical problem under study. Furthermore, the effect of the imposed interface topology on flow redistribution downstream of the spacer grid and hence on computational resource demand was also studied.
2. Methodology
2.1. Experimental
The MATIS-H (Measurement and Analysis of Turbulent Mixing in Sub-channels Horizontal) experiments performed at KAERI (Korean Atomic Energy Research Institute) to provide a set of data on the turbulent mixing due to split-type spacer grid in a rod bundle was employed in the present study. The experimental set up as described by Smith et al. [6] is illustrated in Figure 4.
Experimental data produced downstream of the spacer grid location included axial and lateral velocity components and root mean square velocity fluctuation which were used to measure the extent of predictability of the physical problem by the computational fluid dynamic code STAR-CCM+. Time-averaged data was reported along three line segments, Y1 = 16.56 mm, Y2 = 49.68 mm and Y3 = 82.8 mm as shown in Figure 5 which is modified from the work of Podila et al. [1]. The data was collected at four axial positions of 0.5, 1.0, 4.0 and 10.0 of the hydraulic diameter from the tip of the mixing vane [12] [13].
Figure 4. Schematic diagram of the MATIS-H test facility [6].
Figure 5. Positions for data extraction [12] [13].
2.2. STAR-CCM+ Simulation
Computer Aided Design (CAD) tools embedded in the Computational Fluid Dynamics code STAR-CCM+ were used in the present study to tailor and complete the domain which consisted of a square duct enclosing a 5 × 5 rod bundle arrangement which is supported in turn by a spacer grid with and without split mixing vanes. STAR-CCM+ was also used to generate mesh, set up physics models, solve and conduct post-processing of the results.
2.2.1. Geometry Modeling
The central, wall and corner sub-channels forming the entire domain of interest were modeled individually, duplicated and translated to their respective positions to form the geometry. Spacer grid CAD file made by using CATIA 5 obtained from KAERI was then imported a number of times and the Boolean operation “Subtract parts” applied at each import to subtracts the shape defined by the spacer grid with vanes from the volume of the sub-channel. Upon the Boolean operation, the surfaces of the imported spacer grid were imprinted in the target sub-channel. After identifying the inlet, outlet, and spacer grid boundaries the 3D CAD model was used to create a new geometry part with the named parts forming separate surfaces. The resulting fluid part was then assigned to region with the region mode set to one region per part and the boundary mode to one boundary per part surface. In-place interfaces were automatically formed and correctly defined. The new region has as its boundaries, the part surfaces defining the inlet, outlet, spacer grid and fluid volume. Velocity inlet, pressure outlet and no slip wall boundary types were assigned respectively and a geometry scene created to enable visualization of the final geometry which was used as sub-surface for mesh generation [12] [13].
2.2.2. Mesh Generation
Surface, polyhedral, prism and extruder meshing models were applied. The initial surface was re-tessellated so as to obtain a finer triangulation for generation of a high-quality volume mesh [5]. In the present study, the surface remesher was employed to re-triangulate and optimize the surface quality for volume meshing and also generate the sub-surface for the prism layer mesher. Table 1 shows the surface mesh parameters used in the present study [12] [13].
The polyhedral mesh type was selected considering such factors as, the turn-around time for building the mesh, the desired solution accuracy and convergence rate, the computational memory available, the quality of the starting surface mesh and quality of solution. The prism layer mesh model was used in conjunction with the core polyhedral mesh to generate orthogonal prismatic cells next to wall boundaries. The thickness of each cell layer was calculated using a geometric progression approach based on a constant size ratio. Table 2 shows the specification used for the prism layer in the present study.
Table 1. Surface mesh reference values.
Table 2. Prism layer mesh reference values.
The inlet and outlet boundaries of the completed geometry were extruded using the method of constant rate normal to extend the original domain along the z-axis to the required dimension in Figure 6. After a successful application of the meshing models, the meshed geometry generated is as illustrated in Figure 6. The mesh density and wall y+ are 9.5 million and 11 respectively. Table 3 shows the extrusion mesh reference values [12] [13].
2.2.3. Setting up Physics Models
The physics continuum for single phase liquid (water) used for the present study consisted of the models shown in Table 4.
1) Material Modeling
Single component water which was the material simulated in the continuum was responsible for managing the various thermodynamic and transport properties and physical processes being modeled in the continuum. Table 5 illustrates the properties of water used in the present study [12] [13].
Figure 6. Mesh distribution and extent of quarter domain [12] [13].
Table 3. Extrusion mesh reference values.
2) Turbulence Modeling
The three turbulence models employed for this work are the realizable two-layer k-epsilon model, the non-linear standard k-epsilon model and the SST-k-Omega model [5]. Although it is strictly possible to simulate turbulent flow directly by resolving all the scales of the flow (termed direct numerical simulation), the computer resources required are too large for practical flow simulations [5]. Therefore, a suitable turbulence modeling approach which provides closure of the governing equations in turbulent flows must be employed. Most simulations rely on models that provide closure for the Reynolds Average Navier Stokes (RANS) equations. The challenge in turbulence modeling however is the model of the Reynolds stress tensor in terms of the mean flow quantities. The eddy viscosity model and the Reynolds stress transport models are the two basic approaches employed for modeling the stress tensor which is expressed as [5].
(1)
For the present study, k-epsilon turbulence model which is a two-equation model that solves transport equations of turbulent kinetic energy k and its dissipation rate
was employed. The inherent approach employed to model the stress tensor is therefore the eddy viscosity model which is based on the concept of turbulence viscosity. The most common eddy viscosity model is the Boussinesq approximation which is given by Equation (2) [5].
(2)
where, k is the turbulent kinetic energy and
is the strain tensor computed from the resolved velocity fields
. Although the k-epsilon model also exists in its standard form, the realizable two-layer model was adopted for the present investigation. The realizable two-layer k-epsilon model combines the realizable k-epsilon model with the two-layer model and hence gains added flexibility of an all-y+ wall treatment. An all-y+ wall treatment is a hybrid approach designed to give results similar to the low-y+ wall treatment as
and to the high-y+ wall treatment for
. It will also give reasonable results for intermediate meshes where the cell centroid falls in the buffer layer [5]. The dimensionless wall distance y+ is related to the wall distance y from the first prism layer cell by,
(3)
where,
is wall shear stress,
is density and v is kinematic viscosity.
The realizable k-epsilon model developed by Shih et al. [14] contains a new transport equation for the turbulent dissipation rate
and also, a critical coefficient of the model
, expressed as a function of mean flow and turbulence properties, rather than assumed to be constant as in the standard k-epsilon model. This allows the model to satisfy certain mathematical constraints on the normal stresses consistent with the physics of turbulence. For the realizable k-epsilon model [5], the turbulent viscosity in the transport equation is evaluated by Equation (4).
(4)
where, the coefficient
is no longer constant as with the standard k-epsilon model, but is instead given by Equation (5) [5],
(5)
The coefficients defining the model have the following values [5],
,
,
(6)
(7)
where, the turbulent kinetic energy (k) is derived from a specified turbulent intensity (I) and local velocity (v) by the relation,
(8)
and the modulus of the mean strain rate tensor S and the strain tensor
are expressed as,
(9)
(10)
The two-layer approach, first suggested by Rodi [15], is an alternative to the low-Reynolds number approach that allows the k-epsilon model to be applied in the viscous sub-layer. In this approach, the computation is divided into two layers. In the layer adjacent to the wall, the turbulent dissipation rate
and the turbulent viscosity
are specified as functions of wall distance. The values of
specified in the near-wall layer are blended smoothly with the values computed from solving the transport equation far from the wall. The equation for the turbulent kinetic energy is solved in the entire flow. The realizable models implemented with a two-layer approach, enables resolution of the viscous sublayer when used with fine meshes. The two-layer model is parameterized as a length scale function,
and a turbulent viscosity ratio function,
. The dissipation rate computed from the two-layer model is evaluated as [5],
(11)
The two-layer model blends a one-equation model. According to the Wolfstein one equation model [5],
(12)
where,
and
.
and
The turbulent viscosity is computed from Equation (13) [5],
(13)
where,
and
.
The model coefficients for the cubic constitutive correlation as described in Baglietto et al. [16] were applied for the non-linear standard k-epsilon (NLS-K-E) model.
The under relaxation factors for the turbulence solver that controls the solution of the k-epsilon models (k and
) and the turbulent viscosity solver that controls the update of turbulent viscosity were under-relaxed for the flow solution to be stable and hence converged. The algebraic multigrid parameters were kept at default conditions for the present simulation [12] [13].
3. Results and Discussion
Data was compared over line probe located at the lateral position Y1 and axial position of 1.0 hydraulic diameter from the tip of the mixing vane as shown in Figure 5. Figure 7 shows the location of the line probe for data analysis as presented in this section.
3.1. Validation and Verification of Simulation Results
The extent of predictability of the spacer effect by CFD was determined by comparing trends of the experimental data with the simulation results. Figures 8-11 show the comparison of mean and turbulent flow parameters of the simulation to the experimental for the case of spacer grid with split vane at the spanwise (x-component velocity), transverse (y-component velocity) and axial (z-component velocity) directions respectively. In the legend of Figures 8-11, “SPV” refers to case of spacer grid with split vane, “EXP” refers to trends of experimental data “SIM” refers to trends of simulation data and “RE1” is the
Figure 8. Validation of spanwise velocity data.
Figure 9. Validation of transverse velocity data.
Figure 10. Validation of axial velocity data.
Figure 11. Verification of spanwise velocity data.
Reynolds number at which the simulation was performed and has a magnitude of 3.4 × 104.
In Figure 8, the spanwise velocity distribution from the simulation is compared with that of the experimental data.
As observed in Figure 8, good agreement of the simulation results for the spanwise velocity is obtained. Similar agreement was also realized in the investigation performed by Podila et al. [1], Frank et al. [17], Batta and Class [18] and Lee et al. [19]. Figure 8 showed departure of the simulation results from the experimental especially at the maxima and minima locations. Higher extents of cross-flow are predicted in the gap regions for the simulation results than for the experimental data. The short dash black lines show the gap locations for a pitch of 33.12 mm. Appreciable agreement of the simulation data with the experimental data are made for velocity in the transverse direction at the gap locations as shown in Figure 9 although some level of deviations of the simulation data from the experimental are observed at the maxima and minima.
The axial velocity distribution as depicted in Figure 10 also showed appreciable agreement of simulation data with the experimental observations generally and especially in the gaps. Similar agreements were also reported by Kang et al. [20] and Chang et al. [11].
Results obtained by the application of the realizable k-epsilon model with linear constitutive relation (RTL-K-E) was verified with other turbulence models such as the non-linear standard k-epsilon model (NLS-K-E) as applied in Baglietto et al. [16] with a cubic constitutive relation, and the SST-k-Omega model (SST-K-O). The results are shown in Figures 11-13.
The mesh parameters for the turbulent models in Figures 11-13 are compared in Table 6.
As observed in Figures 11-13, appreciable agreement of the simulation results of the RTL-K-E model with the NLS-K-E model and the SST-K-O models is evident for the mean flow parameters. Holloway et al. [21] in their investigation attributed the success of the SST-K-O model to their near wall treatment
Figure 12. Verification of transverse velocity data.
Figure 13. Verification of axial velocity data.
Table 6. Comparison of mesh parameters.
which relies on damping of the turbulent viscosity in low Reynolds number regions rather than a two-layer wall treatment as implemented in the RTL-K-E model. A better prediction by the SST-K-O model and the non-linear models are also attested to in the works of Holloway et al. [21], Ruifeng et al. [22], Gandhir et al. [23] and Horvath et al. [24]. Consistent with the findings of Yang et al. [25], Zhu et al. [26] established in their analysis appreciable prediction by the standard two-layer k-epsilon model than the realizable two-layer k-epsilon model.
3.2. Effect of Interface Topology on Flow Distribution
3.2.1. Case of Spacer Grid without Mixing Vanes
Figure 14 and Figure 15 show a comparison of curves of axial velocity at a lateral position of Y1 as shown in Figure 5 for the case of flow of water through a square duct that encloses a 5 × 5 rod bundle arrangement and a spacer grid without mixing vanes. For each graph, simulation data were compared for the
Figure 14. Effect of interface topology on axial velocity distribution before spacer grid location [No-Vane].
Figure 15. Effect of interface topology on axial velocity distribution after spacer grid location [No-Vane].
case that involved the imposition of periodic boundary condition (NV-SIM-CPBC) and the other with no imposed periodic boundary condition (NV-SIM) at the inlet-outlet boundaries at axial positions of 0.5 hydraulic diameters Before Spacer Grid (BSG) and 1.0 hydraulic diameter After Spacer Grid (ASG).
Large variations in the trend of axial velocity are noticed between the two interface topology types Before Spacer Grid (BSG) as observed in Figure 14. The deviation has however narrowed after flow through the spacer grid as seen in Figure 15. The upstream effect on flow due to differences in interface topology are however still prominent downstream of the spacer grid without mixing vane, as observed in Figure 15. Thus, much longer segment of the domain is required in order to nullify the entrance effect on the downstream of flow parameters.
The comparison of the interface topology effect on turbulent intensity in the presence and absence of mixing vanes is as shown in Figure 16 and Figure 17. As in Figure 16, large deviations in turbulent intensity profile of the two different topology types is observed at entry of flow to the spacer grid which is depicted by Figure 16. At exit of the spacer grid, although the variation in trend reduced drastically, an appreciable difference is profile especially at the maxima is noticeable in Figure 17. This observation indicates that in the absence of mixing vanes, the effect of the imposed inlet-outlet interface topology type on mean and turbulent flow parameters cannot be neglected immediately downstream of the spacer grid.
3.2.2. Case of Spacer Grid with Split Mixing Vanes
As observed for the case of spacer grid without attached mixing vanes, the inlet-outlet interface topology effect was evident even after flow through the spacer grid. In the presence of mixing vanes however, the deviations observed in the curves upstream were almost completely cancelled downstream of the spacer grid as shown in Figure 18 and Figure 19. Thus shorter segment of the domain than required for the case of flow through a spacer grid without vane is needed
Figure 16. Effect of interface topology on turbulent intensity distribution before spacer grid location [No-Vane].
Figure 17. Effect of interface topology on turbulent intensity distribution after spacer grid location [No-Vane].
Figure 18. Effect of interface topology on axial velocity distribution before spacer grid location [Split-Vane].
Figure 19. Effect of interface topology on axial velocity distribution after spacer grid location [Split-Vane].
to conceal the effect of the inlet-outlet interface topology on flow development downstream of a spacer grid with attached split mixing vanes. The computational resource demand when modeling flow through a spacer grid with split mixing vanes would therefore be appreciably lower than in the absence of mixing vanes.
Comparison of Figure 20 and Figure 21 also show a complete nullification of the interface topology effect on the turbulent parameter in the presence of mixing vanes at exit of the spacer grid.
The independence of the mean and turbulent flow parameter distribution at 1.0 hydraulic diameter downstream of the spacer grid on interface topology for the case of flow of water through a spacer grid with split mixing vanes is attributed to the additional pressure drop at this location due to the mixing vanes. Table 7 shows line averaged pressure drop values across the spacer grid obtained over line probes positioned at spacer inlet and at 1.0 hydraulic diameter from the exit face of the spacer grid. The split vane effects (ratio of pressured drop in the presence of split vane to pressure drop in the absence of split vane) at lateral
Figure 20. Effect of interface topology on turbulent intensity distribution before spacer grid location [Split-Vane].
Figure 21. Effect of interface topology on turbulent intensity distribution after spacer grid location [Split-Vane].
Table 7. Vane effect on pressure drop across spacer grid.
positions Y1, Y2 and Y3 as in Figure 5 for the four Reynolds numbers employed in the present study are also tabulated.
At the lateral position Y1, Table 7 shows an increase in pressure drop due to the mixing vane by 40% at Re1, 44% at Re2, 48% at Re3 and 50% at Re4. These large increases in pressure drops that are also observed at lateral positions Y2 and Y3 led to a reduction in momentum and subsequent deceleration in the flow such that much time was available to attain a more developed redistribution just behind the mixing vanes than for the case of spacer grid without vanes for which no deceleration in the axial flow due to mixing vanes is present.
4. Conclusions
1) In the present study, the suitability of STAR-CCM+ to predict mean and turbulent flow parameters for a new design of spacer grid with split mixing vanes which supported a 5 × 5 rod bundle arrangement and enclosed in a square duct was assessed.
2) Polyhedral meshing model was employed to discretize the flow domain and the segregated flow model which solves the flow equations, one for each component of velocity herein termed as spanwise (x-component), transverse (y-component) and axial (z-component) and one for pressure in an uncoupled manner was used to model the single phase, steady state, incompressible flow in the present study. The segregated fluid isothermal model was applied as the energy model to provide a constant temperature field of 35 degrees for all models that required temperature. The realizable two-layer k-epsilon turbulent model was adopted to model the transported variables of production and dissipation of turbulent kinetic energy.
3) Comparison of simulation results with experimental data showed a good prediction of the mean flow parameters by the models employed in STAR-CCM+ which indicates the suitability of the commercial software for the analysis of the physical problem.
4) Comparison of mean flow parameters before and after the passage of flow through spacer grids with and without mixing vanes showed the effect of the mixing vanes to quickly bring the flow to an equal redistribution downstream of the spacer grid irrespective of the inlet-outlet interface topology. For the case of the spacer grid without mixing vanes, the inlet-outlet interface topology influenced the distribution of the flow downstream of the spacer grid. The observation indicates that a much longer inlet segment of the domain is required to completely nullify the effect of the inlet-outlet interface topology on flow distribution in the absence of mixing vanes than in the presence of mixing vanes which may lead to a relatively higher demand for computational resource. The observation was associated with the additional pressure drop caused by the presence of mixing vanes and the retardation of flow by the mixing vanes resulting from flow interaction of the axial flow by lateral convection.
Acknowledgements
CD-Adapco is acknowledged for providing technical support for this work.
Abbreviations Definition
OECD Organization for Economic Co-operation and Development
NEA Nuclear Energy Agency
MATIS-H Measurement and Analysis of Turbulent Mixing In Sub-Channels-Horizontal
EXP Experimental
Re Reynolds Number
SPV Split Vane
NV No Vane
SST-K-O Shear Stress Transport K-Omega
NLS-K-E Non-Linear Standard K-Epsilon
RTL-K-E Realizable Two-Layer K-Epsilon
BSG Before Spacer Grid
ASG After Spacer Grid
CPBC Cyclic Periodic Boundary Condition
Nomenclature
Model Coefficient
I Identity Matrix
k Turbulent Kinetic Energy
Length Scale
Wall Distance-based Reynolds Number
Modulus of the Mean Strain Rate Tensor
S Strain Tensor
Turbulent Stress Tensor
Lateral Fluctuating Velocity in X-Direction
Lateral Fluctuating Velocity in the Y-Direction
Axial Fluctuating Velocity in the Z-Direction
y Wall Distance
y+ Wall Dimensionless Distance
Greek Letters
Turbulent Dissipation Rate
Turbulent Viscosity
Density
Wall Shear Stress