Effects of Turbulence Models and Grid Densities on Computational Accuracy of Flows Over a Vertical Axis Wind Turbine

Flows through a vertical axis wind turbine (VAWT) are very complex due to their inherent unsteadiness caused by large variations of the angle of attacks as the turbine is rotating and changing its azimuth angles simultaneously. In addition, a turbine must go through a wide range of operating conditions especially the change in blade speed ratio (BSR). Accurate prediction of flows over VAWT using Reynolds-Averaged Navier-Stokes (RANS) model needs a well-tested turbulence model as well as a careful grid control around the airfoil. This paper aimed to compare various turbulence models and seek the most accurate one. Furthermore, grid convergence was studied using the Roache method to determine the sufficient number of grid elements around the blade section. The three-dimensional grid was generated by extrution from the two-dimensional grid along with the appropriate y+ controlling. Comparisons were made among the three turbulence models that are widely used namely: the RNG model, the shear stress transport k-ω model (SST) and the Menter’s shear stress transport k-ω model (transition SST). Results obtained clearly showed that turbulence models significantly affected computational accuracy. The SST turbulence model showed best agreement with reported experimental data at BSR lower than 2.35, while the transition SST model showed better results when BSR is higher than 2.35. In addition, grid extruding technique with y+ control could reduce total grid requirement while maintaining acceptable prediction accuracy.


Introduction
The flow through vertical axis wind turbine has been affected primarily by 3 factors, the viscous effect, the impact change of angle of attack and the blade speed ratio (BSR) (Paraschivoiu 2002).Theses make the flow at a different BSR appear to represent a different regime.To get the accurate computational fluid dynamics (CFD) result using the Raynolds Average Navier-Stoke (RANS) model, attention should be paid to the grid density around the airfoil and turbulence model used.
As a CFD progression, several turbulence models has been developed.Many models were used to simulate a flow through vertical axis wind turbine such as Spallart allmaras (S-A) model, Renormalization group (RNG) model, k-ε model, Shear stress transport (SST) k-ω model and Transition SST model.Spallart Allmaras (S-A) model (Spalart and Allmaras 1992), is the one equation model which based upon turbulent kinematic viscosity transfer equation.The model was designed for aerodynamic problems, especially for the wall bounded flow and the flow involving an adverse pressure gradient.It was designed as a low Reynolds number model, which needs to calculate the viscosity effect in the boundary layer.The gradient of transport variables in the boundary layer is relatively low compared with the k-ε and k-ω models, this makes it less sensitive to a numerical error.The Spalart-Allmaras model has an acceptable accuracy for a twodimensional simulation and takes a relatively short calculation time compared to two-equations models.
The k-ε turbulence model (Jones and Launder 1972) is a semi-empirical model based on the turbulence kinetic energy (k) and turbulence dissipation rate (ε) equations.Assumptions of the model are the flow is fully turbulent and molecular viscosity is neglected.These make it unable to predict separation and turbulent behavior of flow in the boundary layer accurately (Wilcox 1993).The Renormalization Group (RNG) k-ε model (Orzag et al 1993) is the extension of the k-ε model which uses a statistical technique called the Renormalization group theory.The model accounts for the impact of the swirling of fluid (eddy) and includes the effect of viscosity in the low The k-ω model (Wilcox 1998) is a two-equation model modified for a low Reynolds number on the basis of turbulence kinetic energy and specific dissipation rate transfer equations.It gives the more accurate result in the boundary layer region than the k-ε model and is successful in predicting the moderate adverse pressure gradient problem (Menter 1993).Unfortunately, this model fails in predicting the flow with pressure induced separation.Moreover, the ω equation is sensitive to its own value in the free stream outside the boundary layer (Menter et al 2003).This led to the development of the Shear Stress Transport model (SST) which aims to solve the separation and relatively high backward pressure gradient flow (Menter 1993, Menter 1994).
The SST k-ω model is the model that blends the k-ε and k-ω models together by using the blending function.The blending function is the unity at the wall surface and the k-ω model is applied.For the effect of the far-field boundary conditions on the boundary layer, the function is equal to zero and the k-ε model is applied instead.Such a combination allows this model to solve the flow, both in the boundary layer and the fully turbulent flow outside.
The Transition SST model is the extended SST model that has two additional equations, the intermittency and the transition onset criteria equations, which are in the form of the momentum thickness Reynolds number to cover the effect of the bypass transition and low freestream turbulent.Both equations are developed by (Menter 1994), the purpose is to cover the bypass transition and low free-stream turbulent.The Transition SST model, therefore, predicts transition flow better than the SST model.
The wall is the main factor causing eddies or turbulence.It is clear that the mean velocity field is affected by no-slip conditions caused by the surface of the wall.Near-wall turbulence can be classified into 3 layers: 1) The viscous sublayer, the most adjacent layer to the wall.In this layer, the flow is almost laminar and viscosity plays an important role in the momentum as well as mass and heat transfers; 2) The outer layer called the fullyturbulent layer which turbulent plays an important role; and 3) The buffer layer or blending region which is the middle layer that is affected by both viscosity and turbulence.In general, there are two approaches to modeling the boundary layer flow.In the first approach, the viscous sub-layer which is directly affected by the viscosity, and the buffer layer is not resolved but uses a semi-empirical called "Wall functions" to connect the effect of the viscosity between the viscous sub-layer and fully turbulent region.The second approach relies on turbulence modification to calculate the effect of viscosity on the viscous sublayer.Therefore, it is necessary to generate enough grids in the viscous in order to compute the viscosity effect in this area.This method is called the near-wall modeling approach.The k-ε group turbulence model is not designed to calculate the impact of viscous in the viscous sublayer.The simulation needs to use the wall function.In this case, the centroid of the first cell adjacent to the wall should located in the fully turbulent region or log-law layer, the y + value is about 30-100 which can be estimated from equation 1.
For The S-A, SST and transition SST k-ω model, they were designed to calculate in the viscous sublayer and buffer layer.The centroid of the first grid is usually set at y+ ≈ 1.Although there are many turbulence models available, they are usually designed to solve the problem at some specific flow regime.The model that is very accurate at one BSR may not be accurate in other ranges.Hence, there are no turbulence models that can predict the turbine efficiency accurately for all operating speeds of VAWT.There have been many comparative studies on turbulence models conducted but were usually done at a specific BSR.Nobile et al (2013) simulated a flow through an augmented vertical axis wind turbine.Three turbulence models, k-ε, SST k-ω and transition SST k-ω model were chosen and compared, the resulting trend looked reasonable but the values were quite different from the experimental result.Nobile et al (2014) made a comparative CFD study of upright and tilted VAWT and compared three following turbulent models, Spallart-Allmaras model, RNG and SST model at one BSR and found that the SST model gave the best result.In the same year, Almohammadi et al (2015).study a dynamics stall of H-rotor by using SST and SST transition model and point out that Transition SST gave the better solution in observing the dynamics stall behaviours.
Besides the effect of turbulence model that is not clearly depicted, the sufficient number of the grid elements has not been clearly identified.This paper aims to illustrate the grid convergence study to achieve a sufficient grid number around the airfoil and then seek the most accurate turbulence model throughout the range of BSR by focusing on the turbine efficiency.

Grid convergence study
Since the study of grid size in three-dimensional model is more difficult than two-dimensional studies.(Almohammadi et al 2015).Using two-dimensional study methods to guide the creation of a grid in three dimensions makes it easier to manage and control the grid resolution (Paraschivoiu et al 2014).In order to investigate the number of sufficient grid elements around the blade section, the Roache method (Roche 1998), the widely accepted and used in the assessment of numerical uncertainty which based upon Richardson extrapolation theory [Richardson extrapolation] was applied to study the grid convergence in the two-dimensional domain.The quasi-2-Dimensional experiments (Oler et al 1983) were simulated using the two-dimensional grid.The rotor had only one blade with a rotor diameter of 0.61 meters.The airfoil section was NACA0015 with 0.1524 m. of chord length.The turbine was tested in the water tunnel of 5 x 10 x 1.25 m 3 .The rotor speed was controlled to rotate at a constant speed of 0.74918 rad/s and the freestream velocity was adjusted to get the design blade speed ratio.In this case, the operation speed BSR of 2.5 was selected.
The configuration of the tested turbine is shown in Figure 1.Experiment tangential and normal forces were measured using strain gauge that mounted on the handle rod.To satisfy the two-dimensional flow, the tip effect was corrected using the Graham method (Graham 1982) so the tested data from this experiment can be compared to data in the two-dimensional domain.The domain consists of two subdomains, the rotating and stationary domains as shown in Figure 2.During the calculation process, the rotating domain rotated with the rotor speed while the stationary domain was fixed.The flux that is moving through the interface surface between two domains was calculated by an interpolation method.The grid around the airfoil is constructed as an O-grid.
Firs grid centroid was located at 0.014 mm (y+ = 1) and grew with the growth rate of 1.05 before it was adjusted to fit the surrounding area.Figures 3, 4, and 5, show the grid geometry around the airfoil, grid connection around the rotor and grid layout through the domain, respectively.The Spalart-Allmaras model was selected to use as the turbulence model used in this study.The turbulent intensity was set to 5% and the turbulence length scale was 0.07 of the chord length.The boundary condition at the outlet was set to atmospheric pressure.The wall of water tunnel and blade surfaces were set to no-slip wall.The simulation was performed using a segregated pressure-based solver.The pressure-velocity coupling scheme for the solution method was SIMPLE.Spatial discretization for the gradient was the Green-Gauss cellbased scheme.Convergence criteria for RMS was set to 1×10 -6 .The calculation started with the first-order upwind scheme and switched to the second order upwind scheme to achieve accuracy.During the calculation, the tangential force coefficient was monitored, and the calculation stopped when the periodic solution was achieved.
Tangential force coefficients of the blade were defined as While Ft was the force exerted on the airfoil section in the direction tangent to the blade path.r was water density, c was chord length and ¥ U was the freestream velocity.
Rotor torque and power can be found from equations Power coefficient was defined as (Cp) The calculated power coefficients related to the grid resolution are shown in Table 1 In this regard, a threegrid set with a different number of elements around the blade section, fine, medium, and a coarse grid have been created with grid refinement about 1.2.Average grid size can be calculated from Where ∆ > is a cell area of cell i, N is a total element number, and grid refinement factor is All three sets of grids must be on the condition that ℎ @ < ℎ B < ℎ C (index 1 represent the finest grid, 2 and 3 were a medium and coarse grid, respectively.In this case, the grid refinement factor was defined as where ,  B@ was a grid refinement factor of grid 2 and grid 1 and  CB was a grid refinement factor of grid 3 and grid 2. For the straight blade VAWT, the main parameter to use for grid convergence assessment was the power coefficient (Cp).The appearance order, p can be calculated iteratively.
( ) Finally, the convergence criteria were calculated from means Oscillatory divergence

Turbulent model effect study 2.2.1 CFD model
The three straight blades VAWT experimental set up by (Howell et al 2010) were reproduced numerically.The turbine was tested in a low-speed wind tunnel at 1.2 m × 1.2 m square section and 3.0 m of working section length.The turbine blade section was an NACA0022 profile with a 100 mm of chord.The rotor radius of 300 mm and the blade height of 400 mm resulted in turbine solidity and blade aspect ratio of 1.0 and 4 respectively.The blade geometry is shown in Figure 6.For simplicity, computational domains were created by excluding the turbine axis and radial arms.The domain was divided into two parts: stationary and rotational.
Assuming the flow was vertical symmetry, each part was created with half length of the total height in order to save computational time and resources.A total of 190 cells was created around the airfoil in the two-dimensional model and extruded to a three-dimensional model.

Near-wall treatment
For SST and transition SST models, which are the low Reynolds number type models., the first centroid height was generated with y+ of unity for an enhance wall model.For RNG k-e model, wall function was implemented.Note that the "wall function" is semi-empirical formulas used to bridge the viscosity-affected region between the wall and the fully turbulent region.In this case, the viscosityaffected inner region (viscous sublayer and buffer layer) is not resolved.Each wall-adjacent cell's centroid, therefore, should be located within the log-law layer, 30 < y+ < 300.So, the grid was generated using y+ =30 with RNG k-e model.The mesh grew with a rate of 1.12 before merging with the outside mesh layer.Grids were first generated in the 2-dimentional domain and extrude to be 3-dimensional domain along the blade span.The approximate total mesh number was 8 × 10 6 elements.Computational domains and grid are shown in Figures 7 and 8 respectively.

Simulation setup
Boundary conditions, as shown in Figure 9, were set the same as the experiment.The inlet velocity was 5.07 m/s.The outlet was set to outflow with the pressure of 101325 Pa.Blade surface and the rest of the 4 sides of the wind tunnel wall were set to be no-slip walls.The turbulence intensity and length scale were set to 0.01 and 0.1C (chord length), respectively.The rotor speed was varied to get different blade speed ratios.To ensure that the computed results were reasonable, solutions were initially obtained using first order spatial discretization and later switched to the more accurate second order until periodic solutions were achieved.During the calculation process, power and torque coefficients including residuals of each transport variable were monitored.The design target residuals were 10 -8 and the calculations were stopped when periodic convergences were achieved which usually took about 4-5 rotor revolutions.

Grid convergence study
The result of grid convergence study using a two dimensional one blade turbine and grid system in fig.3, 4 and 5 is shown in Table 2.The calculated convergence criterion was 0.1487.The value was between 0 and 1, that means the monotonic convergence was achieved.Appearance order (p) appears to be 6.8726 and the asymptotic power coefficient (  F_HI8>JK8HL ) was equal to 0.1325 or 13.25% while the experiment power coefficient was 0.1082 or 10.82%.That means the CFD solution deviated from the experimental value of 22.47%.The asymptotic power coefficient of 13.25

Stationary domain
Rotational domain means if the grid was kept refined until the solution was constant, the solution of power coefficient will converge to the value of 13.25.Power coefficients calculated from each grid set are shown in Figure 10.The calculated tangential force coefficient from three grid sets are shown in Figure 11.The calculated solution from all sets of the grids from 0 o to 90 o of an azimuth angle exhibited over the predicted value.The solution obtained from the coarse grid was explicitly different from the tested data while solutions from medium and fine grids were more reasonable and closer to the test data than the coarse grid.The maximum torque occurred at about 80 o .After 90 o , all solution sets tend to agree well with the tested data except the results between 250 o and 290 o were less than predicted.
The streamline and vorticity plot show the reasonable rotor wake.In Fig. 12 , the streamline expand after the the flow going through the blade passing path.This result in increasing of velocity in the vortex core area.Fig. 13 show that vortex were shedding and get smaller as it move to the downstream.For a conclusion of grid convergence, after the grid was refined two times, monotonic convergence was achieved.The finest grid with 220 elements and the medium grid with 160 elements around the airfoil gave results that agree well and quite close to the tested data.In this case, using 190 cells around the airfoil may not give too many different solutions but save calculation time.The grid manipulation method above was considered successful because it gave an accurate result with a relatively small amount of total grid elements.

Turbulence model effects study
The result of turbulence model study using the three dimensional grid set of three vertical blade turbine in fig.6,7 and 8.
Figures 14 and 15 show the measured and computed values of cp and Torque as functions of BSR.The power coefficient (cp) was defined as power extracted from the wind divided by the power available in the wind at the same frontal area.Considering the streamline around the turbine blade from Figure 14 to Figure.17, both BSR 2.15 and 2.5 show a similar flow pattern, but BSR 2.15 clearly demonstrates a stronger separation.In the azimuth angle of about 0 to 30, the flow is still attached, no separation.After that, a separation bubble begins to appear on the trailing edge, then expands wider until throughout the chord length at about 150 degrees of the azimuth angle.This occurrence of separation represents the trailing edge stall (Mcculough and Donald 1951), which lift force is continuously reduced as the separation bubble expands to the leading edge but does not drop suddenly like a leading-edge stall.At around 210° vortex is shedding, and the flow returns to reattach.For BSR 2.5, the flow has a similar pattern, but less intense.
The high BSR means the higher rotor speed or lower wind velocity which the flow may be in transition regime.This lets Transition SST be more accurately predicted.
However, it is evident that the SST model yields consistent and close approximation to the values obtained from the experiment throughout the range of BSR from about 1.85 to 2.6 while the Transition SST model gave an underpredicted in BSR range between 1.84 and 2.25 which is a period when strong separation occurred.Therefore, in the case of predicting the overall turbine efficiency without the need for great precision, SST can be a good choice to simulate the flow throughout BSR because of less calculation time was higher than 2.15 onwards, so it provides low accuracy.This may be because the model is not designed to calculate in the viscous sublayer and buffer region, but used the wall function model instead.The wall function is based on the formula achieved from the experiment which conditions are different from VAWT condition.Figures 18  and 19 show the vorticity plot around the rotor compared between the RNG and SST models.The result agrees well with each other, and the flow pattern is reasonable and has a similar pattern.The result indicates that BSR affects the wake length which agrees with Fujisawa et al (2001).The wake shape is not symmetric.At a low BSR, the wake length is long, while it is short at a high BSR.They packed together as a group before shed to the downstream.At the BSR lower than unity, in the azimuth angle from 150° to 210°, the wakes are found to shed backward to the leading edge by the force of freestream.This lets the flow strongly separate due to a high adverse pressure gradient.On the other hand, at the BSR higher than or equal to unity, the wake vortices are shed to the trailing edge because the blade speed is higher than the wind speed.However, in the downstream, the blade hit it own wake, resulting in a drop of lift force and low efficiency

Conclusion
Grid convergence was successfully achieved by using the Roache method.The result indicates that 220 cells around an airfoil with appropriate y+ control is adequate to calculate VAWT efficiency.Based on this number of the grid elements, and y+ control, the method to generate a three-dimensional grid by extruding from a two-dimensional grid is easy to manipulate.The total grid number achieved from this methodology is quite small but can be predicted accurately.The turbulence effect study agrees well with the experiment data.The Transition SST k-w turbulence model gave the most agreeable results at BSR beyond 2.35, below this the SST model gave the best accurate result.It can be seen that turbulence models used have a significant impact on the accuracy of the numerical solutions in VAWT, 20-50% departures from experimental values can be noticed.

Fig 1 .
Fig 1. Oler's one blade VAWT tested in the water tunnel

Fig 2 .
Fig 2. Domain used in calculation

Fig 6 .
Fig 6.Blade geometry The turbine and wind tunnel configuration is shown below

Fig 7 .
Fig 7. Grid and computational domain for the 3 blades VAWT

Fig 10 .
Fig 10.Calculated power coefficients versus grid set

Fig 17 .
Fig 17.Streamline around a middle span for BSR of 2.15 and 2.5, at an azimuth angle of 0° to 60 °For the RNG model, the accuracy of the solution in the range of BSR from 2.15 to 1.8 is quite good, but the BSR

Fig 18 .
Fig 18. Streamline around a middle span blade section for BSR of 2.15 and 2.5, at an azimuth angle of 90° to 150 °

Fig 20 .
Fig 20.Streamline around a middle span blade sectionfor BSR of 2.15 and 2.5, at an azimuth angle of 90° to 150 °

Table 1
Turbine and wind tunnel configuration.