Benchmarking the Performance of the ANSYS-FLUENT Standard k- ε Turbulence Model in Fluid Flow and Heat Transfer Predictions for Complex Flows around Circular Pin-Fins Using Various near Wall Functions

This study compares CFD analyses of the fluid flow and heat transfer phenomena in a popular pin-fin geometry of X/D = 2.5, S/D = 2.5 and H/D = 1 for a range of Re = 5,000 to 30,000 to those from experiment to aid in the benchmarking the performance of the CFD code FLUENT. The CFD analyses use three ANSYS-FLUENT (version 13) near wall treatments available within the code: 1) the Standard Wall Function (SWF), 2) the NonEquilibrium Wall Function (NEWF) and 3) the enhanced wall treatment. Experimental data used in this study were obtained from two papers: 1) by Chyu et al. (1998) for heat transfer predictions and another 2) by Metzger et al. (1984) for pressure loss predictions, both for the same setup. The study also differentiates between the heat transfer occurring by the body of the pin-fin itself and that by the end-wall areas surrounding it. Results from the CFD analyses based on the fourth pin-fin from the inlet (commonly assumed to have a stable flow around it), show very good prediction accuracies of heat transfer coefficients for the pin-fin body itself but rather low accuracies for the end-wall areas (based on heat flux and inlet temperature values). Better accuracies were obtained when using the enhanced wall treatment where pin-fin body heat transfer coefficients were almost identical between the CFD and experimental results. An alternative definition of heat transfer based on the averaged local temperatures around the fourth pin-fin showed that the heat transfer coefficient then (with CFD's capability to establish the local thermal field) is really between 1.5 to 3.5 times that predicted by using the inlet temperature in deducing the local h values. The same accuracies cannot be said about its predictions of pressure loss coefficients where CFD results tended to be lower by 50-100%h.


INTRODUCTION
Benchmarking efforts of CFD codes and models have occurred in the past with the most famous one being that of Freitas (1990).There, many CFD codes' results, each using their own meshes and solution schemes, attempted to solve 5 standard fluid flow cases.The results of this "competition" were devastating since huge differences occurred between the codes solutions even when the same turbulence model was used.A worse situation occurred when modeling flows around a square 180 o bend.There, all the CFD software failed to predict various features of the flow.The flow structures seen through velocity vector results showed completely opposite secondary flow directions of movements between the various codes.Further problems with numerical analyses were reported by Iaccarino (2001), where FLUENT, amongst other codes, showed an inability to correctly predict separation and re-attachment on a grid-independent mesh.The FLUENT analyses using the k-ε, DRSM and two-layer zonal models showed the impossibility of achieving a good comparison to experiment.Since no works could be found since then to document these issues, this study comes to further investigate the accuracy issues in the CFD code FLUENT coupled with a close look at its graphical regional results around the fourth pin inline to establish the flow development in the results provided.
The ANSYS-FLUENT CFD code: This is one of the world's first commercial and wide-spread CFD software used around the world.The code contains many physical treatments, turbulence models and near wall functions.It can solve flows in steady-state scenarios or in transient arrangements.The code also has various special treatments for species and radiation transport, fan and heat-exchanger models to mention a few.The number of turbulence models in version 13 of the code is more than 12 different approaches (depending on the way of looking at it).This is also added to 3 (or more other sub-options) for treating various issues such as the near wall regions.This obviously creates a wide range of solution approaches with little real guidance to know which to use in what case and what type of accuracies to expect from them.Therefore, a case with a complex geometry along with complex heat transfer and fluidflow physics was chosen to test the code's capabilities.

k-ε Turbulence model:
The most important turbulence model for steady state and 3D two equation model suggested by Jones and Launder (1972), which has gone through many enhancements since then.One such enhancement was to fix the flow attachment to the inner side after a 90 o bend (which is obviously in-correct).Other issues include problems with finding the correct flow field for problems with adverse pressure gradients.These problems were fixed in the fluent code around the year 2000.The model has two transport equations (partial differential type) of its own, namely k, the kinetic energy of turbulence and ε, the turbulent dissipation of this turbulence.The model is now the eldest and most used model in all codes and all over the world.This caused the k-ε model to be the model of choice in solving various problems since most other models were either new, too complex or simply could not converge easily.• Non-Equilibrium Wall Function (NEWF): This model was suggested by Kim and Choudhury (1995).The method sensitizes the fully turbulent treatments to be become flexible and extended to account for pressure gradients and variable properties by smoothly blending an enhanced turbulent wall law with the laminar wall law.The wall-neighboring cells are assumed to consist of a viscous sub-layer and a fully turbulent layer.The non-equilibrium wall functions model is recommended for use in complex flows involving separation, reattachment and impingement where the mean flow and turbulence are subjected to pressure gradients and rapid changes similar to the pin-fin case suggested here.By sensitizing the standard wall function model (above) to pressure gradients and by computing the budget of turbulence kinetic energy at the wall-adjacent cells, which is needed to solve the equation at the wallneighboring cells, the two-layer near wall turbulence is calculated using: a = 0.01 and b = 5.

Wall functions in
• Enhanced Wall Treatments (EWT): This method used to be called the Two Layer Zonal Method and was suggested by Kader (1981).
The method also used a one equation model of Wolfshtein (1969) in viscous layer near the wall where: The blending action between the viscosities of both layers uses: where, λ ε is the blending function (zero near wall and 1 away from it) and A is its width of the blending function.Chen and Patel (1988), give the length scales as: where, the quantity u c + is the value of u + at the fictitious "crossover" between the laminar and turbulent region.
Analytical/experimental treatments of pin-fins geometry: Pin-fins are essentially short tube banks, (Incropera and DeWitt, 1990).Applications of this geometry are also wide ranging from aircraft gas turbine blades cooling passages to using them in medium size heat exchangers, devices using pin fins can be as small as simple electronic heat sink devices or even inside computer-micro-channel electronic chips.To this end, it was necessary to find a good set of experimental data for the problem concerned to compare the CFD results obtained here to in the first place.Lyall (2006) provided a list of works published around such geometries.A close look at the literature showed the existence of two sets of data by different authors that can be trusted and to complement one another for heat transfer and pressure drop predictions.Such data were found in the works of Chyu et al. (1998) and that of Metzger et al. (1984), giving h and f values respectively for a common popular case of pinfin geometry at X/D = 2.5, S/D = 2.5 and H/D = 1, Fig. 1.An interesting observation was that both papers used air as the working fluid and their Re range was between 5,000 to 30,000.Following some data manipulation, the following formulae could be obtained (Chyu et al., 1998): Nu = 0.337 Re 0.585 Pr 0.4 Pin Nu = 0.315 Re 0.582 Pr 0.4  End-wall Nu = 0.320 Re 0.583 Pr 0.4 Combined (Metzger et al., 1984): f = 0.5244 Re -0.2For the geometry dimensions concerned.Although the Metzger et al. (1984) paper used all the same nondimensional dimensions for the geometry and Reynolds number used here, the pin-fin diameter used here is Metzger.This provides for four different criteria for the quantitative comparison of prediction accuracy between CFD and experiment.However, the combined Nu formula is of little concern here since its results should be already built into its two other sister formulae for the pin and end-wall regions.This leaves 3 comparison criteria to be used for this quantitative type of analyses.Therefore, the geometry chosen here has strong wall-flow interactions (such as the case in pin-fins as opposed to tube-banks), strong boundary layer creationdestruction sequences, wakes, separation and various curvatures (wall vs. circular pin-fins).It is therefore safe to expect that the length, velocity and time scales occurring in this problem are very complicated with high vortex shedding, eddies and shear stress production.Such a flow should give rise to relatively high levels of k, ε and to the production of k of any Reynolds Averaged Navier-Stokes (RANS) turbulence model.

METHODOLOGY
This study will attempt to compare the experimental data from Chyu and Metzgers' works mentioned earlier to the k-ε model CFD results presented here using three different near wall functions.A geometry created in a CAD package is meshed in the old Gambit1.1 and solved in FLUENT 13 (the latter two software are now part of the ANSYS suite of software).Results from FLUENT are processed using spreadsheet software such as Excel to come up with the comparisons.

The geometry modeled:
The geometry modeled is based on Metzger et al. (1984), setup with a diameter of 10 mm (Metzger et al., 1984) geometry diameter was 0.167inch or 4.2418 mm) and an Re = 5,000-30,000 where S/D = 2.5, X/D = 2.5 and H/D = 1, Fig. 1.Here, X is the longitudinal spacing, S is the side-wise spacing and H is the height of the pin-fin while D is the diameter which will be modeled dimensionally in the CFD code to be 10 mm and all other dimensions to be scaled from this. Figure 2 shows that five complete PFs are located in a single pass of staggered PFs that form a single complete passage for flow with half pins on both sides for symmetry.Hence, the boundary PFs are cut by two-parallel symmetry planes, therefore only 1/2 of their geometry is modeled here.Another symmetry plane passes through H/D of 0.5 to increase the number of cells placed along the height of the pin due to limited computational resources.Modeling of 5 PFs in a row allows for the use of the first three PFs for flow development and the latter 5th for exit (moving the exit away from the inside of the domain).This leaves the 4th PF to extract more accurate data (through good flow development and accurate surroundings) that can be applied to any other geometry with a higher number of pin-fins accurately.

Mesh used:
The mesh used in these analyses is a tetrahedral element mesh created previously in ANSYS Gambit 1.1, Fig. 3.By placing 10 tetrahedral elements across the length of half a PF (modeling the other half effects using symmetry boundary condition across half of the H/D length to save on mesh elements), this effectively provides a 20 tetrahedral elements across the PFs complete H/D length (a higher cell count along the pin's length will be achieved through mesh adaption -Fig.4).The initial uniform will then provide a total tetrahedral mesh cell count of 1,224,536 in half the domain modeled.The maximum skewness of this mesh is 0.5617 (the most important FLUENT mesh quality measure) or a minimum quality of 0.1444, rendering the mesh as a very good initial mesh (will be solutionadapted later) by both industrial and academic meshing standards.

Mesh adaption:
For industrial applicability of this study purposes and due to available computational memory issues, no grid independence exercise was carried out here.Nevertheless, a grid adaption exercise was carried out to the converged k-ε solution at Re = 5,000 to reduce Y+ to 5 and to halve all the important gradients of total pressure, velocity magnitude, static temperature present in the domain.This resulted in the original mesh cell count of 1,224,536 to rise by 399,399 showing an in-between measurements section also Fig. 5: Final adapted surface mesh around the inlet, one side's symmetry plane and pin number 1 (centre void) until pin 1.5 to 1,623,935 tetra elements where the same resulting mesh was used for results throughout all the analyses here including the Re = 5000 originally condition used for adaption in the first place.Figure 4 and 5 shows the effects of the grid adaption process where cells along the end walls and around the sides of the pins have been refined more than once at various places.The adapted cells are also located in areas where high velocity and pressure (and their gradients) tend to take place.

The operating fluid and boundary conditions:
The operating fluid will be air with the following properties used in the FLUENT code is: T wall = 293 /K All initial flow conditions were obtained using inlet values, including suggested k and ε that FLUENT predicts from the turbulence intensity and hydraulic diameter using analytical formulae.
By comparing the CFD results obtained here to those obtained experimentally by Chyu et al. (1998) and Metzger et al. (1984), it is hoped that a better understanding will be reached of the physical problem and a proper assessment of FLUENT's prediction of heat transfer and pressure drop capabilities be reached.Since fluent comes with a number of near wall functions and turbulence models, this study is the first of a number of works attempting to find the better setup to be used in fluent for modeling enhanced heat transfer and pressure drop problems such as the one modeled here containing boundary layer separation and reattachment.

RESULTS AND DISCUSSION
Solver performance: Figure 6 shows the residuals' convergence requirements for all three near wall treatments.The figure shows the effect of Re on the whole range of Re studied is not within the turbulent flow region, nor do transition effects take place here as one would be tempted to assume.The plausible explanation for this may be found in the change of Y+ for near wall cells caused by Re increases wake attenuation effects on the SWF.It would have been also interesting to see if this behavior continued at Re>30,000, or if a sinusoidal behavior may occur, however, this is beyond the scope of this study.
The NEWF solutions on the other hand seems to independent of Re possibly due to the method's blending of laminar and turbulent near wall functions simultaneously leading to a steady behavior with respect to flow Re.This is a comforting behavior but may become unstable at certain conditions from a mathematical point of view.Results for EWT method solutions require less iteration with increasing re.This is a welcomed finding as it gives hope to solving high Re cases with little computational effort.However, this result is also thought to be the result of the correction and blending functions within the EWT method and may not be a reflection of an accurate solution due to the need for higher computational effort needed to solve higher Re problems as given by the Kolmogorov theorems.A better clarification of the better near wall method for solution may become clearer when considering the actual benchmarking of h and f results later.
Flow and pressure fields: Pin-fin CFD results express most of the normal features of flow over tube-banks including pin-fin upstream stagnation points, downstream re-circulation wakes and sides accelerating separation flows, Fig. 7. Velocity magnitude contours also show the flow development at the fourth pin when observing the flows around its neighboring half pins at the symmetrical borders of the geometry.Figure 8 shows the creation of a stagnation (static) pressure region at the frontal impingement area of the pins.
Figure 9 shows the distribution of turbulence around pin4 using turbulence intensity % results.The highest areas of turbulence intensity are areas of highest shear and production of k areas where Ti values can reach 250% at upstream pin locations close to maximum velocity areas.However, relatively mid-level turbulence areas of 100% turbulence and above occur in the impingement areas upstream of the pins and in the wake areas downstream of the pins.The Ti levels immediately behind the pins are themselves not low in reality and stand at 45% and above as can be seen from the figure.
Effective viscosity contour plot shown in Fig. 10 show the trend in which viscous resistance is added to the flow domain above the molecular levels of viscosity observed at the lowest level in colormap distribution seen to the left of the figure (dark blue at 1.789*10 5 ).Areas of high viscosity are located within the wake downstream of the pins providing resistance to flow penetration from pins' sides.Furthermore, wall areas with high velocity and velocity gradients are sources for the production of turbulent energy, Fig. 11.This figure shows that the only areas of any noticeable production of k are really limited to these areas only.This production of turbulent energy is then carried out throughout the domain and is budgeted to the dissipation of k using k and ε transport equations of the model.
It would be easy to become mixed up with Fig. 11 and 12 since they look exactly the same.This is a mere expression of physical essence meaning that the areas of highest production of k are the same areas where ε is the highest.Considering the physics involved, this is very true since the very high shear levels (relatively speaking when considering other fluids) present in these areas, Fig. 13, will not only lead to the generation of turbulent energy but to the excessive "wear" of this energy due to the same cause being shear.This in turn produces the turbulence and its intensities seen here and elsewhere in Fig. 14.
Thermal and heat transfer effects: The latter picture of flow field effects the levels of heat transfer especially near the wall where the thermal field and the thermal and turbulent effective Prandtl number dictates the amount of heat conducted across the boundary layer present and then convected away from the walls by the passing flow streams.Figure 15 shows the contours of static temperature field surrounding pin 4 showing the convective nature of the problem where temperatures follow the main stream of flow convection.Figure 16 shows the effective Pr field at the symmetry mid-plane indicating that higher Pr levels are upstream and downstream of the pin bodies and closer to its molecular value of 0.744 near flow separation points.Heat flux from pin 4 and the end-walls surrounding it can be seen in Fig. 17 where the end-wall seems to be participating heavily in the heat transfer as opposed to the pin-fin with the lower surface area and many low heat flux pockets.This is an interesting finding since another look at the wall shear results in Fig. 18 show that the end-walls impose much lower wall shear on the flow.While not undermining the role of the pin-fin in augmenting heat transfer levels offered by the geometrical arrangement, a clear indication appears that it there could be scenarios where a low packing density of pin-fins maybe better for the heat transfer versus pumping power losses than a high density pin-fin pack.Near wall treatments: Figure 7 shows the heat transfer at pin 4 alone using the three near wall functions.Since we have no knowledge of an analytical/experimental heat exchange here, we will only use this to indicate that the heat transfer predicted with the SWF and NEWF methods seems to be identical while the heat transfer levels predicted by the EWT method at Re >15,000 is above that predicted by the other methods.Since the NEWF seems to under-predict h by a near constant 7-8%, except below Re = 10,000 (where a higher difference occurs), it is possible to apply correction factors here.This could be the same level of correction factor needed to be applied to the Pr level used by the code in its h calculations.A troublesome picture in heat transfer predictions begins to appear when looking at the h levels predicted at the end walls seen in Fig. 9. Here, not only do the results not match (except around Re = 20,000), the trends of increase in h values with increasing Re are completely wrong but also flip around from under- predicting to over-predicting.This is caused by a higher than 1 power index between h and Re when the experimental h values (solid line and black shaded coefficients) seen in Fig. 9 cause a non-linear difference between the two sets of data.The discord between the CFD and experimental results at Re = 5,000 reaches 200%.The huge discrepancy in the power indices indicates a completely physically wrong method in calculating the h values there.It is actually inconceivable that the power indices here are above 1 since heat transfer literature does not show such a power index for any forced convection problem where all power indices have tended to be below 0.8 for example, note the Dittus-Boelter formula (Incropera and DeWitt, 1990).The huge overall discrepancies can be seen again using h EXP /h CFD ratios in Fig. 10 and 11 for the pins and end walls respectively.This is an issue that the code vendors should consider fixing due to its physical significance.Figure 19 shows heat transfer levels at pin 4. While heat transfer predictions for the EWT method is generally higher than the SWF and NEWF results, it is only Fig. 20 that shows the relation between hpin for pin4 to that predicted by CFD where it is clear again that the EWT method is overly predicting the heat transfer while both the SWF and NEWF results are close to one another but lower than the analytical value (solid line).With the agreement between SWF and NEWF methods, the continuous under-prediction and similarity in curvature of these two methods results to that predicted by the analytical methods based on experimentation, it is plausible to think and employ "correction" factors or to recommend a modification to the near wall thermal Prandtl number to fix this.However, this cannot be done since the end wall region's heat transfer coefficient values vary considerably and in a non-orderly manner, Fig. 21. Figure 22 and 23 express the complicated errors in h ratio predictions for pin4 and its end wall regions respectively.The CFD predicted heat transfer using whichever wall function for the pin and end-wall regions is much lower than the experimental (analytical) at Re <10,000 regions.Special attention should be given to these areas due to this.

CONCLUSION
This study benchmarks the performance of various ANSYS FLUENT's k-epsilon turbulence model near wall functions with results showing that the Non-Equilibrium Wall Function being the best in capturing convective heat transfer results.CFD results obtained here showed an under-prediction of heat transfer coefficient values at Re <10,000 with the inconsistencies between such predictions for the pin and for the end-wall regions preventing a unified correction factor or a change to the Pr thermal to be applied here.Also, all wall function methods under predicted the pressure drop by a factor of 2.
fluent: The wall functions serve CFD codes by giving the shear and turbulence boundary condition situation at the walls.Fluent has three such function treatments for its k-ε model given in Fluent's users guide manual (ANSYS Inc. version 12) as: • Standard Wall Function (SWF): The standard wall function came out of the works of Launder and Spalding (1974).In the near wall area: If flow near wall is laminar (Y* < 11.125) then: U* = y* If flow near wall is turbulent (Y* > 11.125) then: 4187 (von-Karman constant), E = 9.793, U p = mean velocity of the fluid at the near-wall node, k p = turbulence kinetic energy at the near-wall node, y p = distance from point to the wall, dynamic viscosity of the fluid kp= turbulent kinetic energy at the first nearwall node P, C p = specific heat of fluid, T p = temperature at the first near-wall node P, T w = temperature at the wall, Pr = molecular Prandtl number, Gk = production of k and ε = dissipation of turbulent kinetic energy.The model assumed Pr t = turbulent Prandtl number to equal 0.85 at the wall.

Fig. 6 :
Fig. 6: Convergence requirement for pin-fin CFD analyses using the k-e turbulence model using 3 different near wall treatments

Fig. 13 :
Fig. 13: Wall shear stress contours using the k-ε turbulence model and standard wall functions around pin 4 [pa]

Fig. 15 :
Fig. 15: Static temperature contours using the k-ε turbulence model and standard wall functions at pin 4 symmetry mid-plane [-]

Fig. 17 :
Fig. 17: Wall heat flux contours using the k-ε turbulence model and standard wall functions at pin 4 symmetry mid-planes [W/m2]

Fig. 20 :
Fig. 20: Heat transfer coefficients for pin 4 versus Reynolds number showing Chyu's experimental results with solid line along with k-e predictions using SWF, NEWF and EWT Fig. 22: Heat transfer coefficient ratio for pin 4 (Exp/CFD)versus Reynolds number with CFD results using the k-e predictions using SWF, NEWF and EWT