Application of Richardson extrapolation method to the CFD simulation of vertical-axis wind turbines and analysis of the flow field

There is still discrepancy regarding the verification of CFD U-RANS simulations of vertical-axis wind turbines (VAWTs). In this work, the applicability of the Richardson extrapolation method to assess mesh convergence is studied for several points in the power curve of a VAWT. A 2D domain of the rotor is simulated with three different meshes, monitoring the turbine power coefficient as the convergence parameter. This method proves to be a straightforward procedure to assess convergence of VAWT simulations. Guidelines regarding the required mesh and temporal discretization levels are provided. Once the simulations are validated, the flow field at three characteristic tip-speed ratio values (2.5 - low, 4 - nominal and 5 - high) is analyzed, studying pressure, velocity, turbulent kinetic energy and vorticity fields. The results have revealed two main vortex shedding mechanisms, blade- and rotor-related. Vortex convection develops differently depending on the rotor zone (upwind, downwind, windward or leeward). Finally, insight into the loss of performance at off-design conditions is provided. Vortex shedding phenomena at the low tip-speed ratio explains the loss of performance of the turbine, whereas at the high tip-speed ratio, this performance loss may be ascribed to viscous effects and the rapid interaction between successive blade passings.


Introduction
Due to their higher power output, horizontal axis wind turbines (HAWTs) have been traditionally the preferred option against vertical axis wind turbines (VAWTs). Nevertheless, VAWTs are increasing their popularity, specially in urban areas, due to their particular characteristics (Tummala, Velamati, Sinha, Indraja, & Krishna, 2016). VAWTs are able to produce useful energy from lower wind speeds and they can work independently of the wind direction. As a result of their lower turbine rotational speeds, their noise levels are lower. In addition, their installation and maintenance operations are much simpler, as the bearings and generator unit may be placed on the ground. All these advantages show the need of further research on VAWTs to overcome the main obstacles for their implantation, which are basically their difficulty to self-start and the complexity of the flow developed through and around the turbine.
There are several methodologies for the prediction of the performance of VAWTs. Three main groups of methodologies may be identified: analytical models, CONTACT Andrés Meana-Fernández andresmf@uniovi.es based on fluid equations developed from the simplification of physical phenomena; numerical models, based on the discretized resolution of the Navier-Stokes equations; and experimental tests consisting in the fabrication and testing of a turbine prototype. In this work, the focus is set on numerical models, as they present better accuracy levels than analytical models. In addition, although experimental tests would be always desirable, the complete description of the flow field by experimental methods requires exhaustive measurement campaigns. At turbine design stages, it would be unpractical to build an experimental prototype for every turbine concept and test it. At these stages, the potential of computer fluid dynamics (CFD) simulations to obtain a complete picture of the flow field represents a great advantage (Akbarian et al., 2018;Ardabili et al., 2018;Mou, He, Zhao, & Chau, 2017).
Nevertheless, when performing CFD analyses, there is sometimes a lack of verification and validation of the employed numerical codes. Some solutions are considered as 'valid' without a proper assessment of mesh and temporal numerical convergence. Besides, the literature shows different procedures to verify the validity of the simulations. Firstly, some authors compare the results of their simulations with experimental results, considering their discretizations valid if they resemble the experiments. Sui, Lee, Huque, and Kommalapati (2015) studied the transitional effect on a turbulence model for a wind turbine blade and Lee, Min, Park, and Kim (2015) and Li et al. (2016) studied the performance of a VAWT, validating their simulations with experimental results. Yang, Guo, Zhang, Jinyama, and Li (2017) validated their study of the tip vortex shedding from a VAWT with experiments as well. Lam and Peng (2016) validated their 2D and 3D simulations to study the wake of a VAWT with PIV measurements, whereas Abdalrahman, Melek, and Lien (2017) compared their results with existing CFD and experimental benchmarks to validate their model for the study of the effect of the blade pitch angle on the performance of a VAWT. The second procedure to assess the validity of the simulations is the study of numerical convergence. This procedure is typically performed by refining the mesh until the solution no longer changes. Although this refinement is performed by several authors in the literature, there are discrepancies that sometimes result in an unnecessary waste of numerical resources. Bhargav, Kishore, and Laxman (2016) used two different meshes and two time step sizes to study the influence of fluctuating wind conditions on a VAWT. Li et al. (2018) studied the convergence of their simulations for the point of maximum power coefficient of the turbine with two different meshes as well. Bianchini, Balduzzi, Bachant, Ferrara, and Ferrari (2017) also used two meshes, but refining the grid only in the rotating part of the domain. Chen, Chen, Huang, and Hwang (2017) analyzed a VAWT at its design point using three meshes and comparing their results with simulations from other authors. Chen and Lian (2015) used three meshes as well to investigate the vortex dynamics in a VAWT. Meng, He, Wu, Zhao, and Guo (2016) also employed three meshes for the simulation of an offshore VAWT. Make and Vaz (2015) analyzed the scaling effects on offshore wind turbines with five different meshes. Balduzzi, Bianchini, Ferrara, and Ferrari (2016) used five different meshes and ten different time step sizes to propose dimensionless numbers for the assessment of mesh and temporal requirements for the CFD simulation of VAWTs. Jin, Wang, Ju, He, and Xie (2018) used seven different meshes to verify grid independence. Lin, Lin, Bai, and Wang (2016) also used seven different levels of refinement near the wall region of the blades to study the effect of modifications in the blade trailing edge on the performance of a VAWT. Wang, Cot, Adolphe, and Geoffroy (2017) recently studied the capacity of wind concentration over a roof using eight different meshes. Finally, some authors combine the verification of grid independence and the validation with experimental results. Subramanian et al. (2017) used two different meshes and time step sizes and validated their results experimentally. Wekesa, Wang, Wei, and Zhu (2016) studied the effect of turbulence on the aerodynamic performance of a VAWT, using two meshes for testing grid independence and validating the simulations with their own experiments. Lei et al. (2017) refined the grid close to the blades, generating three meshes and using the grid with the medium refinement because its results resembled more the experimental results. Qamar and Janajreh (2017) used three 2D VAWT meshes for the grid independence study and validated their simulations with experimental results. Tian, Mao, An, Zhang, and Wen (2017) used three meshes and validated experimentally the performance of a VAWT in the wake of moving vehicles. Réthoré, van der Laan, Troldborg, Zahle, and Sørensen (2014) verified and validated an actuator disk model for wind turbines using four different grid levels and validating the results experimentally. Marinić-Kragić, Vuina, and Milas (2018) performed 40 mesh modifications based on the combinations of different element sizing values along the blade and near the trailing edge. Then, after employing six different time step sizes, they validated their results experimentally. Finally, Rezaeiha, Kalkman, and Blocken (2017b) and Rezaeiha, Kalkman, Montazeri, and Blocken (2017) studied the effect of the shaft on the performance of a VAWT applying the Richardson extrapolation method and validating their results experimentally.
As it may be appreciated in the previous paragraphs, there is still a clear research gap with respect to the mesh and temporal requirements for the proper simulation of a VAWT, regardless of the availability (or not) of experimental data. As remarked by Lockard (2010), it is not unusual to find inconsistent and somewhat disappointing convergence properties on CFD codes, as shown by Vassberg and Jameson (2010). In this context, grid convergence analyses based on the Richardson extrapolation method have been used in different kinds of fluid dynamics related problems (Celik & Li, 2005;Marchi et al., 2016;Roache, 1998;Tengs, Storli, & Holst, 2018). Even some studies applying this methodology to the nominal working point of vertical axis wind turbines may be found in the literature (Almohammadi, Ingham, Ma, & Pourkashan, 2013;Rezaeiha et al., 2017b;Tingey & Ning, 2016;Zadeh, Komeili, & Paraschivoiu, 2014). Nevertheless, the wide operational ranges of a VAWT require the simulation of different tip-speed ratio (λ) values, so it might be difficult to accept that verifying and validating the numerical code for just the nominal λ is enough. In addition, due to the discrepancies between the number of mesh elements proposed by different authors in the literature and even the number of meshes and time step sizes used to assess the convergence of simulations, a straightforward procedure to assess numerical convergence in VAWT simulations becomes necessary for future research.
For all these reasons, the applicability of the Richardson extrapolation method to several points in the useful range of the power curve of a VAWT is studied in this work. A 2D domain of the rotor has been simulated using three different meshes and monitoring the power coefficient values of the turbine as the convergence parameter. Guidelines regarding the mesh near the airfoil boundary layers and the required temporal resolution are provided. Finally, once the simulations have been validated, the finest mesh has been used to perform an analysis of the flow field. The results regarding the flow field have revealed two main vortex shedding mechanisms in the rotor. In addition, insight has been given into the loss of performance of the VAWT when it works at off-design conditions.

VAWT geometry, computational domain and mesh
A 3-bladed low-solidity VAWT with DU 06-W-200 airfoils with a radius of 0.5 m was simulated in a 2D domain, as represented in Figure 2. The DU-06-W-200 airfoil, designed in Delft University, is allegedly supposed to present a good self-starting behavior (Claessens, 2006). Being a straight-bladed turbine, the main tridimensional aerodynamic effects are the blade tip effects, which are detrimental to the turbine performance. Nevertheless, with a sufficient turbine aspect ratio (height/radius), a 2D simulation of the turbine mid-plane is a reasonable choice. In fact, due to the relatively great number of cells required to model the full 3D turbine, most of the analyses in the literature are performed using 2D simulations. A range of tip-speed ratio values between 2 and 5.5 has been selected to compare the results at the nominal working point with two offset points, one in the dynamic stall region (λ = 2.5) and other with the turbine producing a high flow blockage (λ = 5) (Amet, Maitre, Pellone, & Achard, 2009). The details of the turbine, in line with optimum parameters reported in the literature (Meana-Fernández, Solís-Gallego, Oro, Díaz, & Velarde-Suárez, 2018), are collected in Figure 1. A rotor tower has not been included in the design and it has not been modeled, as not all VAWT design concepts include it. If it were desired, the grid requirements for the boundary layers of the blade airfoils may be easily extrapolated to the tower. Figure 2 shows the different regions of the computational domain. An interface between the turbine rotating zone and the fixed one has been placed at a distance of 5.5D (rotor diameters) from the rotor, in order to avoid distortion effects in the transfer of information between adjacent fluid regions. The total size of the circular computational domain is 12D. These values are in agreement with the values reported in the literature (Alaimo, Esposito, Messineo, Orlando, & Tumino, 2015;Almohammadi et al., 2013;Mohamed, Ali, & Hafiz, 2015;Rezaeiha, Kalkman, & Blocken, 2017a;Zadeh et al., 2014).
Finally, in order to perform a mesh convergence study, three different meshes with different levels of discretization have been generated using the software GAMBIT R . The grids have been generated with triangular elements; however, quadrilateral elements have been employed in the near-wall regions (distance to the wall < 2 mm) for a better control of the mesh growth. The position of the first mesh node has been calculated to ensure a y + value less than 1. A view of the mesh around the airfoils is depicted in Figure 2, bottom. In summary, Mesh #1 (fine mesh) has 989,770 elements, Mesh #2 (medium mesh)

Numerical solver
The incompressible unsteady Reynolds-Averaged Navier-Stokes equations have been solved with the commercial package ANSYS FLUENT R , using the k-ω-SST model Menter (1994) for the closure of turbulence, which combines the advantages of k-ε and k-ω models in predicting aerodynamic flows, and in particular in predicting boundary layers under strong adverse pressure gradients (Argyropoulos & Markatos, 2015). The boundary conditions applied may be identified in Figure 2. A constant inlet velocity condition of 9 m/s, a typical design value for VAWTs, and an outlet pressure equal to atmospheric pressure have been set. The domain size, larger than 10 times the turbine diameter, ensures that the boundary conditions will not interfere in the flow developed inside the rotor (Alaimo et al., 2015;Almohammadi et al., 2013;Mohamed et al., 2015;Rezaeiha et al., 2017a;Zadeh et al., 2014). The sliding mesh technique has been applied thanks to the previously defined interface between the rotor and outer domain zones. Finally, the wall boundary condition has been applied to the rotor blades. The time step used to perform the simulations is discussed in the following section.

Grid and temporal convergence analysis
Before proceeding to the analysis of the simulation results, a grid convergence analysis was performed by means of the Richardson extrapolation method (Richardson & Gaunt, 1927). The Richardson extrapolation method, also known as 'h 2 extrapolation', 'the deferred approach to the limit' or 'iterated extrapolation', is a method for obtaining a higher-order estimate of the continuum value (value at zero grid spacing) from a series of lower-order discrete values. As introduced by Roache (1997), a numerical simulation yields a quantity f that can be expressed as: where h is the grid spacing of the simulation. Functions g i are defined in the continuum and thus are independent of the grid spacing. f exact is the continuum value at zero grid spacing.
As stated in Roache (1997), for a second-order method (g 1 = 0), by combining the results f 1 and f 2 from two different grids of spacing h 2 (coarse) and h 1 (fine) and neglecting third-and higher-order terms, an estimate for f exact may be obtained, resulting into the original statement of Richardson and Gaunt (1927) for h 2 extrapolation: Defining the grid refinement ratio as: Equation (2) may be rewritten as: This equation may be generalized to pth order methods (Roache, 1998) as: This method allows not only the estimation of the continuum value, but also provides a practical estimation of the grid refinement error due to the discretization of the simulation domain. For the practical application of the method, the guidelines proposed in ASME (2008) have been followed: • First of all, a representative grid size parameter, relating the cell size and the number of cells, has been defined as follows:  Roache (1998).
In this work, the variable φ used to judge grid convergence is the mean-time power coefficient of the turbine, which has been monitored from the different simulations performed.
• Afterwards, the apparent order of the method has been calculated as: where ε 32 = φ 3 − φ 2 , ε 21 = φ 2 − φ 1 are the absolute errors of the variable of interest φ obtained with the three different meshes. q(p) is a function depending on the refinement ratios between the meshes and the behavior of the solutions obtained as the grid is refined, defined as: where s is the parameter related to a monotonic or oscillatory behavior of the solution as the grid is refined: Negative values of s are an indication of oscillatory convergence. In addition, if either ε 32 or ε 21 is very close to zero, the above procedure does not necessarily work. This breakdown may be ascribed to oscillatory convergence or, in some cases, it could mean that the exact solution has been already attained (Roache, 1997). As it may be observed, Equations (7) and (8) must be solved iteratively. • Once the exact values of p, q(p) and s were obtained, the extrapolated value of the solution (the estimator of the exact solution) has been calculated in a similar way to Equation (5): • Then, the error estimates for the relative error and extrapolated relative error have been obtained as: and • Finally, the grid convergence index (GCI) has been used as an indicator of the mesh convergence level: being F S a security factor for the calculation of this index, which may be set to 1.25 when having three different meshes (Roache, 1998).  As an additional step, if a certain grid convergence level GCI * were desired, the required grid resolution r * with respect to the finer mesh might be obtained as: An example of the values obtained with this method for the three working points with maximum performance, λ = 3.5, 4 and 4.5, is displayed in Table 2. For every case, it was observed that the mean value of the power coefficient of the turbine attained convergence after 6 rotor revolutions, as shown in Figure 3. Thus, the values presented in this work correspond to the seventh revolution of the rotor. Although this value is in discrepancy with results from other authors (values over 20 revolutions have been reported in the literature Rezaeiha et al., 2017b), the necessity of a practical model determined this choice (computational times are around 1 week in a 4-nodes Intel Core i7-52820K at 3.3 GHz and 64 Gb RAM). In addition, the extrapolation towards t → ∞ of the mean C P value using an exponential function gave practically the same value as the value monitored in the last time step. This result is interesting regarding industrial purposes, as it is possible to save simulation time performing just a few rotor revolutions and then extrapolating the result towards infinity if only the C P value is sought.
The results of the Richardson extrapolation method are collected in Figure 4. This figure shows the different values of the power coefficient obtained with the three meshes at different tip-speed ratio values and the extrapolation performed towards h → 0, which would correspond to a mesh composed of an infinite number of elements. It may be observed that, at values near the nominal working point of the turbine λ = 3.5 − 4.5, the mesh convergence levels are really good (below 2% and even 0.043% at the nominal working point). Outside these region, the mesh convergence study is not so good. As previously introduced, if the differences between the values of the magnitudes used to study mesh convergence are very small, the procedure might not reach a trustworthy solution (either there is oscillatory convergence and/or an exact solution has already been attained). The breakdown of the method is clearly happening at a low tip-speed ratio (λ = 2.5), where the oscillatory convergence is clearly observed in Figure 4. This result is consistent with the issues found by Celik and Li (2005) when extrapolating their cases with non-monotonic convergence. Hence, a new simulation with a mesh one level finer than Mesh #1 was performed. The extrapolation method was applied again with this new mesh and meshes #1 and #2 from the previous analysis, as shown in Figure 4 (λ = 2.5). It may be observed that the method converges, avoiding oscillatory behavior and providing a more reasonable value for the power coefficient. On the other hand, the thin attached boundary layers that arise at the highest tip-speed ratios (λ = 5 and 5.5) and that will be shown in the next section reveal the importance of the refinement of the mesh in the boundary layer regions in order to model correctly the flow behavior. Figure 5 shows the CFD prediction of the power coefficient of the turbine with the three different meshes alongside the power curve predicted by an analytical double-multiple streamtube model developed by the authors. Both methods predict the maximum turbine performance at the same tip-speed ratio. Despite the assumptions performed by streamtube models about the nature of the flow (wake not modeled, downwind zone of the rotor assumed to be in the fully expanded wake of the upwind zone, influence of the downwind zone on the upwind zone neglected), the slopes of the predicted power curves with the streamtube model and the finest mesh match. In addition, some comments about the results of the Richardson extrapolation method presented before may be introduced. At first sight, for the lower tip-speed ratios, it may be appreciated that the difference in the power coefficient predicted by the three meshes is almost negligible. By looking at Figure 5 and from the results of the extra analysis with a finer mesh commented before, it may be assured that an accurate enough solution had been already found with the coarsest mesh, being this fact the cause of the breakdown of the extrapolation method for λ = 2.5. For the higher tipspeed ratio values, however, the differences between the power coefficient values of the three different meshes are higher. As previously commented, this might be indicating that the number of cells in the cross-streamwise direction to the airfoil in the coarser meshes is too low.
Nevertheless, as the turbine is unlikely going to work at tip-speed ratio values higher than the nominal one, it may be considered that the fine mesh is accurate enough for a practical study of the turbine.
For further validation of the accuracy of the procedure and the results of the finest mesh, the evolution of the streamwise and cross-streamwise forces on an airfoil during a whole turbine rotation is shown in Figure 6. Values for three different tip-speed ratios are shown (λ = 2.5 (low), 4 (nominal) and 5 (high)). This figure shows also the evolution of the aerodynamic torque of the turbine. The results for all the variables show that convergence is very good for the low and nominal tip-speed ratios, where it becomes difficult to distinguish between the results of the three meshes during the whole rotation cycle. For the highest tip-speed ratio, even when it is possible to differentiate the evolution of the curves for the three different meshes (especially for the aerodynamic torque), the agreement between the results obtained is also relatively good. Therefore, it may be considered that the discretization of the finest mesh is enough to model accurately the flow behavior. Regarding the low tip-speed ratio λ = 2.5, the force oscillations present in the downwind zone of the turbine will certainly cause oscillations in the produced aerodynamic power. This effect must be considered when designing/selecting the generator and electronic control system for the turbine.
Finally, in order to select an adequate temporal resolution, the finest mesh has been used to perform a convergence study of the time step size. Three discretization levels have been tested, corresponding to rotor advancements of 1 • , 0.5 • and 0.25 • per time step. Following the guidelines for the Richardson extrapolation method, a so-called spatial-temporal h ST index has been defined as:  being N s the number of cells, N t the number of time steps per rotor revolution, A i the size of the cell i and ω the rotational speed of the turbine. Figure 7 shows the results of this extrapolation method, which results in a T-GCI (Temporal-Grid Convergence Index) of 0.05% for the medium-to-smallest time step at λ = 4. Hence, all the calculations presented in this study have been performed using a temporal discretization corresponding to a rotor advance of 0.25 • per time step, that is, 1440 time steps per rotor revolution.

Comparison with existing benchmarks
In an attempt to compare the results of this study with the existing literature, a small benchmark has been developed with results from other authors, trying to find VAWTs with similar values of solidity and similar blade airfoils (Bedon, Betta, & Benini, 2016;Delafin, Nishino, Wang, & Kolios, 2016;Sabaeifard, Razzaghi, & Forouzandeh, 2012). Table 3 collects the main characteristics of each VAWT, whilst Figure 8 shows the comparison between the optimal tip-speed ratio of the turbines depending on their solidity. It may be appreciated that the maximum power coefficient is attained at similar values of the tip-speed ratio in the case of turbines with similar solidity values (Bedon et al., 2016;Delafin et al., 2016) and that the turbine studied in this work, with the same airfoil and a similar freestream Reynolds number to the turbine presented in Sabaeifard et al. (2012), fits into the trend of the values obtained by these authors. The turbine from Delafin et al. (2016) Figure 8. Comparison of the optimal tip-speed ratio of the studied turbine with existing results for similar turbines. attains a higher power coefficient than the rest of the turbines compared. It employs a NACA 0015 airfoil, which has shown a higher efficiency than thinner or thicker airfoils from its family (such as NACA 0018) Gosselin, Dumas, & Boudreau, 2013;Meana-Fernández et al., 2018) and has been modeled using a vortex model, which might overpredict results when compared to CFD simulations. All these facts may explain why this turbine presents higher power coefficient values. The differences between the other turbines may be either ascribed to the differences in the Reynolds numbers (it has been observed that an increase in the Reynolds number shifts the power curve up-and leftwards Bausas & Danao, 2015;Meana-Fernández et al., 2018;Paraschivoiu, 2002) and the airfoil used to build the blades. Despite the differences observed between the curves, it may be considered that the results of this study are consistent with the results found in the literature.
In order to provide an additional source of validation, CFD results from Sabaeifard et al. for different VAWT solidities have been plotted in Figure 9 and compared with the power curve of the turbine simulated in this study. It may be appreciated that the turbine fits perfectly into the graph regarding both the shape of the power curve and the optimal tip-speed ratio at which the maximum power coefficient is attained. The smaller value of the peak coefficient may be ascribed to the lower Reynolds number of the turbine of this study, as previously commented.  Figure 9. Position of the power curve of the studied turbine with respect to its solidity.

Numerical description of the flow field
Once the CFD methodology has been verified and compared with existing results from other authors, it has been used to analyze the flow behavior around the turbine at three different tip-speed ratio values: λ = 2.5, corresponding with a working point in the vortex shedding predominant region; λ = 4, corresponding with the nominal working point; and λ = 5, corresponding with a working point past the nominal working point, with highly attached boundary layers and predominant viscous effects. The pressure, velocity, turbulent kinetic energy and vorticity fields have been analyzed.

Pressure on the turbine blades
Figures 10-12 show the chordwise distribution of the pressure coefficient on a blade for λ = 2.5, λ = 4 and λ = 5 respectively. The strong pressure gradients observed corroborate the selection of the k-ω-SST model for turbulence modeling (Alaimo et al., 2015). The maximum pressure differences in the distribution arise in the upwind part of the turbine, the first stage of power extraction from the wind. In addition, differences between the windward and leeward regions of the turbine may be identified, with the blades in the windward zones receiving greater pressure from the incoming wind. The blades have been found to exhibit greater pressures in the region between 30 • and 150 • , which correspond to one-third of the whole rotating path. This fact may explain why 3-bladed turbines are predominant, as with this design only one blade remains in the high-loading zone at every instant. Regarding the results, the worst position of the blade seems to be 60 • . Finally, the effects of the increase of the rotational speed of the turbine are easily identified. At λ = 2.5, the pressure coefficient differences between the pressure and suction sides of the airfoil are only slightly higher in the windward region. The rest of regions do not suffer great shifts in the pressure distribution along the blades, which may be correlated with the low power extraction capacity of the turbine at this operational point. At λ = 4, the pressure coefficient distribution varies accordingly to the blade position, and it may be concluded that significant power is being extracted at several positions during the rotor cycle. Finally, at λ = 5, the pressure coefficient values increase with respect to the nominal working point, but the lower power extraction at this operational point allows to conclude that this rise in pressure does not translate into an effective power generation (maybe due to higher viscous effects on the airfoil surface that arise as a consequence of the higher rotational speed).

Velocity field
Contours of normalized velocity V/V ∞ for a whole rotor cycle are shown in Figure 13 for the three tip-speed ratio values studied. At first sight, a wake of size D is clearly identifiable, similar to the one shed by a cylinder with the same diameter as the turbine. The vortices shed by the blades at λ = 2.5 are also easy to identify. It may be also appreciated that the blockage produced by the rotor in the wind current is very low; hence, it is not surprising that the turbine is not capable of extracting enough power at this operational point. Comparing the contour values at λ = 4 and 5, the greater blockage effect of the rotor at high tip-speed ratio values may be confirmed. This increase on the blockage, however, does not translate into a greater power extraction; hence, the blockage level at λ = 4 seems to be the optimal one to maximize energy harvesting. Figure 14 shows the contours of normalized turbulent kinetic energy (TKE) for λ = 2.5, 4 and 5. This type of contour maps is very useful to identify the main zones of turbulence production and its further convection downstream. The turbulent kinetic energy is defined as:

Turbulent kinetic energy
where u i are the turbulent velocity fluctuations in the flow. Vortex shedding phenomena are clearly visible in the case of λ = 2.5, which start to roll up behind the airfoil before being convected downstream. Leaving aside the obvious differences between the upwind and downwind parts of the rotor (in the upwind part the wind comes "cleanly" onto the blades, whereas in the downwind part blade performance is clearly affected by the   wakes shed from the upwind part), a totally different behavior is present between the windward and leeward zones of the rotor. In the windward zone, the blades move towards the wind at higher relative velocities and generate much narrower wakes. In addition, the air flowing through the turbine convects the vortices before they can start to roll up, resulting in a much cleaner flow pattern such that it is very easy to identify the path followed by the blades during rotation. On the other hand, on the leeward side of the rotor, vortices are shed from the blades and roll up before being convected. This results in more intricate flow patterns, where the new vortices shed from the incoming blades are convected across a sea of vortices shed by the previous blades. The flow becomes blurry, nevertheless, it is still possible to associate each wake to its corresponding blade. This difference between leeward and windward regions is really interesting and can be already detected in Figure 6. In fact, it is the explanation for the oscillatory behavior of the forces on the blade that is appreciated in Figure 6, top (see angles in the downwind-leeward region). Finally, despite the apparent simplicity of the turbulent kinetic energy contours at λ = 4 and 5, there are still some comments that may be made from the comparison between the contours of Figure 14 to try to give insight into the reason why the turbine is more effective at the lower tip-speed ratio. Taking a closer view to the downwind-leeward region, when a blade enters this region it must pass across the turbulent kinetic energy traces left from the previous ones. In the case of λ = 4, four traces, whereas for λ = 5, six traces. This must be translated into a loss of efficiency in the performance of the airfoil in this region. This effect is added to the increasing viscous forces due to the rotational speed of the turbine, and may contribute to explain why further increases in the tip-speed ratio are no longer beneficial to the turbine performance. In addition, on the windward side of the rotor, the 'turbulent kinetic energy wake' is much longer at λ = 5 than at λ = 4. This could be attributed to a poorer airfoil performance in that region, as the turbine rotates so quickly that there is no time for the wakes shed on this side of the turbine to be drifted away from the rotor by the wind before the next blade comes and sheds its own wake. In addition, the high tip-speed ratio of the turbine prevents the airfoil from reaching angles of attack high enough to ensure significant lift generation.

Vorticity
Finally, contours of normalized in-plane vorticity are shown in Figure 15. Vorticity, which is an indicator of the tendency of the fluid to rotate, has been made dimensionless with the rotational speed of the turbine ω r and it is defined by: where u i is the velocity in the i direction.
Vorticity contours show two different mechanisms of vortex shedding and wake development. The first one, of greater magnitude, is vortex shedding from the blades during rotation. Two opposite-sign high-vorticity regions may be identified after each blade trailing edge. The second one arises from the combination of the turbine rotation and the incoming wind velocity. The vortices shed from the blades during their rotation are convected downstream by the incoming wind, forming different patterns depending on the tip-speed ratio. The combination of the turbine rotation and the wind velocity, thus, generates a wake of a size around the turbine diameter D. This wake presents lower values of vorticity, as the vortex mixing process starts as soon as the vortices leave the blades and makes vorticity values drop. At λ = 2.5, the turbine wake is full of unsteadiness due to the retardation in the dissipation of the big vortices shed from the turbine blades. Nevertheless, two different regions with opposite vorticity signs on the leeward and windward zones of the rotor may be identified. Looking at the contours for λ = 4 and 5, the distinction between leeward and windward regions becomes even more evident, with the wake divided in two halves of opposite vorticity. The size of the wake D shed by the whole rotor and its characteristics resemble the wake shed by a cylinder of the same diameter. Following this line of thinking, if VAWT farms were to be simulated, it would be possible to save computational costs by considering them as cylinders if they are working at high enough rotational speeds. Finally, comparing the contours at λ = 4 and 5, the greater number of shifts in the sign of vorticity at λ = 5 may be correlated with a greater generation of vortices. And, although the boundary layers are more attached to the blades, the interaction between their successive passings contributes to generate greater levels of unsteadiness inside the rotor that finally translate into performance losses.

Conclusions
The applicability of the Richardson extrapolation method to the CFD simulation of vertical-axis wind turbines has been verified. A 3-bladed low-solidity VAWT with DU 06-W-200 airfoils was simulated in a 2D domain using U-RANS k − ω SST modeling. A convergence study applying this method to three meshes with different levels of discretization determines the level of uncertainty of the final mesh. When convergence problems arise, an extra analysis with a finer mesh is enough to ensure convergence of the method and the adequate selection of the mesh. It has been concluded that performing the grid convergence analysis at only one point of the turbine power, as typically done in the literature, does not guarantee the same level of accuracy for the whole curve, especially at high tip-speed ratios. Specifically, a discretization of 40 cells in the first 2 mm in the cross-streamwise direction from the airfoil wall and 12 cells/mm in the streamwise direction for the airfoil chord is enough to capture all the relevant fluid phenomena. Additionally, a new spatial temporal h-index for the assessment of the temporal discretization may be defined using the Richardson extrapolation method as well. It was concluded that a rotor advancement of 0.25 • per time step is adequate enough. With these conditions, reasonable agreement was found between the results of this work and existing benchmarks.
The flow behavior (pressure, velocity, turbulent kinetic energy and vorticity) of the turbine shows two main vortex shedding mechanisms: vortex shedding from the blades during rotation and interaction of the turbine rotation with the incoming wind. Vortex convection develops differently depending on the rotor zone (upwind, downwind, windward or leeward), but finally generating downstream a wake of the size of the turbine diameter. The performance of the turbine depends on the tip-speed ratio, with great vortex shedding phenomena at small tip-speed ratios being responsible for the poor turbine performance. At high tip-speed ratios, the combination of viscous effects in the boundary layers and the increased interaction between successive blade passings also decrease the turbine performance.
This study does not consider tridimensional effects, so future work should include the realization of 2.5D or 3D full simulations using the discretization values from the finest mesh. Additionally, experimental tests on the prototype to determine the power curve of the turbine would be helpful to provide a deeper validation of the developed simulations.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This work has been supported by the 'FPU' predoctoral research scholarship provided by the Spanish Ministry of Education, Culture and Sports. The authors also want to acknowledge the support from the Projects 'Desarrollo de una herramienta de diseño optimizado de perfiles aerodinámicos para su utilización en turbinas eólicas de eje vertical' from the University Institute of Industrial Technology of Asturias, financed by the City Council of Gijón, Spain and 'Diseño optimizado de una turbina eólica de eje vertical' from the University of Oviedo Foundation, financed by the company AST Ingeniería and 'Desarrollo y construcción de turbinas eólicas de eje vertical para entornos urbanos' (ENE2017-89965-P) from the Spanish Ministry of Economy, Industry and Competitiveness.