On the accuracy of turbulence models for CFD simulations of vertical axis wind turbines

A comparative analysis of seven commonly-used eddy-viscosity turbulence models for CFD simulations of VAWTs is presented. The models include oneto four-equations, namely the Spalart-Allmaras (SA), RNG k-ε, realizable k-ε, SST k-u, SST k-u with an additional intermittency transition model (SSTI), k-kl-u and transition SST (TSST) k-u models. In addition, the inviscid modeling is included in the comparison. The evaluation is based on validation with three sets of experiments for three VAWTs with different geometrical characteristics operating in a wide range of operational conditions, from dynamic stall to optimal regime and to highly-rotational flow regime. The focus is on the turbine wake, the turbine power performance, and the blade aerodynamics. High-fidelity incompressible unsteady Reynolds-Averaged Navier-Stokes (URANS) simulations are employed. The extensive analysis reveals high sensitivity of the simulation results to the turbulence model. This is especially the case for the turbine power coefficient CP. The results show that the inviscid, SA, RNG k-ε, realizable k-ε and k-kl-u models clearly fail in reproducing the aerodynamic performance of VAWTs. Only the SST model variants (SST k-u, SSTI and TSST) are capable of exhibiting reasonable agreement with all the experimental data sets, where the transitional SST k-u versions (SSTI and TSST) are recommended as the models of choice especially in the transitional flow regime. © 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Two promising and rather unexploited locations for further wind energy harvesting are deep-water offshore and urban/suburban environments [1]. Concerning the former, floating wind farms are currently being investigated [2e8]. Concerning the latter, the potential locations are on the building roofs [9e16], integrated within/between the buildings [17,18] and inside the ventilation systems [19e25].
For the aforementioned locations, vertical axis wind turbines (VAWTs) are highly advantageous over horizontal axis wind turbines (HAWTs) mainly due to their omni-directionality, scalability, less manufacturing/installation/maintenance costs, less noise, higher space efficiency, less shadow flickering, bird/bat safety, and more visual aesthetics [26e29]. However, while several studies have been performed on the impact of various operational [30e33] and geometrical parameters of VAWTs [34e38] and to improve their power performance [39e42], the amount of research is still incomparable to HAWTs. Therefore, the aerodynamic performance of VAWTs in terms of power coefficient C P is currently lower that of HAWTs [43e49]. Note that the aerodynamics of VAWTs is highly complex encompassing simultaneous large excursions of the angle of attack and the relative velocity, which result in several flow features such as dynamic stall, blade-wake interactions, rotational flow effects and flow curvature effects [50e54]. In order to benefit from the advantages of VAWTs for the aforementioned locations, therefore, their aerodynamic performance needs to be further improved.
Due to the high computational cost of the scale-resolving simulations, such as LES [52,55] and hybrid RANS/LES [56e59], the vast majority of the existing numerical studies on VAWTs is based on the unsteady Reynolds-averaged Navier Stokes (URANS) simulations [60,61]. The accuracy of the URANS simulations, however, strongly depends on the computational settings. Previous studies have shown the important impact of the domain size, time step (in terms of azimuthal increment defined as the number of degrees that the turbine rotates per time step) and convergence criterion [62e64]. However, to the best of our knowledge, no systematic investigation has been performed to critically assess the accuracy of the turbulence models in comparison against experiments in order to identify the best-performing models for CFD simulations of VAWTs. This indicates the need for more extensive sensitivity analyses to support future CFD studies of VAWTs.
A detailed review of the literature of CFD studies on VAWTs reveals that while few studies have been performed to provide insight into the impact of turbulence models [65e71], the findings are not conclusive. This is because: (i) the studies analyzed only a few turbulence models; (ii) the predictions of turbulence models were only compared against each other rather than being validated against experiments; and (iii) the focus was mainly on the turbine power coefficient C P rather than the details of the turbine instantaneous loads and moments, turbine wake and blade aerodynamics. Note that as C P is averaged over one or more revolutions, error cancellation could impede a fundamental analysis of the predictions of different models.
Therefore, in the present study: (i) A comprehensive analysis is performed to compare the accuracy of different eddy-viscosity turbulence models for URANS simulations of VAWTs, where seven one-to fourequation eddy-viscosity models together with the inviscid modeling are employed. The models are the Spalart-Allmaras model, RNG k-ε model, realizable k-ε model, SST k-u model, SST k-u with an additional intermittency transition model, kk l -u transition model and g-Re q transition SST k-u model.
(ii) The evaluation is based on validation with three sets of windtunnel measurement by Ferreira et al. [51], Tescione et al. [72] and Castelli et al. [73]. The turbines in the experiments have dissimilar geometrical and operational characteristics.
The best-performing turbulence model is selected as the model which provides the best agreement with all the three sets of measurements. Note that this should be done only after grid convergence analysis has shown that numerical (discretization) errors remain sufficiently low. Therefore, a grid convergence analysis is also included in the present paper. (iii) The analysis focuses on different parameters, namely the turbine wake, the turbine power performance, and the blade aerodynamics.
The wide diversity of the studied turbulence models, the characteristics of the turbines in the three experiments and the analyzed parameters help to scrutinize the performance of various turbulence models and support the generalization of the findings of the present study for accurate CFD simulations of VAWTs. Section 2 presents the description of the three wind turbines and the wind-tunnel experiments. The computational settings and validation studies are mentioned in Section 3. Section 4 presents the details of the employed turbulence models. Analysis of the predictions of the different turbulence models for the turbine wake, the power performance and the blade aerodynamics are discussed in Section 5. Discussion and conclusions are provided in Section 6 and 7, respectively.

Turbine characteristics
In this study, three Darrieus H-type VAWTs with different geometrical and operational characteristics are employed. Table 1 presents the main characteristics of the three turbines. A schematic of the turbines is shown in Fig. 1

Wind-tunnel experiments
Turbine 1: The experiment by Ferreira et al. [51] was performed in the low-turbulence wind tunnel (LTT) at TU Delft, where planar particle image velocimetry (PIV) was employed to measure the flow velocity in the rotor midplane of a one-bladed turbine operating in dynamic stall. The measurement revealed the details of the evolution of the separated/shed leading-edge/ trailing-edge vortex and the corresponding circulation strength.
The tips of the blade were covered with circular plates to avoid the 3D tip effects. The results were phase-averaged over few dozens of turbine revolutions, however the exact value is not reported. Turbine 2: The experiment by Tescione et al. [72] was performed in the open-jet facility (OJF) at TU Delft. Stereoscopic PIV was employed to measure the time-averaged streamwise and lateral velocity components at the rotor midplane in the wake of a twobladed turbine operating in the optimal regime. The high blade aspect ratio, i.e. blade height/chord length ¼ 16.67, was selected to minimize the 3D tip effects at the rotor midplane. The results were phase-averaged over 150 turbine revolutions. Turbine 3: The experiment by Castelli et al. [73] was performed in the low-turbulence wind tunnel at Politecnico di Milano. Low-frequency, i.e. 0.5 Hz, torque and rotor speed measurements were performed to measure the turbine power coefficient for a three-bladed turbine operating within a wide range of tip speed ratios 1.44 l 3.29. The range covered three operational regimes, namely dynamic stall, optimal operation and highly-rotational drag-dominant regimes. The wind-tunnel turbulence characteristics and the number of turbine revolutions over which the experimental C P values were averaged are not reported.

Computational settings
For the three wind turbines, a two-dimensional (2D) computational domain with dimensions of 35d Â 20d is made, where "d" denotes the turbine diameter. The distance from the turbine center to the domain inlet and outlet is 10d and 25d, respectively. The domain width, W, is 20d resulting in a 2D blockage ratio (d/W) of 5%. The domain is divided into a fixed domain and a rotating core to facilitate the turbine rotation. The diameter of the rotating core is 1.5d. The dimensions of the domain are based on the extensive sensitivity analyzes provided in Refs. [62,63]. Fig. 2a shows the domain for Turbine 1. Similar domains are used for the other two turbines.
Note that in the present study simulations are performed using 2D computational domains because: (i) the main objective of this paper is limited to analyze the ability of the different Reynoldsaveraged turbulence models to predict the flow complexities associated with VAWTs in the rotor midplane where the 3D tip effects are insignificant; (ii) two of the employed experiments [51,72] are focused on the turbine midplane to stay clear of and avoid the 3D tip effects; (iii) earlier extensive comparative studies of 2D and 2.5D simulations have shown that 2D simulations can provide a reliable representation of the flow in the turbine midplane [32,74]; and (iv) the high computational cost of 3D and 2.5D simulations makes such an extensive analysis infeasible.
The computational grids for Turbines 1, 2 and 3 consist of 153,303, 402,999 and 572,412 quadrilateral cells, respectively. The grid is refined near the walls where the maximum and averaged y þ values are 3.0 and 1.0, respectively. A grid convergence analysis has been performed and presented in detail in Refs. [36,63]. Fig. 2bed shows the different regions of the grid for Turbine 1. Similar grid topology and resolution are used for the other two turbines.
A uniform velocity is set at the domain inlet with a zero static gauge pressure at the outlet and symmetry condition on the domain sides. A no-slip condition is used on the turbine walls. A non-conformal sliding grid interface is employed between the rotating and fixed domains. The approach-flow, i.e. inlet, and the incident-flow total turbulence intensities for the three turbines are given in Table 1. The incident value represents the value at the turbine incidence calculated using an empty domain where such a value is typically lower than the inlet value due to the turbulence decay in the computational domain [63,75e77]. The turbulence length scale at the domain inlet is set equal to the turbine diameter for each case ( Table 1).
The four-equation g-Re q (transition SST) model is used, as the reference turbulence model, for the validation studies presented in Section 3.2. ANSYS Fluent 16.1 [78] is employed for the simulations. The 2D incompressible unsteady Reynolds-Averaged Navier-Stokes (URANS) equations are solved with second-order (spatial and temporal) discretization and the SIMPLE scheme for pressurevelocity coupling [79]. The details of the employed turbulence models are discussed in Section 4. Steady RANS solution is used to initialize the transient simulations. In the transient simulations, an azimuthal increment, dq, of 0.1 with 20 iterations per time step is employed. Each transient simulation is continued for 20 turbine revolutions to reach statistical convergence. The presented results are, thus, for the 21st turbine revolution. The employed azimuthal increment and the convergence criterion are based on extensive time-step sensitivity analysis and convergence study presented in Refs. [62,63].
In the present study, the turbine plane of rotation is split based on the azimuth q into the following four quartiles: (i) upwind, 45 q < 135 ; leeward, 135 q < 225 ; downwind, 225 q < 315 ; and windward, 315 q < 45 .

Validation
Three sets of validation studies have been performed for the three turbines presented in Table 1.
In the first validation study (see Fig. 3a), the normalized strength of the circulation of (clockwise) leading-edge separated/ shed vortex for a turbine blade in dynamic stall in transitional Re c of 50,000 at a low freestream TI of 0.015% is compared with the experimental data by Ferreira et al. [51]. Note that the CFD results correspond to a single turbine revolution. The quantified comparison is shown in Table 2. Given the flow complexity and the reported experimental uncertainty, an overall good agreement between the present CFD results and the experimental data is observed, where for most azimuthal angles, the absolute deviation between the CFD and the experiment is within the reported experimental uncertainty ( Table 2). The highest deviation is observed at q ¼ 90 , where CFD overestimates the experimental value. In general, the CFD results seem a bit shifted to the left (smaller q) compared to the experimental data. The observed deviations could be attributed to: RANS limitations to model the complex dynamic stall phenomenon: Previous studies have revealed that Reynoldsaveraged modeling predicts an earlier occurrence of the dynamic stall phenomena [59]. Phase-averaging of the experimental data: The experimental data is phase-averaged over dozens of turbine revolutions while the computational cost of such phase-averaging is prohibitively high, thus, the CFD results correspond to a single turbine revolution.
Method for the calculation of circulation strength: The experiment reports the size of the region to integrate the vorticity field in order to calculate the circulation strength as a main source of the reported uncertainty. In the present study, the integration regions are selected to cover the separated/shed vortices as will be described in Section 5.3. Possible differences in the size of the integration windows could be a source of the deviation. In the second validation study, the time-averaged (over the last turbine revolution) normalized streamwise and lateral velocities along the lateral direction, À1.5 y/R 1.5, in the turbine wake at different downstream locations, x/R ¼ 1.5e4.0, are compared with the PIV measurements by Tescione et al. [72]. The results are presented in Table 3. Average deviations of 6.8%e16.0% and 2.2%e2.8% are found for the streamwise and lateral velocity components. Fig. 3b shows the comparison for the streamwise velocity at x/ R ¼ 2.0. The observed deviations could be due to: (i) the differences in data averaging, where the experimental data are phase-averaged over 150 turbine revolutions while the CFD results are timeaveraged over one turbine revolution; (ii) geometrical simplification, where the less-aerodynamic turbine components, e.g. spokes and connections, are ignored in the CFD modeling; and (iii) the limitations of the two-dimensional URANS modeling. This validation study is presented in detail with a more extensive discussion of possible reasons for the observed deviations in Refs. [35,63].
In the third validation study (see Fig. 3c), the turbine power coefficient is compared with the experimental data by Castelli et al. [73]. The comparison of the CFD and the experimental results are shown in Table 4, where the maximum absolute deviation in C P , i.e. 23.2%, corresponds to the highest l of 3.08. The observed deviations could be due to: (i) the experimental uncertainties, where several important operational and geometrical parameters, e.g. the turbulence characteristics of the tunnel and the location of the bladespoke connection, are not reported; (ii) geometrical simplifications, where the less-aerodynamic turbine components, e.g. spokes and connections, are ignored in the CFD modeling; and (iii) the deficiencies of 2D URANS. This validation study with a more extensive discussion of possible reasons for the observed deviations are presented in Refs. [32,36,62].
The three validation studies, using dissimilar experimentallyreported parameters for turbines with different geometrical and operational characteristics, further increase the confidence on the accuracy of the present numerical results.

Turbulence models
Seven one-to four-equation eddy-viscosity models are employed where three of the models also include additional transition modeling to predict the laminar-to-turbulent transition ( Table 5). In addition, the inviscid modeling, as a fast and lowcomputational cost method, is also evaluated. Note that the Reynolds stress models (RSM) are not considered in the present study due to their inferior computational economy, robustness and convergence deficiencies, compared to the eddy-viscosity models for unsteady flow simulations [80]. The computational cost of RSM, for example, could be an order of magnitude higher than the twoequation models for airfoil studies [81].   Table 3 Comparison between CFD and experiment [72] for Turbine 2.
x/R A short description of the employed turbulence models is briefly presented: 1) The one-equation Spalart-Allmaras model (SA) [82]: This model has been originally developed for aerospace applications of wall-bounded flows. The model has been successful in modeling flows over airfoils with adverse pressure gradient. Later, the model has also become popular for turbomachinery applications. It solves a modeled transport equation for the kinematic turbulent viscosity, m t . Turbulent kinetic energy, k, is indirectly calculated based on a combination of both vorticity and the strain rate tensor. In the present study, the alternative formulation by Ref. [83] is used that decreases the turbulent viscosity in the region where the vorticity exceeds the strain rate, e.g. vortical flows, and, thus, improves the performance of the model in predicting the impact of the rotation in damping the turbulence inside a vortex.
2) The two-equation RNG k-ε model (RNG) [84]: This model, which is in similar form to the standard k-ε model, is derived using a statistical method, i.e. the renormalization group (RNG) theory, in order to improve the flow predictions for rapidlystrained and swirling flows. The RNG model calculates the turbulent Prandtl numbers analytically, rather than using constant values as in the standard model. The model also accounts for the low-Reynolds number effects in contrast to the standard model, which is basically a high-Reynolds number model. In the present study, the enhanced wall treatment [85,86] including the influence of pressure gradients is employed for the near-wall modeling.

3) The two-equation Realizable k-ε model (RKE) [87]:
This model is developed to improve the flow predictions of the standard k-ε model for flows with strong streamline curvature, rotational and vortical features. It includes a modified transport equation for the dissipation rate, ε, and is derived based on the exact transport equation of the mean-square vorticity fluctuation. The modified equation makes the model 'realizable', i.e. consistent with the turbulent flow physics, which is not the case for the standard and RNG k-ε models. In addition, an alternative formulation for the turbulent viscosity is used. In the present study, the enhanced wall treatment [85,86] including the influence of pressure gradients is employed for the near-wall modeling.
4) The two-equation SST k-u model (SST) [88]: This model combines, through a blending function, the accuracy and robustness of the standard k-u model [89] in the near-wall region with the freestream insensitivity of the standard k-ε model [90], in a transformed formulation, away from the wall. In addition, the transport of the turbulent shear stress is included in the formulation of the turbulent viscosity to improve the predictions of the flow separation onset for flows on smooth surfaces with adverse pressure gradients, such as flows over airfoils. In addition, the turbulent production limiter by Menter [91] is included to avoid excessive turbulent kinetic energy prediction in stagnation regions. The limiter is known not to influence the shear layers.

5)
The three-equation SST k-u with additional intermittency transition model (SSTI) [78]: An additional transport equation for the intermittency, g, is coupled with the SST k-u transport equations to enable the SST model to predict the laminar-to-turbulent transition. The SSTI model is a correlation-based model, implying that empirical correlations for the critical momentum-thickness Reynolds number, Re qc , sensitized by the freestream turbulence intensity, Re, pressure gradient, etc., are employed to trigger the model. The model is a successor of the four-equation transition SST k-u model (TSST) [92], which eliminates the need for the second transport equation for the momentum-thickness Reynolds number, Re q , employed in the TSST k-u model. Therefore, the model is less computationally expensive, less complex and includes less number of correlations compared to the TSST k-u model. While the coefficients have been calibrated for the classical boundary layer flows, fine-tuning for specific applications can also be conveniently performed. Similar to the other engineering transition models, the model is only applicable for the wall-bounded flow, thus, will predict the free shear flows as fully turbulent. The turbulent production limiters by Menter [91] and Kato and Launder [93] are included to avoid excessive turbulent kinetic energy prediction in stagnation regions.
6) The three-equation k-k l -u transition model (KKL) [94]: The k-k l -u transition model is developed based on the standard k-u model [89], and solves an additional transport equation for the laminar kinetic energy, k l , together with the k-u transport equations in order to predict the transition onset and length. The additional equation for k l is to predict the low-frequency large-scale velocity fluctuations within the boundary layer, which can signal the pre-transition processes. The model accounts for the natural and bypass transitions. As opposed to the TSST and SSTI models, this model is more physics-based rather than correlation-based.

7)
The four-equation transition SST model (TSST) [92,95]: The transition SST k-u model, also known as g-Re q model, couples two additional transport equations, one for g and one for Re q , with the SST k-u transport equations to enable the SST k-u model to predict the transition onset and length. The intermittency equation is used to trigger the transition locally while the momentum- and is capable of predicting the natural, separation-induced and bypass transitions for wall-bounded flows. Similar to the other engineering transition models, the model will predict the free shear flows as fully turbulent. The turbulent production limiters by Menter [91] and Kato and Launder [93] are included to avoid excessive turbulent kinetic energy prediction in stagnation regions.

Turbine wake
Figs. 4 and 5 show the contour plots of time-averaged (over the last turbine revolution) normalized (by freestream velocity) streamwise and lateral velocity components for Turbine 2 calculated using the different turbulence models compared against the experimental data by Tescione et al. [72]. The experimental contour is phase-averaged over 150 turbine revolutions. The turbine is operating within the optimal range at a moderate tip speed ratio of 4.5. The comparison implies that: The overall wake expansion observed in the experiment is underpredicted by the different turbulence models. This can be inferred from the wake width at x/R ¼ 3.0 using the vertical dashed lines shown in Figs. 4e5. The characteristics of the turbine near wake, e.g. the traces of the blades' wake and the shed vortices observable at the most windward position, are not well predicted by the inviscid, the SA and the KKL models (Fig. 5). The wake of the turbine shaft, highlighted in Fig. 4 as a comparatively lower velocity region, at y/R z 0 and 0 x/ R 1.0, is experimentally found to be deflected towards the positive y-direction due to the counter-clockwise rotation of the shaft. This is not well predicted by the inviscid and the KKL models (Fig. 4).
To provide a more detailed analysis, Figs. 6 and 7 present profiles of the time-averaged (over the last turbine revolution) normalized streamwise and lateral velocity components along the lateral direction, À1.5 y/R 1.5, at different downstream locations, x/R ¼ 1.5, 2.0, 2.5, 3.0, 3.5 and 4.0, in the wake of Turbine 2 compared with the experiment by Tescione et al. [72]. The experimental data is phase-averaged over 150 turbine revolutions. Tables 6 and 7 list the average deviations of the time-averaged (over the last turbine revolution) normalized streamwise and lateral velocity components at different downstream locations in the wake.
For given y/R values, considerable variations can be seen in the predicted normalized streamwise and lateral velocity by the different models. Note that the predicted overall wake profiles are in reasonable agreement with the experiment and the average deviations, between the predictions and the experiment, are comparable for different models. In addition, all the models clearly fail to predict the experimentally-observed wake asymmetry, which could be due to the geometrical simplification in the CFD simulations and/or the limitations of the RANS modeling. Nevertheless, none of the models stands clearly superior to the others. The similar predictions of the different models in the wake of the VAWT could be associated with the transition to bluff-body aerodynamics, which occurs in the wake of the turbine [100], meaning that the turbine wake, neglecting the asymmetry, somehow resembles the generic case of flow over circular cylinders, where the different models are known to show reasonable predictions of the mean flow [101]. Fig. 8 illustrates the contour plots of the normalized velocity magnitude calculated using the different turbulence models. The length of the turbine wake is approximately the same for the different turbulence models, while the inviscid modeling predicts a considerably shorter wake, where the wake length is defined as the length along the centerline in which V =U ∞ ¼ 0.7. The predicted wake length values by the inviscid, SA, RNG, RKE, SST, SSTI, KKL and TSST models, shown using dotted lines in Fig. 8, are approximately 7d, 9d, 8.5d, 9d, 9.5d, 9d, 8.5d and 9d, respectively, where d denotes the turbine diameter.

Turbine power coefficient
The turbine power and thrust coefficients, C P and C T , are the key parameters for the analysis of the aerodynamic performance of wind turbines. Therefore, its accurate prediction is of paramount importance. Fig. 9 shows the power and thrust coefficient values, C P and C T , for the three turbines calculated using the different turbulence models. For Turbine 3, the predicted C P values are also compared against the experimental data by Castelli et al. [73] within a wide range of l, i.e. 1.44e3.29. Note that this range of l covers the three regimes where the turbine operates in dynamic stall (l < 2.33), optimally (2.33 l 2.68), and in a highlyrotational flow (l > 2.68). Table 8 presents the deviations between the CFD results and the experimental data [73]. The following observations can be made for Turbine 3 (Fig. 9a and e): There are significant variations for both C P and C T . This is especially the case for C P . The inviscid modeling results in a very large overprediction of the C P values, compared to the experiment, for the whole range of l. This overprediction is not systematic and reduces within the turbine optimal regime. The overprediction is because such a modeling inherently neglects the viscosity, therefore, cannot predict the essential flow features driving the turbine aerodynamic performance, i.e. the boundary layer development and separation on the blades. The C P predictions using the ε-based models, RNG and RKE, have significant deviations from the experimental data where neither of the two models stand superior. The inaccuracy is because the ε-based models suffer from lack of sensitivity to adverse pressure gradient and overprediction of the turbulent kinetic energy   to their more accurate near-wall formulation [91,102].
To further analyze the performance of the turbulence models, the C P values for Turbine 1 and Turbine 2 are also presented in Fig. 9b and c. Similar to Turbine 3, the predictions using the inviscid and the ε-based models are quite different than those using ubased models, where the differences are more pronounced for higher l. Although the predications of the four u-based models themselves are almost clustered together, there are some variations between them at some values of l. This necessitates a more detailed analysis of the instantaneous moment to elucidate where the differences originate from. The variations in the predicted C T values by the different models, shown in Fig. 9eeg, are less pronounced than for the C P values. ii. The downwind quartile, 225 q < 315 , corresponding to the region where the blades interact with the shed vortices from the blades further upstream. The difference predictions of the models in this region are also strongly connected to the predictions of the blade aerodynamics further upstream in the upwind, which will also be discussed in Section 5.3.
Focusing the analysis on the u-based models, which already exhibited the least deviation from the experimental C P values, reveals that the instantaneous moment calculated using these models, namely SST, SSTI, TSST and KKL, are also noticeably different and, therefore, the respective C P values for at least some of them might have been influenced by error cancellation. This necessitates further analysis of the blade aerodynamics for a complete understating of the performance of these turbulence models.

Blade aerodynamics
Detailed comparison of the C P values predicted using the different turbulence models, for a wide range of l, against the experimental data in Fig. 9a already demonstrated the inaccuracy of the inviscid, SA, RNG and RKE models for CFD simulations of VAWTs. Only the u-based models showed a reasonable agreement with the experiment. Further investigation of the turbine instantaneous moment revealed considerable variations in the predicted values even between the four u-based models which necessitated further analysis. Therefore, in this section, the analysis steers towards comparing the VAWT blade aerodynamics against the experimental data by Ferreira et al. [51]. For the remainder of the paper, the discussion only focuses on the u-based models (i.e. SST, SSI, 'TSST', 'KKL') while for completeness the results of the other turbulence models listed in Table 5 are also presented. Fig. 11 shows the evolution of the instantaneous (clockwise) leading-edge separated/shed vortex compared against the PIV results by Ferreira et al. [51] for Turbine 1, operating at l ¼ 2.0 in dynamic stall. Note that the plots exclude the trailing-edge separation, which is separately presented later in this section. Fig. 12 illustrates the contour plots of the dimensionless zvorticity at selected critical azimuthal positions to clarify the growth of the unsteady separation predicted by each model. Fig. 13 shows the normalized strength of the circulation of (clockwise) leading-edge separated/shed vortex compared to the experimental data in Ref. [51] where the deviations are reported in Table 9. The values exclude the trailing-edge separation.
Note that the ε-based models completely fail to predict the formation of the dynamic stall vortex (DSV) near the leading edge while instead predict a trailing-edge stall (see Fig. 12c), which is  inconsistent with the experimental observations shown in Fig. 11a. For the same reason, Fig. 13 also presents no value corresponding to the ε-based models. Based on Figs. 11e13, the following observations can be made: i. The flow prediction by the 'KKL' model is not in line with the experimental observations. Although this model predicts the formation of a leading edge vortex, it clearly fails to predict its shedding. This is apparent from Figs. 11h and 12d where a leading-edge vortex is formed on the blade (e.g. at q ¼ 55 in Fig. 12d), which grows in size by increasing q and eventually the flow reattaches without shedding the leading-edge vortex. The inconsistency of the KKL predictions compared to the experiment is also clearly shown in Fig. 13 where the strength of the circulation of leading-edge separated vortex suddenly drops for q > 108 due to the reduction in size and eventual reattachment of the leading-edge vortex. ii. The best agreements with the experiment are provided by the SST model and its variants (SSTI and TSST) where the essential flow features, namely the emergence of the dynamic stall vortex, its growth, shedding and propagation downstream in the rotor plane, are well predicted (see Fig. 11f, g, i and Fig. 12eeg).
In addition, the predicted strength of the circulation of leadingedge separated/shed vortex by the models is in line with the experimental data and for the majority of the studied cases, the deviation with the experimental data is within the reported experimental uncertainty (see Fig. 13). Among the SST variants, the SST model, rather than its transitional SSTI and TSST models, might initially seem to have the best agreement with the experiment in terms of the spatial location of the shed leading-edge and trailing-edge vortices (see Figs. 11f and 14f), while such vortices predicted by the SSTI and TSST models are positioned at comparatively lower y/c values than the experiment (see Fig. 11g, i and 14g, i). However, as the SST predictions in Figs. 11f and 12e show, right after the shedding of the dynamic stall vortex at q z 80 a secondary vortex (SV) is also formed and shed, which is later merged with the dynamic stall vortex forming into a larger vortex (see Fig. 12e at q ¼ 100 ). Therefore, the shed vortices shown in Fig. 11f for q ¼ 108 , 133 , 158 and 233 correspond to a merged 'dynamic stall vortex e secondary vortex' (DSV-SV), rather than only the dynamic stall vortex (DSV), and this is less consistent with the experimental observations. Note that the dynamic stall vortex refers to the well-defined primary vortex shed during the dynamic stall process and the vortices shed afterwards are referred to as secondary/tertiary vortices.
As already mentioned, the shed leading-edge and trailing-edge vortices predicted by the SSTI and TSST models are positioned at comparatively lower y/c values than the experiment. This is explained as the two models predict an earlier separation and shedding of the dynamic stall vortex than the experiment. The earlier shedding of the dynamic stall vortex results in the vortex detaching from the blade at a lower y/c, which therefore, also ends up in a lower y/c position. Considering the transitional regime of  the flow at Re c of 50,000 and TI of 0.015%, the early prediction of the separation and the shedding of the dynamic stall vortex is thought to be due to one or multiple of the following reasons: (i) the limitations of the RANS modeling for such complex flow phenomena, as also reported in the literature [59]; (ii) the relatively high blockage in the wind tunnel, i.e. blockage ratio of z18%, which might have resulted in higher flow velocity and thus Re c than the reported values; (iii) higher incident-flow turbulence intensity than the value reported further upstream in the wind tunnel; (iv) notperfectly-smooth blade surface in the experiment. Fig. 15aeh illustrate the spatiotemporal contours of surface pressure and skin friction coefficient coefficients along the blade suction side over the turbine half-revolution for Turbine 1 at l ¼ 2.0. Fig. 15i-l and 15m-p illustrate the spatiotemporal contours of skin friction coefficient along the blade suction side over the turbine half-revolution for Turbine 2 at l ¼ 4.5 and Turbine 3 at l ¼ 1.68, respectively. The figures are helpful to further highlight the differences in the prediction of the models regarding the details of the laminar separation bubble (LSB), the dynamic stall vortex (DSV), the secondary vortex (SV) and the trailing-edge roll-up vortex (TEV) introduced using the indications on the figure. Note that during the dynamic stall process, as the angle of attack increases, a laminar separation bubble would typically form near the airfoil leading edge. Further increase of the angle of attack could lead to growth and a sudden bursting of the bubble. By the bursting of the laminar separation bubble, suction along the blade leading edge would abruptly collapse and a well-defined vortical structure, termed as dynamic stall, would form. The dynamic stall vortex would grow in size while traveling downstream along the blade and would finally shed from the blade. By detachment of the dynamic stall vortex from the blade surface, secondary and tertiary vortices might also form.
The following conclusions can be made from Fig. 15: For Turbine 1, the KKL model predicts the formation and the growth of the laminar separation bubble. However, it cannot predict the shedding of the dynamic stall vortex. Note that this can lead to a significant inconsistency with the experimental observations as shown in Figs. 11h and 12d. For Turbine 1, the SST predictions are different from those of its transitional SSTI and TSST models mainly with respect to characteristics of the dynamic stall vortex, the secondary/tertiary vortices, the trailing-edge vortex and the incipient trailing-edge separation. There exists noticeable similarity between the predictions by the SSTI and TSST models. Similar observations are also made for Turbine 3. For Turbine 2, considering the transitional flow regime with Re c of z150,000 and TI of 5%, a laminar separation bubble is known to be present on the blade [106,107]. Note that the SST model is unable to predict the laminar separation bubble (see Fig. 15j).
Neglecting the laminar separation bubble, also impacts the extent of the incipient trailing-edge separation, which is found to be larger for the SST model than that for the SSTI and TSST models.
Therefore, it can be concluded that the transitional SSI and TSST models are to be preferred over the SST model for CFD simulation of VAWTs due to their ability to predict the laminar-to-turbulent transition, which occurs on the blade. Note that although, there exists differences between the predictions of the SSI and TSST models, e.g. on the laminar separation bubble characteristics (see Fig. 15k and l), so far they are the only eddy-viscosity turbulence model to exhibit reasonable agreement with the experimental data when comparing both detailed blade aerodynamics as well as the turbine power performance (Section 4).

Discussion
Although the present analysis yielded the SST-based models as the most accurate eddy-viscosity turbulence models for the CFD simulations of VAWTs, it also revealed several shortcomings of such models in the prediction of the relevant flow features, e.g. the dynamic stall. Further investigation via detailed comparison against high-spatiotemporal-resolution experimental data, which currently e to the best of our knowledge e do not exist, and the more advanced scale-resolving simulations is needed in order to provide recommendations to improve the SST models' predictions.
While the focus of the study is on the accuracy of the model predictions, two other important issues regarding the turbulence models, i.e. computational cost and robustness, might also be of the interest for the reader. Regarding the robustness, problems with convergence were experienced for the RNG and KKL models, where both were hard to converge. Regarding the computational cost, focusing on the SST-based models as the models of interest for this application and considering the SST k-u model as a reference, the additional equation for the SSTI model and the two additional equations for the TSST model increased the computational cost by z 14% and z 30%, respectively.
In the present study, the employment of the curvature  correction for the turbulence models was not considered, where this can be interesting to be investigated. In addition, a detailed comparison of the 2D prediction of the SST models with the 2.5D predictions can also be recommended to clarify the threedimensional modeling effects.

Conclusions
A critical analysis of seven eddy-viscosity turbulence models, namely the Spalart-Allmaras (SA), RNG k-ε, realizable k-ε (RKE), SST k-u, SST k-u with an additional intermittency transition model (SSTI), k-k l -u transition model (KKL) and g-Re q transition SST model (TSST), for URANS simulations of vertical axis wind turbines (VAWTs) is presented. A comparison with the inviscid modeling is also included in the comparison.
The analysis is based on evaluation against three sets of experiments using turbines with dissimilar geometrical and operational characteristics. In addition, the analysis is focused on three different parameters, namely the turbine wake, the turbine power performance, and the blade aerodynamics. Such a wide diversity further supports the generality of the findings.
The best-performing turbulence model is selected as the model which provides the best agreement with all the three sets of measurements. Note that this should be done only after grid convergence analysis has shown that numerical (discretization) errors remain sufficiently low. Therefore, grid convergence analysis was included in the present paper.
The results show that: The SST model variants (SST k-u, SSTI and TSST) exhibit reasonable agreement with the three experiments in predicting the details of the VAWT's vortex evolution in the dynamic stall, the wake mean streamwise and lateral velocity profiles in the optimal operational regime and the turbine power coefficient within a wide range of tip speed ratios. Therefore for URANS simulations of VAWTs, the SSTbased models are found to be the most accurate among the models evaluated. Note that within the transitional flow regime (approximately Re c < 250,000), the transitional SST versions (SSTI and TSST) will be the models of choice due to their ability to account for the laminar-to-turbulent transition on the blade surface. The recommendations can also extend to the employment of the SST-based models in the RANS region of the hybrid RANS/LES models. The inviscid, SA, RNG, RKE and KKL models fail to provide accurate predictions of the aerodynamic performance of VAWTs. This is mainly attributed to the following reasons: (i) The inviscid modeling neglects the viscous boundary layer.
Given the essential influence of the boundary layer development and the unsteady separation for the prediction of the for Turbine 1 at l ¼ 2.0.
aerodynamic performance of VAWTs, the inviscid modeling leads to a large non-systematic overprediction of turbine C P ; (ii) The SA model is unable to account for the highly-rotational features of the flow around the VAWTs, which results in the unrealistic prediction of the turbine C P ; (iii) The RNG and RKE models cannot predict the occurrence of the experimentally-observed dynamic stall on the turbine blade while instead they predict a gradual trailing-edge stall, which results in substantial differences in the blade aerodynamics, flow topology and turbine power coefficient; (iv) The KKL transition model cannot predict the shedding of the dynamic stall vortex, which therefore significantly changes the predicted blade aerodynamics and the turbine power performance.

Acknowledgement
The authors would like to acknowledge support from the European Commission's Framework Program Horizon 2020, through the Marie Curie Innovative Training Network (ITN)