Received 23 November 2015; accepted 4 January 2016; published 7 January 2016
The coolant pressure loss across the rod bundle, as well as across the reactor coolant circuits, is an important factor both for design and safety. The friction pressure losses are often estimated using the equivalent hydraulic diameter (Dh = 4A/P, where A is the cross sectional area and P the wetted perimeter) and a friction factor (λ) valid for circular pipes, even for more complex geometries. It is well known that this way is approximated, i.e. in the rod bundle the flow passage has sharp corners and the flow in the corners is slower than in the core of the channel, so the friction calculated on the mean velocity and the hydraulic diameter is not a good representation of the true flow. Since the seventies several researches have been performed on the pressure drops in different rod bundle shapes. Rehme  developed a method to predict the turbulent friction factor in non-circular channels if the geometry factor of laminar flow is known. The geometry factor (K) is defined as K=Reλ, where Re is the Reynolds number (based on the hydraulic diameter and cross sectional mean velocity) and λ is the laminar friction factor. Cheng and Todreas  developed a friction factor formula for laminar and turbulent flow in different array shapes by fitting the results of Rehme. Lee  developed an analytical method to predict the friction factor for infinite bare rod square and triangular arrays of low pitch to diameter ratio in turbulent flow from the normalized wall shear stress profile. Another approach, proposed for a triangular pipe, was formulated by Ahmed  : it consists in using the friction correlations developed for the circular tube in triangle tube flow calculations replacing the hydraulic diameter with a new characteristic length scale. With this approach the characteristic length has not a generic formulation and has to be developed for each new shape. More recently Gotts  developed a method to calculate the hydraulic resistance multiplier, defined as the ratio of the friction coefficient of a noncircular passage over the friction coefficient of a circular tube of the same hydraulic diameter, to use in the pressure drop calculation. This method in principle is applicable to an arbitrary shape.
In the thermo-hydraulic system codes (e.g. RELAP5), the equivalent diameter is used to model the rod bundle and a dedicated function for friction factor can be activated, but when the analysis is restricted to a bundle sub-channel no friction factor corrections are available. Moreover, friction factor corrections may not be available for new geometries, nor can easily be estimated analytically. In these cases, the equivalent diameter combined with friction factors valid for circular pipes is largely used as first tentative. The objective of this work is to assess the inaccuracy of this approach (when applied, for instance, to a square rod bundle) by comparing the analytical and semi-empirical predictions with results from two different CFD codes: CFX 14 and NEPTUNE_ CFD 1.0.8. NEPTUNE_CFD is a thermal-hydraulic two-phase CFD tool designed to two-phase flow configurations encountered in nuclear reactor power plants. It is developed within the framework of the NEPTUNE project, financially supported by Commissariat à Énergie Atomique (CEA), Électricité de France (EDF), Institut de Radioprotection et de Sûreté Nucléaire (IRSN), and AREVA-NP   . Ansys CFX is a CFD multipurpose commercial code  used to design in a diverse range of industries.
The work involves 2 steps: the CFD code assessment step and the error estimation step. In the first step, the accuracy of the two CFD codes, both in laminar and turbulent regime, is evaluated comparing the friction pressure drop along a smooth circular tube obtained by CFD simulations with the analytically predicted values. In the second step, the errors done in predicting the pressure drop across a square bundle with two different analytical approaches that involve the hydraulic diameter are estimated comparing the analytically predicted values with the CFD simulated ones taken as reference.
The experience gained through the years at the university of Pisa, the current Research group in San Piero, see  - , such as the compliance with OECD guidelines  , have formed the necessary basis of the work.
2. CFD Code Assessment
The CFD code assessment step consists in evaluating the accuracy of two CFD codes: Commercial code CFX14 and research code NEPTUNE_CFD 1.0.8. The results obtained with the CFD simulations of a smooth circular tube are compared with the results obtained with analytical calculations: in the laminar case the exact analytical solution is available (Poiseuille’s flow), in the turbulent case we have only an approximate analytical solution: the friction factor is expressed by semi-empirical correlations.
2.1. Laminar flow
For the laminar flow in a smooth circular tube, in the fully developed region an exact expression for the friction factor defined as per Darcy-Weisbach Equation  exists, it can be easily obtained from momentum balance considerations:
In Equation (1): Vm is the mean velocity (m/s), ρ is the water density (Kg/m3), D is the pipe diameter (m), L is the pipe length (m),
is the Darcy-Weisbach friction factor. In Equation (2) Re is the Reynolds Number of the flow.
This study has considered a smooth circular tube with an internal diameter of 0.006 m, length 2 m, filled with water at 293 K and 105 Pa and an inlet uniform velocity of 0.1 m/s. The Reynolds Number results to be 629 and the pressure drop per length unit 80.62 Pa/m. In Table 1 the results obtained with the two CFD codes are compared with the analytically calculated value.
In the laminar case the CFD codes error with respect to the analytical solution is far below 1% for both codes.
A sensitivity analysis on grid refinement has been performed for both codes. The fully developed zone is reached at z ≈ 0.2 m as expected from empirical correlation Equation (3)  .
The pressure drop per length unit has been calculated in the fully developed region. The velocity profile in a section where the flow in the fully developed has been compared with the analytically calculated profile, see Equation (4) and Equation (5)  , founding a good agreement.
2.2. Turbulent flow
In the turbulent case the pressure drop in the fully developed region can be analytically estimated with Equation (1) combined with an approximate formula for the friction factor. The Colebrook’s empirical equation  is the most common equation used for the friction factor in the turbulent case. Moody  states that the accuracy of this equation for Re > 2000 is well within the experimental error (about ±5% for smooth pipes and ±10% for rough pipes). This relatively large error in the practice is accepted because it is covered by uncertainties in tube size tolerances, form losses, undeveloped flow, etc. The Colebrook’s equation is transcendental and can be solved only iteratively. For smooth tube there are several expressions that approximate the Colebrook expressions, for the analytical calculations two different equations are used: Equation (6), MacAdams  , and Equation (7), Fang  .
The McAdams relation gives an f that is indistinguishable from Colebrook in its range of applicability and gives significant errors when used outside its range   . The discrepancies between Fang and Colebrook formulas is <=0.05%   . Applying to same geometry of the laminar case an inlet uniform velocity of 10 m/s, the Reynolds Number results to be 63 × 103 and the pressure drop per length unit 1.6 × 105 Pa/m using the MacAdams relation, 1.58 × 105 Pa/m using the Fang’s relation. In Table 2 the CFD results obtained with CFX 14 are compared with analytically calculated value.
Table 1. CFD code accuracy for laminar flow.
aThe above results are obtained with the same grid resolution for both codes.
Table 2. CFD code accuracy for turbulent flow.
In the turbulent case the analytical solution is approximate because the friction factor is estimated using an empiric formula. The simulated value is assumed here as the reference for the error calculation because the CFD simulation relies on a RANS 2-equation turbulence modelling following a Low Reynolds approach (i.e. by resolving the turbulent boundary layer instead of applying a law-of-the-wall). Relatively accurate results are thus expected. The Mac Adams relation shows the best agreement with the Low Reynolds simulation.
A sensitivity analysis on grid refinement has been performed. The fully developed zone is reached at z ≈ 0.25 m as expected from empirical correlation Equation (8)  .
The pressure drop per length unit has been calculated in the fully developed region. The velocity profile in a section where the flow in the fully developed has been compared with the analytically calculated profile, see Equation (9) and Equation (10)  , showing a good agreement in the range of validity of the formula.
3. Error Estimation Step
The error estimation step aims at estimating the error done calculating the pressure drop due to friction (in the fully developed flow region) across a square bundle using the Darcy formula combined with an equivalent hydraulic diameter and a friction factor valid for circular equivalent pipes. A fictitious 2 × 2 square rod bundle was considered for the purpose, the geometrical characteristics of which are reported in Figure 1 and Table 3.
The estimation is done by comparing the analytical results with CFD simulation results for the square rod bundle that are taken as reference. The analytical results are obtained with two different approaches: the “Whole bundle approach” and the “Sub-channel approach”.
In the “Whole bundle approach” the pressure drop in the fully developed region is analytically calculated with Equation (1) substituting the circular tube diameter with the equivalent hydraulic diameter of the whole bundle Dh, and using the friction factor of a tube with a diameter equal to the equivalent hydraulic diameter (fDh).
Figure 1. Rod bundle.
Table 3. Rod bundle geometrical characteristics.
The “Sub-channel approach” case consists of:
・ dividing the bundle into sub-channels; taking the symmetry of the geometry into account, three types of sub- channel have been identified, as illustrated in Figure 2: 1-corner, 2-side, 3-central.
・ calculating the hydraulic diameter for each sub-channel (Dhi);
・ solving a system of equations based on the following hypotheses:
o Hp.1: the pressure is uniform over the bundle cross section, then the pressure drop across the three sub- channels is the same;
o Hp.2: mass conservation;
o Hp.3: fully developed flow;
o Hp.4: the shear stresses at the wet interfaces among pairs of sub-channels are ignored.
The pressure drop per unit length is then calculated with Equation (1) by substituting the circular tube diameter with the equivalent hydraulic diameter of the sub-channel and using the friction factor of the equivalent hydraulic tube of the sub-channel. Different average velocities are assumed for each sub-channel. The “Sub- channel approach” wants to improve over “Whole bundle approach” by representing in more detail the internal structure. However, since shear stresses at the interfaces are neglected (Hp4) this approximation is expected to give good results when fluid velocity at the wet interfaces is similar among different sub-channels. In the turbulent case when the core velocity is more uniform, the “Sub-channel approach” is expected to perform better than in the laminar case when it seems more appropriate to consider the channel as a whole.
3.1. Laminar flow
For the laminar flow, in the fully developed region from the analytical calculations we obtain:
・ Whole bundle approach: (ΔP/L)th =80.62 Pa/m.
・ Sub-channel approach: (ΔP/L)th = 55.4 Pa/m.
The two analytical approaches (which are both approximated, as they rely on a friction factor that is valid for circular pipes only) provide noticeably different results. In Table 4 the analytically calculated values are compared with the CFD result.
The CFD result is in between the two analytical estimates. It can reasonably be assumed as more accurate than the two estimates. The best agreement is obtained by the whole bundle approach (11% relative error). Table 5 reports the differences between the mean velocities calculated with the analytical “sub-channel approach” and the simulated value. Figure 3 shows the velocity distribution in the bundle cross section.
In each sub-channel there is a different mean velocity: the highest velocity is reached in the interior sub channel, the lowest in the corner one. The sub-channel approach overestimate the velocity in the central sub- channel while underestimates the velocity in the side and corner sub-channel. The error done by the analytical calculation is higher in the side sub-channel type.
Figure 2. Sub-channels.
Table 4. Analytical calculations accuracy respect to CFD code result.
Figure 3. Laminar case: velocity map in the bundle section.
Table 5. Analytical calculations accuracy respect to CFD code result.
3.2. Turbulent flow
For the turbulent flow, in the fully developed region the analytical estimation using the MacAdams for the friction factor gives:
・ Whole bundle approach: (ΔP/L)th =1.6 × 105 Pa/m.
・ Sub-channel approach: (ΔP/L)th = 1.42 × 105 Pa/m.
The two analytical approaches provide similar results. In Table 6 the analytically calculated values are compared with the CFD result.
Similarly to the laminar case, the two approximated approaches provide noticeably different results. The pressure drop obtained with CFD code, 1.423 × 105 Pa/m, is between the two analytical estimated values, which are 1.6 × 105 Pa/m for the whole bundle approach and 1.42 × 105 Pa/mfor the sub-channel approach, but it is far closer to the sub-channel estimated than to the whole-bundle estimate. There are not enough bases to indicate whether the analytical (sub-channel) estimate or the CFD result is to be expected as the more accurate. In the turbulent flow regime the best agreement is obtained with sub-channel approach and the low Reynolds simulation, the difference with respect to CFD is −0.14%.
Table 7 reports the differences between the mean velocities calculated with the analytical “sub-channel approach” and the simulated value. Figure 4 shows the velocity distribution in the bundle cross section.
As for the laminar case in each sub-channel there is a different mean velocity: the highest velocity is reached in the interior sub channel, the lowest in the corner one. The sub-channel approach overestimates the velocity in the central and corner sub-channels while underestimates the velocity in the side sub-channel, but the absolute value of the error is similar in all the sub-channels, around 3% and smaller than the one found in the laminar case.
Figure 4. Turbulent case: velocity map in the bundle section.
Table 6. Analytical calculations accuracy respect to CFD code result.
Table 7. Analytical calculations accuracy respect to CFD code result.
From the code assessment step, it results that the error done by both CFX and NEPTUNE_CFD codes is:
・ in the laminar case, far lower than 1%, with respect to the exact analytical solution;
・ in the turbulent case, the best estimate available is 2%.
In the turbulent case the analytical solution is approximate because the friction factor is calculated using an empirical relation. The simulated value is assumed here as the reference for the error calculation, The Mac Adams’ relation has the best agreement with the Low Reynolds simulation, the error being lower than −0.6%, while for Fang’s relation the error is −1.9%. The empirical information is affected by a relatively large uncertainty (5% for smooth pipes). It is reasonable to expect that the CFD result to be more accurate than that, but there is no rigorous basis to state that unless quality data is available for validation. For these reasons the best estimate error available is assumed to be the relative difference with Fang.
From the error estimation step, it results that:
・ in the laminar case, the whole bundle approach has the best agreement with the simulated results and the error is 11%. The sub-channel approach gives −23%;
・ in the turbulent case, the best agreement is obtained with the sub-channel approach. The difference with respect to the low Reynolds CFD simulation is −0.14%. The whole bundle approach gives 12%;
Based on the validation step, the best estimate available for the CFD result accuracy is 2% (=relative difference with Fang). Such value is also an estimate of the error made when the turbulent friction pressure losses are estimated by the sub-channel approach.
Todreas  demonstrates that the friction factor for non-circular geometries in laminar flow could not be obtained by transforming circular tube results using the equivalent diameter concept. This result is expected because in the laminar case the molecular shear effects are significant throughout the entire flow cross section so that the governing equations have to be solved for the specific geometry. Todreas also reports the corrected friction factor for different geometries. From Kays  data Todreas reports the product of laminar friction factor and Reynolds number for fully developed flow in an annular channel as a function of the ratio of the two annular radii showing that the value can change by up to 50%. From Sparrow and Loffler  work Todreas reports the data for parallel flow in a triangular array as a function of pitch to diameter ratio of the bundle (p/D) showing high variations of the product of friction factor and Re number with the geometry (from 40 to 320). Rehme  reports a complete set of laminar results for the friction factor as a function of Re number measured by Gunn and Darling  for different rod bundle geometries. Rehme’s article reports data for a square rod bundle with p/D = 1.31 and W/D = 1.16, which is a geometry similar to our bundle configuration Our bundle has a p/D=1,3 and a W/D = 0.65. The Gunn and Darling configuration is the best approximation of our set-up found in literature. From the digitalized figure, fitting the data and the tube approximation, it is possible to estimate an error between experimental data and the tube approximation of around +3.6%. Our error has the same positive sign with a higher value because of the smaller W/D parameter.
In the turbulent case Todreas  estimates that the hydraulic diameter concept can be more accurate in predicting the friction factor because in this case the velocity gradient is principally near the wall and thus the flow channel geometry does not have such an important influence on the friction factor as it has in the laminar case. In the case of rod bundles, taking as reference the experimental data of Deissler and Taylor  for square and triangular rod bundles (p/D =1.12; 1.20; 1.27), Todreas reports that the Mac Adams relation for Re > 104 provides a result that lies within the scatter of the data, but the data show a dependence on p/D that cannot be reproduced by the equivalent diameter concept and the circular tube correlation. Le Tournea  tested rod bundles of square lattice with p/D = 1.12 and 1.20 and of triangular lattice with a p/D = 1.12. These data fall within a band between the smooth tube prediction and a curve of 10% below that for Re in the range 3 × 103 to 3 × 105. Trupp and Azad  tested airflow through triangular array bundles: for Reynolds between 104 and 105 their data at p/D = 1.2 were about 17% higher than the circular tube data. The best approximation to our geometry also in this case is that reported in the Rehme  work. From the digitalized figure, fitting the data and the tube approximation with a power trend-line, it is possible to extrapolate an error between experimental data and the tube approximation of around +3.8%. Our error has the same positive sign with a higher value because of the smaller W/D parameter  .
In the recent years different researchers have used CFD codes to estimate the pressure drop across a rod bundle and the code results have sometimes been validated against experimental data with relatively good results. Smith  developed a methodology for CFD modelling that provide confidence in the numerical results. Capone  performed a CFD analysis of fuel bundle assembly with spacer grids and mixing vane for a 5 × 5 bundle configuration. Then the simulations were compared with the Texas A&M University experimental data with a satisfactory agreement. In and Oh  developed a loss coefficient correlation to estimate the pressure drop across a PWR fuel grid spacer with the mixing vane using a simplified CFD analysis. Yan  performed a benchmark of the CFD modeling of pressure drop for the inlet region of a PWR fuel assembly; the CFD results were compared with the experimental data finding an acceptable agreement.
An accuracy assessment of the friction pressure loss estimation based on Darcy formula combined with an equivalent hydraulic diameter and a friction factor valid for circular pipes has been performed for a square rod bundle in laminar and turbulent flow conditions. Two analytical estimation methods have been considered, the whole bundle and the sub-channel approach. The whole bundle approach gives an error in the order of 11% for both flow regimes. The sub channel approach gives a larger error in the laminar case while in the turbulent case gives the best estimation with an error of few percent. The error estimation has been done comparing the analytical results with the CFD simulations results taken as reference. The CFD results assessment has been performed comparing the CFD results with the analytical estimations in the case of a circular tube and the error is of few percent both in laminar and turbulent case.
The performed assessment covers a limited variety of configurations and thus cannot be generalized (e.g. to different bundle geometries); however, it provides indications that may be helpful when developing thermal- hydraulic analyses by system codes. The CFD tools are confirmed to constitute a useful and reliable support to the thermal-hydraulic analysis of nuclear reactor systems.