The Effects of Airfoil Thickness on Dynamic Stall Characteristics of High‐Solidity Vertical Axis Wind Turbines

The flow physics of high solidity vertical axis wind turbines (VAWTs) is influenced by the dynamic stall effects. The present study is aimed at investigating the effects of airfoil thickness on the unsteady characteristics of high solidity VAWTs. Seven different national advisory committee for aeronautics (NACA) airfoils (0008, 0012, 0018, 0021, 0025, 0030, 0040) are investigated. A high fidelity computational fluid dynamics (CFD) approach is used to examine the load and flow characteristics in detail. Before the study is undertaken, the CFD simulation is validated with experimental data as well as large eddy simulation results with sound agreement. The investigation demonstrates that increasing the airfoil thickness is actually beneficial not only for suppressing the dynamic stall effects but also to improve the performance of high solidity turbines. Interestingly this is accompanied by a slight reduction in thrust component. The strength and radius of the dynamic stall vortex decrease with increasing airfoil thickness. The airfoil thickness strongly influences the pressure distributions during dynamic stall process, which is driven by the suction peak near the leading edge. The knowledge gained might be used by blade engineers for designing future turbines and for improving the accuracy of engineering models.


Introduction
The increased fossil fuel depletion in past decades has led to a shortage in electricity production causing the increased demand to sustainable power sources. Wind energy has been identified as one of the most promising sources for renewable energy industry. There are two main categories of wind turbines according to its rotational axis: Horizontal axis wind turbines (HAWTs) and vertical axis wind turbines (VAWTs). The former has been developed significantly during the last decades in terms of research and commercialization. Unlike HAWT, since its first proposal by

DOI: 10.1002/adts.202000204
Darrieus [1] back in 1931, lift driven VAWT has faced a hard stage in its development phase. The main factors that influenced the discontinued development of VAWT are the low tip speed ratio and the difficulty in controlling the rotor speed, [2] not to mention the low starting torque. However, VAWT offers some potential in terms of smallscale and domestic applications. [3] In this range of application, the mounting system of VAWT seems not to be as complex as HAWT and therefore the mounting cost can be reduced, since the heavy generator equipment can be mounted on the ground.
Not only the development from the construction and structural sides, studies concerning the aerodynamic property of VAWT have increased considerably in the last years. The fact that VAWT experiences dynamic stall has become the huge challenge for improving the aerodynamic performance and the structural stability of the turbine. Dynamic stall is a common phenomenon that appears mainly on rotating rotor blades and it becomes one of the most limiting factors of the aerodynamic performance. [4] It is also likely unpreventable during the operation of VAWT at low tip speed ( ≤ 3). [5] A lot of studies have been carried out through various experiments and simulations. [6][7][8][9][10][11] The rotor blades of VAWT have a critical role to convert lift into power. Due to the appearance of the dynamic stall, the rotor blades experience a drop in lift and simultaneously an increase in pressure drag, both of which reduce the rotor torque. [12] Because of this phenomenon, the rotor blade has to withstand stronger loads which consequently can reduce the turbine lifetime but also potentially its performance. Therefore a robust design of rotor blade is required. [13] The topic of airfoil design is never far from discussing the airfoil geometry in which special attention is paid to the profile thickness. Fundamentally, different stall properties are expected from airfoils with different thickness. [14] At lower Reynolds number, thicker profiles are preferable. [15] Danao et al. [16] have found an improvement in the overall performance of the vertical axis wind turbine when using an airfoil with slight camber such as the LS0421, in contrast to the less favorable results when using the NACA5522 airfoil which has a 5% camber. Furthermore, it was found that higher values of torque in both the downwind and upwind regions are generated by blades with a camber along the blade path, whereas profiles with an inverted camber produce power mainly in the upwind region. [16] Aside from the camber and relative thickness, it has been found that a large leading edge radius and a sharp trailing edge also encourage optimum performance. [17] Islam et al. [17] studied various symmetric NACA 4-digit airfoils and found that they are not suitable for application in VAWTs and neither is asymmetric airfoils designed with aviation purposes in mind. Wang et al. [18] evaluated several NACA airfoils applied on medium solidity VAWTs. Their study showed that the power level depends on the airfoil thickness. They recommended the airfoil thickness to be around 15-18%. The optimized airfoil effect on VAWT performance was studied by Liang and Li [19] based on the NACA 15% airfoil. Very recently, Rogowski et al. [20] investigated nine different NACA airfoils on VAWT performance, but their work was focused only on large turbines operating at low to medium angles of attack. In this sense, the performance is driven solely by the aerodynamic characteristics of the airfoil itself and the vortex dynamics is not of importance. Several other authors have made contributions to the studies concerning the airfoil effects on VAWT within the last 3 years, for examples Song et al. [21] , Liang and Li, [22] Mazarbhuiya et al., [23] and Rezaeihaa et al. [24] However, most of the studies were focused on low to medium solidity turbines and on the rotor power.
Because there is no generally applicable rule when designing the airfoil for VAWT applications, research on the influence of airfoil thickness on its aerodynamic properties is increasing. Furthermore, past studies were mostly dedicated to the performance evaluations but detailed analyses on the dynamic stall characteristics on high solidity vertical axis wind turbines are rarely found. The present study is aimed at investigating the effects of airfoil thickness on the dynamic stall characteristics on high solidity vertical axis wind turbines. Seven airfoils with varying relative thickness from the NACA airfoil family are chosen for the present investigations. The main novelties of the paper lie on the detailed evaluations of the airfoil thickness effects for high solidity VAWTs associated with the flow physics and quantification of the dynamic stall vortex in the wake. To the best of the authors' knowledge, this aspect has not been done in the past. The paper is structured as follows. The numerical setup and simulation strategy are presented in Section 2. The simulation approach will be validated against experimental data and other numerical results in Section 3.1. Then, detailed examinations on the effects of airfoil thickness in terms of integrated and azimuthal loads will be given in Section 3.2. Section 3.3 discusses the evolution of the dynamic stall vortex for each airfoil and presents the quantifica- tion of the vortex strength. To gain further insights, the evolution of the pressure coefficient will be given in Section 3.4. At last, all the results and discussion will be concluded in Section 4.

Studied Turbine Geometry
The investigated rotor consists of three blades employing a symmetrical NACA airfoil cross-section rotating in counter-clockwise direction. Seven different relative thicknesses are assessed in the present paper employing the NACA 0008 (t∕c = 8%), NACA 0012 (t∕c = 12%), NACA 0018 (t∕c = 18%),NACA 0021 (t∕c = 21%),NACA 0025 (t∕c = 25%), NACA 0030 (t∕c = 30%) and NACA 0040 (t∕c = 40%), where t is airfoil thickness, (m) and c is blade chord, (m). Figure 1 illustrates the geometry of the airfoil section. Note that the 2Dsimulations represent the mid-span section of the blade where the occurrence of the spanwise flow is neglected. The full 3D results may provide better accuracy (especially concerning the tip loss and spanwise flow effects) but running such large computation for many test cases is proven to be difficult. Studies have demonstrated that 2D computations are able to provide reliable results in comparison with measurement data with reasonable computational expenses. [25][26][27][28] The rotor has a radius of 0.4 m, a chord length of 0.2 m with a zero pitch angle. The turbine has a very high rotor solidity of = Nc∕R = 1.5. This makes the turbine suitable for the low tip speed ratio regime with a good starting torque. The general overview of the turbine is listed in Table 1. In this study, a constant wind speed of 8 m s −1 was applied, chosen in accordance with the measurement campaign used for validating the CFD setup. [29] The coordinate system (CS) of studied VAWT is defined as shown in Figure 2.
www.advancedsciencenews.com www.advtheorysimul.com The inertial CS defines the global non-rotating system of the VAWT. Another system is the local CS on the blades which rotates around the rotor shaft as its center point. The angle of attack of the blade is defined as the angle between the relative velocity vector V R and the chord line.

Computational Mesh and Numerical Setup
The computational fluid dynamics approach is a method to model the physical phenomena of the fluid flow and can be seen as a combination of interdisciplinary fields such as physics, numerical mathematics, and computer sciences. For flow simulation with CFD, turbulence is the limiting factor for accuracy and applicability. One way to simplify the solution variables of the Navier-Stokes equations is given by the Reynolds-averaged Navier-Stokes (RANS) approach, which pose the time-averaged form of the Navier-Stokes equations. Many turbulence models were developed based on turbulent viscosity relation by means of the Boussinesq hypothesis, commonly referred to as eddy viscosity turbulence models (EVTMs). Examples of the EVTMs are the Spalart-Allmaras (SA) based model (Equation (1)), k − based model (Equation (2)), k − model (Equation (2)) as well as the Reynolds stress based model (Equation (7)). Among these models, the Menter shear-stress-transport (SST) k − turbulence model [30] is widely used for many engineering problems. This model has been proven to be accurate for flows involving massive separation. Therefore, this model will be used for the present investigation. This is in line with past studies concerning the choice of turbulence models for VAWT simulations. [25,[31][32][33][34][35][36] McLaren et al. [35,36] demonstrated that the SST model was superior in predicting the deep stall characteristics of airfoil, being one of the main importance for high solidity VAWT. They also indicated that the SST model was able to reconstruct the power curve as compared with the measurement data. A similar conclusion was obtained from Danao et al. [34] for the SST based turbulence models. Bianchini et al. [25] provided a very interesting study concerning the feasibility of the SST based turbulence models, with and without considering transitions. They demonstrated that both models provide good consistency for a wide range of flow cases. He et al. [32] also reached a similar conclusion as compared with the 3D large eddy simulation (LES) and experimental data. Having reviewed these available research studies, the usage of the SST k − model in the present paper is in line with past works in VAWT aerodynamics. The FLOWer code applied in the present study solves the compressible RANS equations applying a finite volume formulation on a block-structured mesh, making it possible to simulate complex shapes and flows. [37] The FLOWer code was continuously extended by the Institute of Aerodynamics and Gas Dynamics (IAG) at the University of Stuttgart for the simulation of wind turbines. The Chimera approach is used in the study to simplify the generation of high quality mesh components of complex structures. The code employs a central cell-centered Jameson-Schmidt-Turkel (JST) [38] method for flux computation resulting in a second order spatial accuracy on smooth meshes. The method makes use of a central space discretization with artificial dissipation being calculated in relation to the grid cell aspect ratio. [37] This approach is robust and well suited for parallel application. [37] Furthermore, the time integration for the unsteady calculations is conducted under the use of an explicit 5stage Runge-Kutta scheme, which can be accelerated by methods of local time stepping, resulting in a second order accuracy in time. The multigrid level 3 is further used to accelerate convergence. The simulations were carried out assuming a fully turbulent boundary layer condition. The physical time step size used is equivalent to 1 • of blade revolution with 35 sub-iterations per time step. The numerical setup is summarized in Table 2. The FLOWer code has been validated at IAG for VAWT computations under the turbulent wake state and under strong dynamic stall conditions. [10,11,39,40] In general, a good simulation model is required to obtain reliable numerical solutions. Since the structures of a VAWT are not so simple, it will not be convenient to model the VAWT using one single grid. In the present study, the structures of the VAWT will be reconstructed into several simple structures based on its purpose, such as the background mesh, the refined area and the blade (airfoil) itself. According to preceding studies carried out by Castelli et al., [41] Song et al., [42] and Bangga et al., [11] some detailed rotor components like bolts and member bars that connect blades www.advancedsciencenews.com www.advtheorysimul.com with the main shaft can be omitted in the simulations without sacrificing the simulation accuracy. 2D geometry was applied in the study as the blade section is exactly the same along the blade length, neglecting the influence of the tip vortices.
The FLOWer code supports the usage of the overlapping mesh (chimera) technique. This allows several grid components to be built independently, making the grid generation process easier without sacrificing the accuracy and quality of the mesh. The computational grid consists of five main components, namely coarse background (black color), near wake refinement (red color), far wake refinement (green color), ring refinement (blue), and rotor and shaft (purple color) meshes as shown in Figure 3. The background domain is large enough to contain the wake in all directions. A very recent study [39] has shown that this plays an important role for accuracy if the rotor loading is large.
The near and far wake refinement meshes have an equidistance cell size. The near wake refinement has a resolution of Δ∕R = 0.025, where R is the turbine radius (m), which is far finer than previous studies in. [11] The far wake refinement has a resolution of Δ∕R = 0.5. The ring refinement domain has a circular shape of a radius of 1.7R, with a resolution of around 3% of the chord. Based on a grid study, the blade mesh is composed of 256 cells in circumferential direction. The maximum cell size in the rotor region is less than ≈3% of the chord. Furthermore, the condition y + < 1 is maintained with 64 layers in the boundary layer region. The Chimera intersection area was defined near the outer area of the rotating zone for at least four cells overlap. The communication between the individual grids is established by interpolating data from the overlapping grids. The total number of cells in all zones is about 776 704.
Grid study has been performed by increasing the cell number in the rotor and in the wake area by uniformly varying the cells number. This was done in two steps; 1) first by varying the number of grid cells on a static airfoil case to evaluate the lift characteristics ( Table 3); 2) the suitable airfoil grid is then applied for the rotor computations and the background mesh  was refined subsequently ( Table 4). The grid convergence index (GCI) [43] is quantified to show the basis of the grid selection.
The error between each investigated grid with the "extrapolated" grid independent value [43] is also shown. One can see that all investigated airfoil meshes have errors well below 5%. This indicates that any grid will be suitable to be used. Therefore, G1 is chosen for further investigation to limit the number of grid cells. A slightly different result is obtained in Table 4. It can be seen that refining the wake mesh improves the prediction considerably from 7.59% to around 2.04% for the finest mesh. Therefore, to ensure that the CFD solutions are reliable, GW3 will be used throughout the paper.

Validation of the Numerical Results
The simulations were performed for ten rotor revolutions until the wake fully develops. Some authors [20,28,44] recommended at least 20-30 revolutions for reaching the convergence. These studies, however, were done for low to medium solidity turbines operating at relatively large tip speed ratio (e.g., 2.5, 3.5, and 4.5 in ref. [44]). To the authors' experience, for turbines at low tip speed ratios, this number can be reduced as shown in our previous papers. [11,40] Therefore, a simple estimation is made in this paper to evaluate why such differences occur. First, consider the tip speed ratio as = (ΩR)∕U ∞ , where Ω is the turbine rotational speed and U ∞ , freestream wind speed. From this equation one can calculate the rotation period of the turbine as T = (2 R)∕( U ∞ ). Therefore, the distance traveled by the fluid flow past through the turbine may be approximated as: with N Rot being the number of rotor rotations. This relation is plotted in Figure 4. It can be clearly seen that the distance traveled by the fluid flow in axial direction scales greatly with . For the same number of rotation, x traveled for the smaller value is far greater than than the larger one. As an example, for = 1.2 and N Rot = 10, x traveled has reached as far as 62.83 relative to the rotor radius, while for = 4.5 and N Rot = 20 it is only as low as x traveled = 27.92R. This explains the deviations often found in several simulation studies. [11,20,28,40,44] Figure 5 presents the predicted driving moment coefficient for the last four rotor revolutions for blade-1. One can see that a perfect periodicity cannot be achieved because the rotor operates at a small tip speed ratio, where the dynamic stall and flow separation effects are prominent. This is especially true for the thicker airfoil where, by its nature, is prone to exhibit flow separation during its operation. This effect is especially dominant in the downwind phase because the separated flow from the blade interacts with the wake emanating from the other blade and from its own during the upwind phase. The integrated results are then obtained by averaging the last three revolutions. The power coefficient can be calculated as: and where Ω is the rotational frequency of the rotor, C Moment is the driving moment coefficient, and A defines the cross section area of 2R ⋅ H. For two dimensional case H equals to unity.
To assess the accuracy of the present CFD computation, the experimental and LES results from ref. [29] are employed for comparison. The experiments [29] were carried out in a large scale open circuit wind tunnel with the test section of 1.2 m × 1.2 m and 1.4 m long. The turbulence intensity level at the inlet was less than 0.8%. The LES study carried out in ref. [29] was performed using the commercial code ANSYS Fluent with a similar setup as in the experiment. The results of the present study are displayed in Figure 6 in terms of the power coefficient for several tip speed ratios. It can be clearly seen that the present simulation is in a sound agreement with the reference data. Since the averaged power is obtained from the last three rotor rotations, there is a chance that the generality could be affected by the evaluation period. Therefore, this aspect is further investigated in Figure 7. Here one can see that using the last three rotor rotations is reasonable since the obtained power coefficient does not change significantly. The power coefficient still reduces slightly for the NACA 0040 due to the effects of flow separation in the downwind region (influencing the periodicity). However, this reduction of the power level is still at a reasonable level and this should not affect the generality of the obtained conclusions.

Load Characteristics
This section presents the effects of airfoil thickness on the load characteristics. The global integrated loads are presented in Figure 8 in terms of power and thrust coefficients, with the later defined as: where C Thrust is the thrust coefficient and F Thrust is the thrust force. It can be seen that C Power increases significantly with increasing airfoil thickness from t∕c = 8%, 12% to 18%. Then the power coefficient remains at a similar level until it slowly decreases again for t∕c = 40%. The thrust coefficient, on the other hand, generally shows a steady reduction with increasing airfoil thickness. One can see that the power and thrust generated by the turbine are not always in line with the best-fitted trend. The observed deviations are caused by the dynamic stall effects. It will  be shown in the subsequent discussions for the azimuthal loads, flow field, and pressure distributions, that dynamic stall could yield an increased normal force coefficient without an increase in power due to a significant drag augmentation. This occurs especially for thin airfoils having strong suction peak near the leading edge area. This is the reason of the power reduction for the NACA 0008 and 0012 airfoils. However, by having the right amount of thickness, the increase of lift (i.e., normal force component) due to the dynamic stall effects is not followed by the massive increase in drag. This aspect becomes the main source of the additional power increase for the NACA 0018 and 0021 airfoils. Further increase of the thickness yields the reduction of the dynamic stall effects and thus the performance of the rotor is only determined by the aerodynamic characteristics of the airfoil itself.
To gain deeper insights about the reason for such behavior, the azimuthal loads of blade-1 are presented in Figure 9, with C Moment being the driving moment, C T being the tangential force defined in Figure 2 and C M being the pitching moment relative to the blade-mounting location. These forces are defined as: where C T , tangential force coefficient; F T , tangential force; C N , normal force coefficient; F N , normal force; C M , pitching moment coefficient; and F M , pitching moment. One can see that the maximum of driving moment coefficient (C Moment ) increases with airfoil thickness up to t∕c = 21%, and afterwards gradually decreases. Interestingly, the location of the maximum point also shifts from blade-1 azimuth angle, ≈ 50 • (for NACA 0008), to ≈ 70 • (for NACA 0012) up to ≈ 115 • (for NACA 0018, 0021, and 0025). Only by increasing the thickness further by 40%, the location of the maximum value shifts back to a smaller azimuth angle at ≈ 95 • . This observation is also supported by the tangential force distribution in Figure 9. The location of the maximum driving moment (M) (and tangential force) indicates the location of maximum airfoil efficiency for a certain rotor configuration and aerodynamic polar, that is, the maximum ratio of the lift to drag force. For the azimuth angle larger than this position, the performance of the blade gradually decreases. In the downwind region within 180 • < < 360 • , the blade generally no longer contributes to the rotor power, or even creates negative power output. It can be observed that the NACA 0008 and NACA 0012 airfoils have a much smaller azimuth regime with positive C Moment than the other airfoils. This aspect causes a great increase of the power coefficient by a slight increase in the relative thickness with respect to the NACA 0018 airfoil in Figure 8.
In contrast with the driving moment and the tangential force coefficients, the maximum peaks of the normal force and pitching moment coefficients for all airfoils are located at an approximately the same position at ≈ 70 • and ≈ 90 • , respectively. Note that the normal force is proportional with the lift force. Similarly, the pitching moment is mainly driven by the lift force. Therefore the maximum normal force and pitching moment indicate the peaks of the dynamic stall cycle over one rotor revolution. This force component is also the main driver for the thrust force presented in Figure 8. One can see clearly that the normal force consistently decreases with increasing airfoil thickness. This indicates that the vortex lift effect of the dynamic stall cycle reduces with increasing airfoil thickness. As a result, this observation explains a consistent reduction of the thrust force with t∕c in Figure 8. Similarly, as the vortex lift contribution reduces with t∕c, the pitching moment also becomes more positive with airfoil thickness.
Unlike all the other airfoils, NACA 0040 experiences a spike in C Moment and C T in the downwind regime at ≈ 230 • , which may indicate a different flow regime for this airfoil. To obtain a better overview concerning this aspect, the streamlines surrounding blade-1 are plotted in Figure 10. It can be seen clearly that the flow fields for three different airfoils are different (NACA 0008, NACA 0021, and NACA 0040). These are particularly true concerning the incoming flow angle and on the vortical structure attached on the inner side of the blade. The inflow angle changes from positive to negative with increasing airfoil thickness, while it remains around zero for the NACA 0021 airfoil. This directly affects the force direction acting on the blade. Furthermore, the size and strength of the vortical structure near the trailing edge are also decreasing with airfoil thickness.

Flow Field Development
In this section, the development of the flow field surrounding the rotating blade will be investigated in detail for several important azimuth angles. The analyses are presented in Figure 11  For the smallest investigated azimuth angle in Figure 11, one can see that the flow is attached on the airfoil surface for thin to medium airfoil thickness. For the NACA 0040 airfoil, light flow separation is observed near the trailing edge. This creates a periodic wake of the von Kármán type to travel downstream. Interestingly, for this extremely high solidity turbine, the vorticity field is strongly deflected clockwise right after the flow passes the trailing edge. This causes a strong flow curvature effect. This highlights the importance of modeling the decambering effect in engineering models for high solidity turbine because the blade section sees a large variation of the angle of attack along the chord, that is, as documented in refs. [39,40]. A similar characteristic is also observed in Figure 12.
Further increasing the azimuth angle to = 80 • , a strong leading edge separation bubble is clearly observed for the thinnest investigated airfoil, NACA 0008. This bubble creates a strong suction effect causing a large increase of the normal force component, while at the same time also increases the drag force significantly. As a result, the normal force is significantly large in comparison to the other airfoils, whilst the generated power is the smallest. On the other hand, trailing edge separation increases in intensity with increasing airfoil thickness. Further increasing the azimuth angle to 108 • , one can see that a strong leading edge separation bubble also occurs on the inner side of    the NACA 0012 airfoil. Having a look back at the azimuthal load distributions in Figure 9, one can see that the driving moment of the NACA 0012 airfoil also already drops at this location. In contrast, the normal force still remains high for this airfoil. This confirms that the leading edge bubble of the dynamic stall cycle is the main cause of the power loss for a vertical axis wind turbine under the influence of dynamic stall.
For the higher azimuth angle at = 140 • , now the NACA 0018 airfoil is also characterized by the leading edge vortex bubble in combination with trailing edge separation. At this position, the driving moment already drops for all airfoils. For the other thicker airfoils, the leading edge bubble is also produced, but it is suppressed significantly with increasing airfoil thickness. Finally, all airfoils are fully separated as the blade reaches the azimuth angle of 180 • . However, one can see clearly that the strength of the detached vortex near the airfoil surface is visually weaker for the thicker airfoil.
To better quantify the vortex characteristics, a deeper analysis is done at the azimuth angle of 180 • . The dynamic stall vortex is evaluated in terms of the vortex core radius (r core ), absolute maximum vorticity at the vortex center ( Z , core), and integrated circulation (Γ core ). The latter is calculated by: with A Z being a circular area enclosing the vortex core according to Holzäpfel. [45] Furthermore, the center of the vortex is defined as the location with the maximum absolute vorticity according to Vollmers. [46] By this definition, then one is able to calculate the vortex core diameter by quantifying the distance between the minimum and maximum x−velocity components. This procedure has been implemented and tested for a horizontal axis wind turbine application in ref. [47]. The results are plotted in Figure 17.
It can be observed in Figure 17 that the vortex core radius gradually reduces with increasing airfoil thickness up to t∕c = 25%. Further increasing the relative thickness yields a fairly constant core radius. In contrast, the maximum vorticity at the vortex core increases with airfoil thickness up to t∕c = 21% before it constantly decreases. The increased maximum core vorticity with thickness for this specific regime is strongly related with the age of the vortex. From the above discussion, one realizes that the thinner airfoil generates the leading edge vortex at a much earlier azimuth angle. Therefore, the energy of the vortex core is already dissipated into the surrounding area. This is confirmed by the integrated circulation in Figure 17 (where it also covers the surrounding regime, not only the vortex center). Now one can see clearly that the vortex strength is almost constant up to t∕c = 21%. Further increase in the airfoil thickness causing a sharp change of the curve to be observed. The reason lies in the fact that the leading edge vortex is significantly suppressed by the airfoil thickness as becoming evident from Figure 16.

Pressure Coefficient Development
In the present section, evaluations of the pressure coefficient distribution over the azimuth angle for several investigated airfoils are conducted, defined as: where p, static pressure; p ∞ , freestream pressure; and , air density. Figure 18 presents the pressure distribution on the outer side of the blade for blade-1. This part becomes the pressure side of the blade for 0 • < < 180 • and switches to the suction side for the rest of the azimuth angles. Opposite is true for Figure 19 which shows the pressure distribution on the inner part of the blade. At = 0 • , the pressure level shows some variation on the outer side of the blade, but less on the inner part except for the thickest airfoil (40%) which shows a strong flow acceleration. Referring back to Figure 9, one can see that most airfoils generate almost the same level of the driving moment except for the NACA 0040 airfoil. This explains the source of the difference. Interestingly, usually the suction side is more sensitive to changes than the pressure side. Opposite is true in the present observation at small angle of attack (note that for this azimuth angle, the geometric angle of attack should be close to zero). This flow characteristic can be attributed to the effects of flow curvature for high solidity turbines. Recent studies highlighted the importance of considering this aspect on engineering models. [39,40] In these studies, the correlation was modeled using the chord-to-radius ratio, while the airfoil thickness was neglected. The present study highlights that this might be insufficient providing the relative thickness is large enough.
By increasing the azimuth angle to 40 • , the suction side shows less variation than the pressure side except for the NACA 0040 airfoil for x∕c > 0.1. Despite that, a very strong suction peak is observed for the thinnest airfoil (NACA 0008), which indicates a very strong adverse pressure gradient. Further increase of causes this strongly accelerated flow to detach, creating an intense leading edge vortex bubble. As observed at = 80 • , the peak shifts to x∕c ≈ 0.15, while the pressure level increases at the vortex location (becomes more positive). Starting from this azimuth location until the blade returns to the zero angle of attack at = 180 • , the pressure variation on the outer side of the blade is no longer significant. In contrast, the pressure coefficient strongly varies on the inner part of the blade, which dominates the load fluctuations.
Having a look at Figure 19, one notices a certain pattern governs the dynamic stall process. First, the suction peak appears causing a strong adverse pressure gradient just in the vicinity of the leading edge of the airfoil. Then, flow separation takes place at this location abruptly that is indicated by the increase of the pressure level. Despite that, the dynamic motion of the airfoil causes the airfoil to see "delayed-response" of the incoming flow, see for example explanation of dynamic stall in ref. [48]. The delayed response causes the flow just downstream of the separation area near the leading edge to remain attached on the airfoil surface, creating the separation bubble widely known as the leading edge vortex. The more negative the suction peak, the stronger the effects become. As one can see, the suction peak clearly weakens with increasing airfoil thickness. Furthermore, because the leading edge radius is very large for the thickest airfoil, leading edge separation is not observed at all. In fact, for this thick airfoil, the reduction of the turbine performance is mainly driven by trailing edge separation.

Conclusions
The effects of airfoil thickness on the dynamic stall characteristics of high solidity vertical axis wind turbines have been thoroughly investigated in the present paper. The studies were conducted employing a computational fluid dynamic approach. Seven different airfoils with varying thickness from the NACA airfoil family have been evaluated systematically. Several conclusions can be drawn from the paper: • The turbine power coefficient increases with airfoil thickness significantly until t∕c ≈ 25 − 30%. Then it reduces gradually due to trailing edge separation effects.
• A gradual reduction of thrust is observed with increasing airfoil thickness up to t∕c ≈ 40%. The reason lies in the reduction of the suction peak effects. • A large normal force coefficient during dynamic stall does not always positively contribute to large power production because drag also increases considerably, causing a significant reduction in the airfoil efficiency. As a result, the tangential force drops while the normal force still increases. • The leading edge vortex bubble is the main cause of the power loss for high solidity vertical axis wind turbines. • The effects of dynamic stall decrease with increasing airfoil thickness. • The leading edge vortex strength and radius reduce with increasing airfoil thickness. • At small azimuth angles, the airfoil thickness has a stronger influence on the pressure level on the outer side of the blade than the inner side. The effects can be attributed to the flow curvature effects. • At larger azimuth angles, as flow separation already takes place, the inner side of the blade plays a further significant role. In general, the dynamic stall process is initiated by the occurrence of the suction peak associated with the strong adverse pressure gradient. The effect becomes weaker with increasing airfoil thickness.
From the findings obtained in the present work, it can be concluded that the airfoil thickness plays an important role on the rotor performance. Using thin airfoils for high solidity rotor is not desirable because the dynamic stall effects cause massive reduction in rotor performance. On the other hand, thick airfoils suffer from trailing edge separation even at small angles of attack. Therefore, it is recommended to utilize medium airfoil thickness ranging from 18% to 30% relative to the chord length for obtaining the optimal rotor performance.