Towards accurate CFD simulations of vertical axis wind turbines at different tip speed ratios and solidities_ Guidelines for azimuthal increment, domain size and convergence

The accuracy of CFD simulations of vertical axis wind turbines (VAWTs) is known to be significantly associated with the computational parameters, such as azimuthal increment, domain size and number of turbine revolutions before reaching a statistically steady state condition (convergence). A detailed review of the literature, however, indicates that there is a lack of extensive parametric studies investigating the impact of the computational parameters. The current study, therefore, intends to systematically investigate the impact of these parameters, on the simulation results to guide the execution of accurate CFD simulations of VAWTs at different tip speed ratios (λ) and solidities (σ). The evaluation is based on 110 CFD simulations validated with windtunnel measurements for two VAWTs. Instantaneous moment coefficient, Cm, and power coefficient, CP, are studied for each case using unsteady Reynolds-averaged Navier-Stokes (URANS) simulations with the 4-equation transition SST turbulence model. The results show that the azimuthal increment dθ is largely dependent on tip speed ratio. For moderate to high λ, the minimum requirement for dθ is 0.5° while this decreases to 0.1° at low to moderate λ. The need for finer time steps is associated to the flow complexities related to dynamic stall on turbine blades and blade-wake interactions at low λ. In addition, the minimum distance from the turbine center to the domain inlet and outlet is 15 and 10 times the turbine diameter, respectively. It is also shown that 20–30 turbine revolutions are required to ensure statistically converged solutions. The current findings can serve as guidelines towards accurate and reliable CFD simulations of VAWTs at different tip speed ratios and solidities.

It is widely recognized that the accuracy and reliability of CFD simulations of VAWTs can be very sensitive to the computational settings. For example, earlier studies have shown the significant importance of azimuthal increment [6,[47][48][49][50][51] and domain size [6,50,52,53]. However, the results of these studies are case-specific and limited in scope where the conclusions did not (aim to) provide generally valid recommendations and guidelines to support future CFD studies of VAWTs. Therefore, there is no consensus on the employed computational settings for the Unsteady Reynolds-averaged Navier-Stokes (URANS) simulations of VAWTs in the literature. This is especially the case for azimuthal increment, domain size and convergence criterion (defined as number of turbine revolutions before reaching a statistically steady state condition), as an extensive review of CFD studies of VAWTs indicates that: • The range of the azimuthal increment (which is defined as the number of degrees that a turbine rotates per time step) varies from 0.03° [6,48] to 10° [6,54,55].
• There is no clear consensus on the convergence criterion. Some studies [56] considered convergence from the 4th turbine revolution while other studies considered this from 100th turbine revolution [6].
To the best of our knowledge, the first study to derive the minimum requirements for computational settings of CFD simulations was performed by Rezaeiha et al. [6]. The results of that study showed that: • Employment of an azimuthal increment of 5°would lead to 17% underprediction of turbine C P .
• Using a computational domain with 5D distance from the turbine center to domain inlet results in more than 5% overprediction of turbine C P.
• Sampling data starting after the 5th turbine revolution would result in more than 20% overprediction of turbine C P .
It should be noted that the study by Rezaeiha et al. [6] was only performed for a low solidity (σ = 0.12) and moderate tip speed ratio (λ = 4.5). Those parameters were selected to minimize the flow complexity associated to dynamic stall [8], blade-wake interactions [10] and flow curvature effects [12]. Nevertheless, the more complex flows for lower tip speed ratios might bring more challenges to CFD simulations, which might influence the minimum requirements for computational settings and parameters. To the best of our knowledge, there is a lack of sensitivity analysis for minimum requirements for computational settings and parameters for CFD simulations of VAWTs at different tip speed ratios and solidities. The current study, therefore, intends to systematically investigate the impact of azimuthal increment, domain size and convergence criterion on the simulation results at different tip speed ratios and solidities to guide the execution of accurate CFD simulations of VAWTs.
The outline of the paper is as follows: the geometrical and operational characteristics of the wind turbine and the computational settings are described in Section 2.1-2.3. A reference case for the sensitivity analysis is defined in Section 2.4. Solution verification and validation studies are presented in Section 3. The sensitivity of the minimum requirements for azimuthal increment, domain size (distance from the turbine center to the domain inlet and outlet) and convergence criterion to tip speed ratio and solidity is extensively discussed in Sections 4-7, respectively. The discussion and the conclusions are provided in Section 8-9.

Geometrical and operational characteristics of the reference case
An H-type vertical axis wind turbine with 2 blades of symmetric NACA0018 airfoil section is used. The turbine has a diameter (D) of 1 m and a shaft diameter of 0.04 m. The turbine and the shaft rotate counter-clockwise with the same rotational speed. The freestream velocity (U ∞ ) is 9.3 m/s with a total turbulence intensity of 5%. The same freestream velocity as in the validation study for the reference turbine (see Section 3) is used. The turbine chord-based Reynolds number (Re c ) is 176,000. The turbine for the reference case has a solidity (σ) of 0.12, a chord length (c) of 0.06 m and a turbine rotational speed (Ω) of 83.8 rad/s (800 rpm), leading to a tip speed ratio (λ) of 4.5. Tip speed ratio and solidity are defined using Eqs. (1) and (2). The main geometrical and operational characteristics of the turbine for the reference case are shown in Fig. 1 and described in Table 1. Note that these characteristics are the same as those used in the earlier study by Rezaeiha et al. [6] in which guidelines for the minimum domain size and azimuthal increment for a low-solidity (σ = 0.12) VAWT operating at a moderate tip speed ratio of 4.5 were provided.   [26]. In this study, the blade-aspect-ratio for the reference turbine is 16.67.
Earlier studies have shown that for CFD simulations of VAWTs, a blockage ratio of less than 5% (a minimum domain width of 20D) is required to minimize the effect of side boundaries on the results [6]. For 2D simulations, the blockage ratio is defined as D/W where D and W correspond to the turbine diameter and domain width. As the blockage ratio is defined based on the full diameter of the turbine, it is assumed to be not significantly affected by the turbine tip speed ratio and solidity. Therefore, in the present study, the domain width of 20D is employed. In addition, the impact of the diameter of the rotating core d c (the rotating region of the domain which is connected to the fixed domain using a sliding grid interface as shown in Fig. 2) was found to be very negligible (< 0.04%) on the calculated turbine performance [6]. In the present study, therefore, the diameter of the rotating core is 1.5D. The rotating core is joined to the surrounding fixed domain using a sliding grid interface. The distance from the turbine center to domain inlet (d i ) and outlet (d o ) for the reference case are 5D, 25D, respectively. Four quartiles are defined on the plane of rotation of the blades [58], namely windward 315°≤ θ < 45°, upwind 45°≤ θ < 135°, leeward 135°≤ θ < 225°and downwind 225°≤ θ < 315°. Azimuthal angle θ is defined based on the position of the top blade (which is currently positioned the most windward) as shown in Fig. 2.
The computational grid is made up of quadrilateral cells everywhere. The maximum cell skewness in the grid is 0.55 and the minimum orthogonal quality is 0.21. The maximum y + values on the surfaces of the blades and the shaft are 4 and 2, while the average y + values are 1.37 and 1.02, respectively. The total number of cells is approximately 400,000. Fig. 3 shows the computational grid. The boundary conditions are uniform velocity inlet, zero gauge pressure outlet, symmetry sides and no-slip walls. The turbulent length scale at the domain inlet and outlet is selected based on the length scale of the interest, i.e. turbine diameter, 1 m.

Other numerical settings
Incompressible unsteady Reynolds-averaged Navier-Stokes (URANS) simulations are performed using the commercial CFD software package ANSYS Fluent 16.1. The computational settings are based on an earlier solution verification and validation study [6]. Secondorder schemes are used for both temporal and spatial discretization. The SIMPLE scheme is used for pressure-velocity coupling.
Turbulence is modeled using the 4-equation transition SST turbulence model, also known as γ-θ SST turbulence model [59], which is an improvement over the two-equation k-ω SST turbulence model [60,61]. The model solves two additional equations for intermittency γ and momentum-thickness Reynolds number Re θ to provide a more accurate prediction of laminar-to-turbulent transition for wall-bounded flows where such prediction is critical to correctly model the boundary layer development and to calculate the loads on the walls [62,63]. Further details of the γ-θ SST turbulence model are provided by Menter et al.    [59] and Langtry et al. [63]. Modeling of the laminar-to-turbulent transition using intermittency-based models has also been proven successful by Suzen et al. [64,65], Walters and Leylek [66], Cutrone et al. [67] and Genç et al. [68][69][70][71][72][73] for flow on airfoils, which further provides confidence for application of these models. Note that the ω-based turbulence models benefit from a y + -insensitive formulation which means that the model, based on the y + of the cells near the wall, will either opt for the low-Reynolds modeling, to resolve the buffer layer and viscous sublayer near the wall, or for wall functions [61,74,75]. The unsteady simulations are initialized with the solutions of steady RANS simulations and continued for 20 turbine revolutions. The azimuthal increment (dθ) employed for the unsteady calculations for the reference case is 0.1°. Each time step consists of 20 iterations to have the scaled residuals < 1 × 10 −5 . The instantaneous moment coefficient is sampled at the 21st turbine revolution and the values are employed to calculate the turbine power coefficient.

Reference case
The computational settings for the reference case is described in Table 2. The sensitivity analysis is performed for domain size, azimuthal increment and convergence criterion (defined as number of turbine revolutions before reaching a statistically steady state condition) by applying systematic changes to the reference case. In Sections 4-7, one of the computational settings is varied, while all others are kept the same as in the reference case. In the current study, tip speed ratio is modified by altering the turbine rotational speed and solidity is modified by altering the blade chord length.
Because the flow velocity coming to the blade, the experienced velocity, is the vector sum of the freestream velocity and the rotational velocity, therefore, modifying the experienced velocity via changing either of the two is expected to have similar results on the derived minimum requirements. In the present study, the experienced velocity is modified by changing the rotational velocity. Note that changing the experienced velocity will change Re c and can also alter λ. While the impact of λ on the derived minimum requirements is investigated in the present study, the investigation of the dependence of the derived minimum requirements on Re c is proposed for future research.

Solution verification and validation
Solution verification and two sets of validation studies are performed. Because these studies have been published as separate papers [6,39], only the headlines are briefly repeated here. The verification studies were performed for the reference case, described in Section 2.4. The following observations were made [6]: • A grid sensitivity analysis was performed using three grids which were uniformly refined by √2. Grid Convergence Index (GCI) [76] was calculated based on C P values using a safety factor F s of 1.25. For the coarse-medium grid pair, GCI coarse and GCI fine was 6.1 × 10 −3 and 3.5 × 10 −3 , which represent 1.48% and 0.85% of the respective C P values.
• A time step sensitivity analysis was performed where an azimuthal increment of 0.5°was found to be sufficient. Further refinement of the time step to 0.1°had negligible (< 1%) effects on C P values.
• An iterative convergence study was performed where a statistically steady state solution was obtained after 20-30 turbine revolutions.
The difference between C P values at 20 and 30 revolutions with that of 100 revolutions was found to be 2.41 and 1.06%, respectively.
• A comparative study was performed where the results of the 2D domain were compared with that of a 2.5D with span of 6 cm (1 blade chord length). The comparison showed a systematic difference of approximately 6% in the C P values [38].
The first validation study [6], which was performed for the reference case, compared the CFD results against the experimental data by Tescione et al. [26]. The comparison of the time-averaged normalized streamwise and lateral velocities in the near wake of the turbine at different downstream locations from x/R = 2.0-3.0 showed a deviation of 8.6-11.8% for the streamwise component and 2.3-2.5% for the lateral component, respectively.
The second validation study [39], which was performed for a 3bladed VAWT, compared the CFD results against the experimental measurements by Castelli et al. [77]. The comparison of turbine C P for 8 different tip speed ratios 1.44 ≤ λ ≤ 3.29, shown in Fig. 4a, confirmed that the CFD results were in line with the experimental data. In addition, the figure shows that substantial improvement was achieved compared to the numerical results by Castelli et al. [77]. The good agreement observed for different tip speed ratios shows that the simulations could provide a reasonable prediction of the turbine performance also for lower λ (≤2.33) where the occurrence of dynamic stall results in a complex flow over the turbine blades. The occurrence of Table 2 The computational settings for the reference case. dynamic stall for λ ≤ 2.33 is confirmed from the instantaneous moment coefficient on the blades as shown in Fig. 4b. The possible explanations for the observed deviation between the present CFD results and the experimental data are extensively discussed in Refs. [6,39].

Azimuthal increment
The minimum required time step of the URANS simulations is driven by the time scales of the unsteady complex phenomena occurring in the flow. For the case of wind turbines, it is customary to correlate the time step with the turbine revolution and express it as azimuthal increment dθ. Azimuthal increment is then defined as the number of degrees that the turbine rotates per time step. As the complexity of the flow on VAWTs varies with the operating conditions and geometrical parameters, i.e. λ and σ, therefore, a sensitivity analysis is performed for dθ by applying systematic changes to the reference case. In this section, the sensitivity of the minimum requirements for the dθ to different tip speed ratios and solidities are investigated where the details of the computational parameters for all the cases are presented in Table 3. Fig. 5 indicates the instantaneous moment coefficient C m for the last turbine revolution versus azimuth θ for tip speed ratios 1.5 ≤ λ ≤ 5.5 calculated using different azimuthal increments 0.05°≤ dθ ≤ 1.0°. The turbine solidity is 0.12, which is the same as the reference case. It can be seen that for relatively high tip speed ratios of 4.5 and 5.5 the difference between the results for dθ ≤ 0.5°is negligible (Fig. 5a-b). The main reason is that moderate to high tip speed ratios of 4.5 and 5.5 correspond to a regime where the flow is mostly attached and the variations of the angle of attack on the blades during the turbine revolution are limited below the static stall angle. Consequently, the flow physics are less complex compared to those at lower λ values due to the absence of dynamic stall and blades' wake-blade interactions. A closer look at Fig. 5a and b indicates small deviations for dθ = 1.0°on the downwind side around θ = 270°, corresponding to the wake of the shaft [40]. This confirms that even at relatively high λ values, the interactions between the vortices shed from the shaft and the blades passing in the shaft wake (shaft's wake-blade interaction) can make predictions of the flow more challenging. In this case, finer azimuthal increments are required. For low to moderate tip speed ratios λ < 4.5 ( Fig. 5c-f), where the variations of the angle of attack on the blades get larger and consequently dynamic stall and blade-wake interactions would occur, the results are highly dependent on azimuthal increments. The difference is mainly observed in two regions: upwind (45°< θ < 135°) and downwind (225°< θ < 315°). The former is associated with the occurrence of the dynamic stall, which happens at such azimuthal angles and results in a sudden drop in C m . The latter corresponds to the interactions of the vortices shed from the blades in the upwind and the shaft with the blades passing downwind (blade-wake interactions) and results in strong fluctuations in C m [40,79,80]. The accurate predication of the flow complexities associated with the two phenomena, therefore, requires a fine azimuthal increment. The lines corresponding to dθ = 0.1°and 0.05°almost overlap for all the cases implying that the solution is no longer dependent on dθ. Fig. 6a shows how dθ affects the turbine power coefficient for different values of λ. For λ ≥ 4.5, it can be seen that dθ has limited effects on the turbine C P , though it can lead to a better prediction of the shaft's wake-blade interaction, as shown in C m curves in Fig. 5a and b.

Dependency on tip speed ratio
The increased deviation between the results of simulations using fine (0.05°-0.1°) and coarse (0.5°-1.0°) dθ with reducing λ (< 4.5) is also apparent from the predicted turbine C P shown in Fig. 6a. The results start to deviate at λ = 3.5. The deviation in terms of C P increases as λ reduces to 2.5 and then stays almost constant when λ further reduces to 1.5. This suggests that at 2.5 < λ ≤ 3.5, the separation is still growing on the blades by reducing λ, and there is a light stall. For 1.5 ≤ λ ≤ 2.5, the blades have already gone through a deep stall, and further reduction in λ from 2.5 to 1.5 might not bring more complexities in the flow. This is in line with the trend observed for C m in Fig. 5cf where the typical signatures of the dynamic stall (sudden drop in loads and follow-up fluctuations) are pronounced for 1.5 ≤ λ ≤ 2.5 ( Fig. 5e-f) but not present for 2.5 < λ ≤ 3.5 ( Fig. 5c and d). The figure also confirms that the predicted turbine C P using dθ of 0.1°and 0.05°a re the same for all λ values from 1.5 to 5.5. It can be concluded that for the solidity of 0.12, dθ = 0.1°is sufficient to capture the complex phenomena of dynamic stall on the turbine blades and further refinement in time step would not have any significant effects on the results. Fig. 6b shows the difference between the calculated C P obtained using different values of dθ and that obtained using the finest dθ. For λ > 3.5, it can be seen that for dθ = 0.5°, ΔC P is less than 1%. Therefore, an azimuthal increment of 0.5°is found to be the minimum requirement for moderate to high λ (> 3.5). The figure shows that ΔC P drops below 1% using dθ = 0.1°for all λ values, except for λ = 1.5, which is associated with the very low C P values at such λ. For λ = 1.5, this is further confirmed by the fact that ΔC P calculated using dθ of 0.05°and 0.025°does not show any further reduction. Based on the above discussions, the minimum requirements for azimuthal increment as a function of tip speed ratio is presented in Table 4. Two operating regimes are considered: (i) low to moderate λ and (ii) moderate to high λ. For the given turbine with a solidity of 0.12, the former regime corresponds to 1.5 ≤ λ ≤ 3.5 while the latter corresponds to 3.5 < λ ≤ 5.5. For low to moderate λ, the minimum requirement is dθ = 0.1°while for moderate to high λ this value can go up to dθ = 0.5°. Fig. 7 shows the instantaneous moment coefficient C m for the last turbine revolution versus azimuth θ for different solidities 0.06 ≤ σ ≤ 0.24 and azimuthal increments 0.05°≤ dθ ≤ 0.5°. Fig. 8 shows the predicted C P for different solidities calculated using different values of dθ. The tip speed ratio is 4.5, which is the same as the reference case. To reduce the computational cost, the azimuthal increment of 1.0°is no longer considered as the minimum requirement derived in Section 4.1 is 0.5°and dθ = 1.0°is thus considered to be too coarse. It can be seen that for 0.09 ≤ σ ≤ 0.24 (Fig. 7a-d), no significant difference is observed for different values of azimuthal increments. The C m values in the downwind near θ ≈ 270°show that as solidity decreases (airfoil chord-length decreases), the effect of the shaft's wake interaction with blades passing downstream becomes more pronounced. This is due to the increase in the relative length scale of vortices shed from the shaft to the airfoil chord. The more blade-wake interactions result in very small deviations in C m in the shaft wake  (θ ≈ 270°) between the lines corresponding to dθ = 0.5°and dθ < 0.5°f or σ = 0.09. However, this results in negligible (< 1%) difference in calculated C P between the simulations with dθ = 0.5°and the finest dθ (0.05°) for 0.09 ≤ σ ≤ 0.24 (Fig. 8). For a very low solidity of 0.06, significant differences in C m (Fig. 7e) and C P (Fig. 8) are observed for dθ = 0.5°and dθ < 0.5°. The difference in C m is mainly apparent in two regions: 90°< θ < 160°and 225°< θ < 300°where the former is associated with the azimuthal angles where dynamic stall might occur, and the latter corresponds to the region where the consequent blade-wake interactions would happen (see Section 4.1).

Dependency on solidity
Contours of instantaneous normalized lateral velocity at θ = 100°f or turbines with solidities of 0.06 and 0.09 are presented in Fig. 9. The promotion of trailing edge separation and increase of the laminar separation bubble length due to lower Re c for σ = 0.06 is clearly observed. This confirms that the observed deviation in C m and C P for σ = 0.06 is a result of the Reynolds number effect. Although at λ = 4.5 dynamic stall is avoided on the turbine blades (see Fig. 5b), the presence of strong Reynolds number effects could result in an earlier occurrence of stall on turbine blades at the same λ but lower chord-based Reynolds number Re c [81]. The Re c could be a result of the reduction in the airfoil chord length (which is equivalent to decrease in turbine solidity for a given n and D), which is also the case here. At λ = 4.5, Re c for σ = 0.12 (c = 0.06 m) is 170 × 10 3 while this reduces to 80 × 10 3 for σ = 0.06. Therefore, regardless of tip speed ratio and solidity, low Reynolds number regime (Re c < 10 5 ) would impose the finer minimum requirement of dθ = 0.1°to accurately predict the flow.
The impact of solidity on the minimum requirement for dθ is further investigated at a lower tip speed ratio of 2.5. The simulations are performed for dθ = 0.05°and 0.1°. The instantaneous moment coefficient C m for the last turbine revolution versus azimuth θ is shown in Fig. 10. The comparison shows insignificant differences in C P (< 1%) for dθ = 0.1°, which was shown to be the minimum requirement for dθ at such λ, and dθ = 0.05°. This means that in case the temporal resolution of the simulation is adequate to accurately predict the dynamic stall and blade-wake interactions, the minimum requirement for dθ is independent of solidity at a given Re c . Therefore, at a given Re c , tip speed ratio remains the dominant parameter to set the minimum requirement for dθ while a low Re c (< 10 5 ) might need the finer dθ of 0.1°. It should be noted that in the current study, solidity is modified by changing the blade chord length, however, modifying the solidity by changing the number of blades would most likely have negligible effect on the provided minimum requirements for dθ. This is because changing the number of blades would neither affect Re c nor the variations of angle of attack of the blades and consequently the dynamic stall behavior on blades while it only influences the frequency of occurrence of the blade-wake interactions. Nevertheless, more research is required to further investigate this.

Domain size: distance to the inlet
The operation of the turbine results in an induction field upstream of the turbine where the incoming flow is decelerated to lower velocities than the freestream. Rezaeiha et al. [6] showed that at λ = 4.5 and σ = 0.12 this induction field extends to more than 7.5D upstream. Therefore, the computational domain should extend far enough upstream not to cut this induction field. As the turbine induction field might change for different tip speed ratios and solidities, the dependency of the minimum requirement for the distance from the turbine center d i to the inlet to λ and σ needs to be investigated. Note that the minimum requirements for domain width and the diameter of the rotating core are considered to be unaffected by the tip speed ratio and solidity (see Section 2.2). Dependency on distance to the outlet will be investigated in Section 6. A list of the main computational parameters of the cases investigated in the present section is given in Table 5.

Dependency on tip speed ratio
The instantaneous moment coefficient C m for the last turbine revolution versus azimuth θ for tip speed ratios 2.5, 3.5 and 4.5 and for domains with different distance from the turbine center to the inlet 2.5D ≤ d i ≤ 15D is shown in Fig. 11. The tip speed ratios are selected as typical moderately low and moderately high values relevant for VAWTs. The turbine solidity is 0.12, which is the same as the reference case. It is seen that the simulations using domains with d i < 10D overestimate C m . This is more pronounced for higher tip speed ratios ( Fig. 11a-b), especially during the upwind quartile 45°≤ θ < 135°. This is associated with cutting the induction field, which results in    higher incoming velocity at the turbine incident and consequently higher C m on blades. For higher values of d i (10D and 15D), however, the C m distribution is almost the same for all tip speed ratios. Fig. 12 shows the difference between the predicted turbine C P for different tip speed ratios calculated using domains with different d i compared to the longest domain with d i = 15D. Substantial differences can be observed for domains with d i < 15D. The figure illustrates that |ΔC P | increases when λ becomes larger. This shows that the solution is more sensitive to the domain size at higher λ values. The non-monotonic value of |ΔC P | occurring for λ = 2.5 using domains with d i of 5D and 7.5D could be attributed to some error cancellations in C m curve. This suggests that comparison of the instantaneous values of C m is more reliable than the averaged values of C P for such an analysis.
The value of |ΔC P | for domains with d i = 10D with respect to the domains with d i = 15D is 1.9% and 1.3% for λ = 2.5 and 4.5, respectively. The above discussions show that the minimum requirement for the distance from the turbine center to the domain inlet d i in order to have |ΔC P | < 1% is 15D regardless of the turbine tip speed ratio.

Dependency on solidity
The instantaneous moment coefficient C m for the last turbine revolution versus azimuth θ for solidity of 0.12 and 0.24 calculated using domains with different distance from the turbine center to the inlet 5D ≤ d i ≤ 15D is shown in Fig. 13. The solidities are selected as two typical moderately low and moderately high values relevant for VAWTs. The tip speed ratio is 4.5, which is the same as the reference case. Note that to reduce the computational cost, the d i = 2.5D is no longer considered as the minimum distance as used in Section 5.1 as it was already shown to be too small. It can be seen that the simulations using domains with d i < 10D overestimate C m compared to the domains with d i ≥ 10D. This is more significant for solidity of 0.24 (Fig. 13b) and during the upwind quartile (45°≤ θ < 135°). Similarly, the observed overestimation of C m for domains with d i < 10D is attributed to cutting the turbine upstream induction field. For higher values of d i , the differences are negligible.   Fig. 14 shows the difference between the predicted turbine C P for different solidities calculated using domains with different d i values compared to the widest domain with d i = 15D. Significant differences can be observed for domains with d i ≤ 10D. This is in line with the trend observed for C m in Fig. 13. The figure illustrates that C m variations with domain size increases when σ becomes larger. This shows that the solution is more sensitive to the domain size at higher σ values.
The observations above conclude that the derived minimum requirement for d i in Section 5.1, i.e. d i = 15D, is valid regardless of the turbine solidity. This value will ensure that the predicted turbine performance, at different tip speed ratios and solidities, is minimally affected by d i and the domain sufficiently contains the induction field upstream of the turbine.

Domain size: distance to the outlet
The operation of the turbine results in the generation of a wake downstream of the turbine where the flow has lower velocity and higher turbulence intensity than the freestream. Rezaeiha et al. [6] showed that at λ = 4.5 and σ = 0.12, the turbine C P is underpredicted with more than 2% for domains with d o < 10D, where d o is the distance from the turbine center to the domain outlet. Therefore, the computational domain should extend far enough downstream to allow the turbine wake to sufficiently develop. As the turbine wake length might alter for different tip speed ratios and solidities [82], a systematic analysis of the minimum requirement for d o for different tip speed ratios and solidities is needed. A list of the main computational parameters of the cases investigated in the present section is given in Table 6. Fig. 15 presents the instantaneous moment coefficient C m for the last turbine revolution versus azimuth θ for tip speed ratios of 2.5, 3.5 and 4.5 obtained using domains with different distance from the turbine center to the outlet 6D ≤ d o ≤ 55D. Fig. 16 shows the difference between the C P for different tip speed ratios calculated using domains with different d o compared to the longest domain, i.e. d o = 55D. Note that two tip speed ratios are selected as typical moderately low and moderately high values relevant for VAWTs. The turbine solidity is 0.12, which is the same as the reference case. It can be seen that for all λ values the sensitivity of C m to d o is considerably less compared to that observed for d i . In addition, Fig. 15 illustrates that C m variations with domain size increases when λ becomes larger. In the downwind quartile (225°≤ θ < 315°), a small difference between the results of the domains with d o < 10D and domains with d o ≥ 10D is observed. This results in 2.2% deviation in the predicted C p using domains with d o = 6D and d o = 55D for λ = 4.5 (Fig. 16). Moreover, Fig. 16 also exhibits higher |ΔC P | values when λ becomes larger. This shows that the solution is more sensitive to the domain size at higher λ. Note that a similar trend was also observed for d i (see Section 5.1). For all the studied tip speed ratios, the difference between the predicted turbine C P using the domains with d o ≥ 10D and the domain with d o = 55D is less than 1%. Therefore, d o = 10D is considered as the minimum requirement for d o regardless of the tip speed ratio.   Contours of the normalized velocity magnitude illustrating the turbine wake for different tip speed ratios are shown in Fig. 17. The length of the turbine wake, defined as the length in which U/ U ∞ = 0.97, as a function of tip speed ratio (σ = 0.12) is presented in Fig. 18. The two figures imply that the CFD results could still be sufficiently accurate, though the computational domain might partially cut the turbine wake at x/D ≥ 10. This could be associated with strong velocity gradients occurring at x/D < 10, which is most particularly observed for higher λ values. This also explains the higher sensitivity of the predicted turbine C m (Fig. 15a) and |ΔC p | (Fig. 16) to d o for high λ values.

Dependency on tip speed ratio
The decrease in the turbine wake length with increasing λ was already observed for horizontal axis wind turbines [82]. The detailed analysis in Ref. [82] has shown that by increasing λ the distance between the shed vortices in the wake decreases, which leads to amplification of instabilities in the wake. This eventually leads to an earlier breakdown of the wake structure. Contours of instantaneous normalized lateral velocity for two tip speed ratio of 3.5 and 5.5 are shown in Fig. 19. The figure illustrates the trace of the blades' wake (dashed line) and the shed vortices from the blades on the leeward region traveling downstream (dotted circles). Two adjacent zones with different signs for lateral velocity are used to qualitatively identify the shed vortices. A comparison of Fig. 19a and b shows that the relative distance between the dashed lines and the center of the dotted circles has significantly increased for λ = 5.5 suggesting that the aforementioned mechanism is responsible for the reduction in the VAWT wake length by increasing λ.

Dependency on solidity
The instantaneous moment coefficient C m for the last turbine revolution versus azimuth θ for solidities of 0.12 and 0.24 calculated using domains with different distance from the turbine center to the outlet 6D ≤ d o ≤ 25D is shown in Fig. 20. The difference between the predicted turbine C P for different solidities calculated using domains with different d o compared to the longest domain with d o = 25D is shown in Fig. 21. The solidities are selected as two typical moderately low and moderately high values relevant for VAWTs. The tip speed ratio is 4.5, which is the same as the reference case. Note that to reduce the computational cost, the case with d o = 55D is no longer considered because the minimum requirement derived in Section 6.1 is 10D and d o = 55D is already predicted to be too large. It can be seen that there is an observable difference in C m between the values calculated using domains with d o < 10D and d o ≥ 10D while the lines corresponding to d o = 10D and 25D overlap. This difference results in 2.3% deviation between the calculated C P using domains with d o of 6D and 25D while this deviation is less than 1% for domains with d o of 10D and 25D (Fig. 21). Therefore, the already derived minimum requirement for d o (in Section 6.1) can also be employed regardless of the turbine solidity. This value (d o = 10D) will ensure that the predicted turbine performance, at different tip speed ratios and solidities, is not significantly affected by d o and the domain sufficiently allows the turbine wake to develop so that cutting the wake will have < 1% effect of the predicted turbine C P .
The turbine wake length as a function of solidity at λ = 4.5 is shown in Fig. 22. As already discussed in Section 6.1, partial cutting of the wake at d o ≥ 10D will have < 1% effect on the predicted turbine performance.
A similar explanation can be presented for the decrease in the wake length with increasing solidity. Increasing solidity (with increasing the blade chord length) will lead to increase in the length scale of the shed vortices in the turbine wake, which leads to the smaller relative distance between the vortices at the given λ. This will amplify the instabilities in the wake and eventually will result in the earlier breakdown of the wake structure. A comparison of the trends observed in Fig. 22 and Fig. 18 implies that the length of the turbine wake is more sensitive to tip speed ratio than the solidity.

Convergence criterion
For transient simulations, it is importance to initiate the data sampling only after the solution has reached a statistically steady state condition. For the case of wind turbines, it is customary to correlate the convergence of the unsteady simulation with the turbine revolution and to express it as the minimum number of turbine revolutions that the simulations needs to continue before reaching a statistically steady state condition. Fig. 23 schematically presents the definition of the statistical steady state condition employed in the present study. As the complexity of the flow for VAWTs varies with the operating and geometrical conditions, i.e. λ and σ, therefore, a sensitivity analysis is performed for the convergence criterion by applying systematic changes to the reference case. In this section, the sensitivity of the minimum requirements for the convergence criterion to different tip speed ratios and solidities are investigated. The details of the main computational parameters of the cases that are investigated are presented in Table 7.   The convergence of the results as a function of the turbine revolution is studied for different tip speed ratios. Fig. 24 shows the time history of C P during 50 turbine revolutions for two different tip speed ratio of 2.5 and 4.5. The tip speed ratios are selected as typical moderately low and moderately high values relevant for VAWTs. The turbine solidity is 0.12, which is the same as the reference case. It can be seen that after 20 turbine revolutions a statistically steady state condition is observed for all cases, as the C P variation with respect to 50th revolution has become < 2%. The difference in C P between two subsequent changes at 20 revolutions onwards is less than 0.2%. This is in line with the results of the earlier study by Rezaeiha et al. [6], which shows that 20-30 turbine revolutions is sufficient as the convergence criterion for a moderate tip speed ratio of 4.5 with a low solidity of 0.12. The results of the present study confirm that this convergence criterion can be generalized for other tip speed ratios. This suggests that 20-30 turbine revolutions are sufficient to allow the development of the flow over a VAWT so that statistically steady results are obtained irrespective of tip speed ratio. Note that at low tip speed ratios due to the inherent three-dimensional flow complexities associated to dynamic stall on blades and blade-wake interactions and the limitations of 2D modeling, the calculated C P might have some small oscillations between the subsequent turbine revolutions. Nevertheless, after approximately 20th revolution, the difference between the C P , for λ = 2.5, obtained after 20 revolutions and that after 50 revolutions remains below 2%.

Dependency on solidity
The convergence of the results as a function of the turbine revolution is studied for different solidities. Fig. 25 shows the time history of power coefficient C P during 50 turbine revolutions for two different solidities of 0.12 and 0.24. The solidities are selected as two typical moderately low and moderately high values relevant for VAWTs. The tip speed ratio is 4.5, which is the same as the reference case. It can be seen that after 20 turbine revolutions a statistically steady state solution is observed for all cases, as the C P variations with respect to 50th revolution is negligible (< 2%). The difference in C P between two subsequent changes at 20 revolutions onwards is less than 0.2%. The results of the present study confirm that this convergence criterion can be generalized for other solidities. This suggests that 20-30 turbine revolutions are sufficient to allow the development of the flow over a VAWT so that statistically steady state results are obtained for the two studied cases with moderately low and moderately high solidities.
The CFD simulations are based on two validation studies with wind tunnel measurements. Note that the presented validation study is performed using the very fine dθ of 0.1°, which is finer than (or equal to) the derived minimum requirements for all λ values.
The minimum requirements are derived based on limiting the effect of each studied parameter (azimuthal increment, distance from turbine center to domain inlet and outlet and number of turbine revolutions) on turbine C P to less than 1%. However, one should note that this does not mean that the total error of the CFD simulation associated with such computational parameters is 1%. A simple numerical error analysis    would reveal that the aforementioned errors can sum up to higher values which together with the other numerical errors; i.e. physical approximation error (e.g. 2D assumption as a representative of the turbine mid-plane), round-off error, iterative convergence error, spatial discretization error; would propagate to the final predicted turbine C P . Therefore, one should not underestimate the importance of such minimum requirements to minimize the numerical error.
On the other hand, the findings of this study are limited to URANS simulations. For scale-resolving simulations (SRS), such as hybrid RANS-LES and LES, different requirements, specifically for convergence criterion and azimuthal increment, would be necessary which demands further investigation. Moreover, for SRS, the difference between consecutive turbine revolutions would be significant which, therefore, C m and C P values would need to be averaged over multiple turbine revolutions to be a representative of the turbine averaged performance. The number of turbine revolutions for such averaging would require a dedicated study.
In the current study, solidity was modified by changing the blade chord length at a given number of blades (n = 2). Although the provided minimum requirements are most likely not affected by the number of blades, a detailed investigation of that is proposed for future research.
Moreover, chord-based Reynolds number Re c was also found to be influential for determination of dθ. Therefore, a dedicated investigation of the impact of Re c on the provided minimum requirements for dθ is proposed for future research.

Conclusions
A systematic sensitivity analysis is performed using URANS simulations to provide the minimum requirements for azimuthal increment, domain size and convergence criterion for CFD simulation of VAWTs at different tip speed ratios and solidities. The evaluation is based on validation with wind-tunnel measurements for two VAWTs. The following conclusions are obtained: (1) Azimuthal increment: The minimum requirement for azimuthal increment is significantly dependent on the flow regime over the turbine blades. As the variations of the angle of attack of the turbine blades are mainly driven by the tip speed ratio, this parameter is the dominant factor to set the minimum requirements for azimuthal increment. The derived minimum requirements for azimuthal increment are: • Low to moderate tip speed ratios: dθ = 0.1°• Moderate to high tip speed ratios: dθ = 0.5°T hese values will ensure that the flow complexities such as dynamic stall and blade-wake interactions are accurately predicted. For the studied turbine with a solidity of 0.12, the former regime corresponds to 1.5 ≤ λ ≤ 3.5 while the latter refers to 3.5 < λ ≤ 5.5. In general, the low to moderate tip speed ratio regime corresponds to the region where the variations of the angle of attack of the blades exceed the static stall angle and dynamic stall and the blade-wake interactions are thus present. Conversely, the moderate to high tip speed ratios correspond to the regime where the variations of the angle of attack of the blade are below static stall angle, and the flow is mostly attached. Therefore, one should have a prior knowledge of the following two important factors to distinguish between the aforementioned operating regimes of tip speed ratio: • Variations of the angle of attack for the operating regime of the turbine to be studied.
• Static stall angle of the employed airfoil on the turbine blade.
In addition, it is important to note that any change in geometrical and operational characteristics of the turbine which could induce a different flow separation behavior over the blades would require the finer dθ of 0.1°. The change in geometrical parameters includes, but not limited to: • Changing the airfoil shape • Changing the blade chord length (solidity) • Changing the blade surface roughness • Introducing pitch angle to blade • Adding passive/active flow control devices on blades.
The change in operating conditions includes operating the turbine at: • Low chord-based Reynolds number (Re c < 10 5 ) • Low freestream turbulence intensity (TI < 5%) Therefore, for any reason, if a user is not certain of the regime of tip speed ratio of the turbine, an azimuthal increment of 0.1°is always the safe choice to be employed for URANS simulations of VAWTs.
(2) Domain size: The minimum distance from the turbine center to the domain inlet needs to be 15D. This can be a safe choice to let the domain sufficiently contain the turbine upstream induction field regardless of the tip speed ratio and solidity.
The minimum distance from the turbine center to the domain outlet needs to be 10D. This value will ensure that the turbine wake is sufficiently developed inside the domain so that it allows accurate prediction of the turbine performance regardless of the tip speed ratio and solidity.
(3) Convergence: The minimum number of turbine revolutions to ensure that the simulations have reached a statistically steady state condition is 20-30 turbine revolutions. This convergence criterion can be employed regardless of tip speed ratio and solidity.
The aforementioned minimum requirements can be employed as guidelines to ensure the accuracy of CFD simulation of VAWTs.