A Large-Eddy Simulation Study of Vertical Axis Wind Turbine Wakes in the Atmospheric Boundary Layer

Abstract: In a future sustainable energy vision, in which diversified conversion of renewable energies is essential, vertical axis wind turbines (VAWTs) exhibit some potential as a reliable means of wind energy extraction alongside conventional horizontal axis wind turbines (HAWTs). Nevertheless, there is currently a relative shortage of scientific, academic and technical investigations of VAWTs as compared to HAWTs. Having this in mind, in this work, we aim to, for the first time, study the wake of a single VAWT placed in the atmospheric boundary layer using large-eddy simulation (LES). To do this, we use a previously-validated LES framework in which an actuator line model (ALM) is incorporated. First, for a typical threeand straight-bladed 1-MW VAWT design, the variation of the power coefficient with both the chord length of the blades and the tip-speed ratio is analyzed by performing 117 simulations using LES-ALM. The optimum combination of solidity (defined as Nc/R, where N is the number of blades, c is the chord length and R is the rotor radius) and tip-speed ratio is found to be 0.18 and 4.5, respectively. Subsequently, the wake of a VAWT with these optimum specifications is thoroughly examined by showing different relevant mean and turbulence wake flow statistics. It is found that for this case, the maximum velocity deficit at the equator height of the turbine occurs 2.7 rotor diameters downstream of the center of the turbine, and only after that point, the wake starts to recover. Moreover, it is observed that the maximum turbulence intensity (TI) at the equator height of the turbine occurs at a distance of about 3.8 rotor diameters downstream of the turbine. As we move towards the upper and lower edges of the turbine, the maximum TI (at a certain height) increases, and its location moves relatively closer to the turbine. Furthermore, whereas both TI and turbulent momentum flux fields show clear vertical asymmetries (with larger magnitudes at the upper wake edge compared to the ones at the lower edge), only slight lateral asymmetries were observed at the optimum tip-speed ratio for which the simulations were performed.


Introduction
Vertical axis wind turbines (VAWTs) offer some advantages over their horizontal axis counterparts and are being considered as a viable alternative to horizontal axis wind turbines (HAWTs). The research on VAWT technology started in the 1970s, and while the main focus of the performed studies has been on the overall turbine performance (quantities such as power and torque) and on the mechanical loading on the blades, relatively few studies have attempted to analyze the wake of a VAWT (for a comprehensive and chronological review of the studies on VAWTs before 2000, see Paraschivoiu [1] (Chapters 4-7)). Having a thorough understanding of VAWT wakes is especially crucial in designing VAWT wind farms, where downstream turbines can potentially be located in the wake of upstream ones, and consequently, the performance of the whole wind farm could be significantly affected by the wake flow characteristics. Among the experimental works investigating VAWT wakes, one can find a relatively larger number of studies that have focused only on the near wake region (e.g., [2][3][4][5]), compared to those that have considered also the far wake region (e.g., [6][7][8]). Nevertheless, from a wind farm design point of view, it is the far wake behavior of the flow that has more relevance and importance.
In the numerical flow simulation domain, the studies performed on the flow through VAWTs can be divided into two main categories: (1) the simulations in which the blades of the turbine (and consequently, the boundary layer around them) are resolved; and (2) the simulations in which the blades are modeled by an actuator-type technique, which uses immersed-body forces to take into account the effects of the blades on the flow. While the first approach (for instance, the work of Castelli et al. [9]) can be highly valuable to calculate the loading on the blades and the flow characteristics inside the rotor and in the near wake, to simulate the far wake of VAWTs and especially VAWT wind farms, the second approach is deemed to be more feasible and attainable [10,11]. The use of actuator-type techniques for VAWTs dates back to the 1980s, when Rajagopalan and Fanucci [12] for the first time modeled the VAWT rotor by a porous surface, swept by the blades, on which time-averaged blade forces are distributed and continuously act on the flow (which has also been called the actuator swept-surface model [11]). An extension of this work to three dimensions was made by Rajagopalan et al. [13]. Later on, Shen et al. [14] introduced the actuator surface model and employed it to obtain the flow field past a VAWT in two dimensions. More recently, Shamsoddin and Porté-Agel [11] used large-eddy simulation (LES) coupled with both the actuator-swept surface model (ASSM) and the actuator line model (ALM) to simulate the flow through a VAWT placed in a water channel and compared the resulting wake profiles with experimental data.
Acknowledging the fact that any given real VAWT is likely to be working in the atmospheric boundary layer (ABL) and benefiting from the helpful experience gained from the extensive research on HAWT wakes, it is imperative to study in detail the characteristics of the wake of VAWTs placed in boundary layer flows, especially if VAWT farms are to be envisaged as a viable source of power in future energy outlooks. Having this in mind, the present study is a step in this direction and attempts to use a previously-validated LES framework, in which an actuator line model is incorporated, to analyze the wake of a typical straight-bladed VAWT in a relatively long downstream range. Moreover, before the wake study, using the same framework, the power production performance of the VAWT for different combinations of blade chord lengths and tip-speed ratios is studied to find the optimum combination for the aforementioned wake analysis. To the best knowledge of the authors, this study is the first attempt to characterize the wake of a VAWT in ABL using LES.
The LES framework is presented in Section 2, and the numerical setups and techniques are described in Section 3. Next, the results for both the power production parametric study and the wake analysis are presented and discussed in Section 4. Finally, a summary of the study is given in Section 5.

Large-Eddy Simulation Framework
In the LES framework used for the simulations of this paper, the filtered incompressible Navier-Stokes equations (for a neutrally-stratified ABL) are solved. These equations can be written in rotational form as: where the tilde represents a three-dimensional spatial filtering operation at scale ∆, u i is the filtered velocity in the i-th direction (with i = 1, 2, 3 corresponding to the streamwise (x), spanwise (y) and vertical (z) directions, respectively),p * = p ρ + 1 2 u i u i is the modified kinematic pressure where p is the filtered pressure, τ ij = u i u j − u i u j is the kinematic subgrid-scale (SGS) stress, f i is a body force (per unit volume) representing the force exerted by the flow on the turbine blades (observe the minus sign), F p is an imposed pressure gradient and ρ is the constant fluid density. In this paper, u, v and w notations are also used for the u 1 , u 2 and u 3 velocity components, respectively. Regarding the parametrization of the SGS stresses, in these simulations, the Lagrangian scale-dependent dynamic model [15] is used.
To parameterize the VAWT-induced forces on the flow (i.e., to model the term f i /ρ in Equation (2)), an actuator line model is used. According to the ALM, each blade of the turbine is represented by an actuator line on which the turbine forces, calculated based on the blade-element theory, are distributed. This method has the advantage of being capable of tracking the rotation of the blades at each time step. For a detailed explanation of the application of the ALM for VAWTs, the reader can refer to Shamsoddin and Porté-Agel [11] (Section 2.2).

Numerical Setup
In this section, the techniques used to numerically solve Equations (1) and (2), as well as the configuration of the performed numerical experiments are presented.
The LES code, which is used to realize the simulations in this study, is a modified version of the code described by Albertson and Parlange [16], Porté-Agel et al. [17] and Porté-Agel et al. [18]. The computational mesh is a 3D structured one, which has N x , N y and N z nodes in the x, y and z directions, respectively. The mesh is staggered in the z direction in a way that the layers in which the vertical component of velocity (w) is stored are located halfway between the layers in which all of the other main flow variables (u, v, p) are stored. The first w-nodes are located on the z = 0 plane, while the first uvp-nodes are located on the z = ∆z/2 plane.
To compute the spatial derivatives, a Fourier-based pseudospectral scheme is used in the horizontal directions, and a second-order finite difference method is used in the vertical direction. The governing equations for conservation of momentum are integrated in time with the second-order Adams-Bashforth scheme.
The pressure term in Equation (2) is not a thermodynamic quantity, and it only serves to have a divergence-free (i.e., incompressible) velocity field. Therefore, by taking the divergence of the momentum Equation (2) and using the continuity Equation (1), we can solve the arising Poisson equation for the modified pressure,p * , using the spectral method in the horizontal directions and finite differences in the vertical direction.
The boundary conditions (BCs) in the horizontal directions are mathematically (and implicitly through using the spectral method) periodic. For the bottom BC, the instantaneous surface shear stress is calculated using the Monin-Obukhov similarity theory [19] as a function of the local horizontal velocities at the nearest (to the surface) vertical grid points (z = ∆z/2) (see, for instance, Moeng [20], Stoll and Porté-Agel [21]). For the upper boundary, an impermeable stress-free BC is applied, i.e., ∂ u 1 /∂z = ∂ u 2 /∂z = u 3 = 0.
Since the study of the flow through a single turbine is desired, we need to numerically enforce an inflow BC to practically override the implicitly-imposed periodic BC in the x direction. For this purpose, a buffer zone upstream of the VAWT is employed to adjust the flow to an undisturbed ABL inflow condition. The inflow field is obtained by saving the instantaneous velocity components in a specific y-z plane in a similar precursory simulation of ABL over a flat terrain (with the same surface roughness) with no turbine on it. The use of this technique, i.e., using an inflow boundary condition in a direction in which the flow variables are discretized using Fourier series, has been shown to be successful in the works of Tseng et al. [22], Wu and Porté-Agel [23] and Porté-Agel et al. [24].
To implement the ALM, values of the airfoil's lift and drag coefficients (C L and C D , respectively) as a function of Reynolds number (Re) and angle of attack (α) (i.e., C L(D) = f (Re, α)) are needed. This information was obtained from the tabulated data provided by Sheldahl and Klimas [25]. Moreover, the dynamic stall phenomenon, which is known to have a considerable effect on the performance of VAWTs [14,26], is accounted for using the modified MIT model [27]. A detailed explanation of the implementation of the dynamic stall model is provided in Appendix A. Figures 1 and 2 show the geometrical specifications of the VAWT and the computational domain in which it is placed. The turbine rotor is made of three straight blades and has a diameter (D) of 50 m and a height of 100 m. The blades' airfoil is selected to be the symmetrical NACA 0018 airfoil, which is widely used for VAWTs. It is attempted that these chosen turbine specifications are representative of those of real VAWTs with a nominal capacity of 1 MW (this fact will be reaffirmed by the results of the simulations). For example, a curve-bladed (or Φ-rotor) VAWT of similar size and capacity (96 m high and an equatorial diameter of 64 m) with two NACA 0018 blades of a 2.4-m chord length was operational as part of Project Éole in Cap Chat, Quebec, Canada, between 1987 and 1993 [28]. This turbine was designed to deliver a maximum power of about 4 MW (at high winds and high rotational speeds), and its maximum measured power of about 1.3 MW is hitherto one of the greatest measured power outputs for a VAWT ([1] Section 7.3.4). The buffer zone occupies about 12% of the domain length. The domain dimensions are L x = 1200 m (=24D), L y = 600 m (=12D) and L z = 400 m (=8D) in the streamwise, spanwise and vertical directions, respectively. The blockage ratio of the turbine in the computational domain is 2.08%, which is well below the value of 10%, which is reported by Chen and Liou [29] as the threshold below which it is acceptable to neglect the blockage effect. Regarding the computational mesh, the number of grid points in each of the three directions is N x = 360, N y = 180 and N z = 240. The code has been shown to yield grid-independent results provided that a minimum number of grid points is used to resolve the rotor [11]. In this study, we chose a resolution (15 points in each horizontal direction covering the rotor area) that falls within the grid-independent range. The time resolution for all of the simulations is 0.0155 s. For the wake study simulation, the total physical time of the simulation is 90.4 min, and for mean velocity and turbulence statistics results, we have time-averaged the quantities in question over the final 77.5-min time span.  Figure 3 shows mean and standard deviation profiles of the inflow streamwise velocity. As mentioned earlier, the inflow field is generated by using the flow field of a precursory simulation of the neutrally-stratified ABL on a flat terrain. The surface roughness, z o , and the friction velocity, u * , used in this precursory simulation are 0.1 m and 0.52 m/s, respectively. In Figure 3a, it can be seen that the mean streamwise velocity profile approximately follows the log law in the surface layer. The mean inflow streamwise velocity at the equator height of the turbine (i.e., z = 100 m in this case), U eq , and the turbulence intensity of the inflow at the same height (σ u /U eq ) are 9.6 m/s and 8.3%, respectively. It should be noted that the above-mentioned inflow field is used for all of the simulations of this paper (i.e., both Subsections 4.1 and 4.2).

Results and Discussion
In this section, the results of the simulations are presented and discussed. First, we examine the turbine's energy-extraction performance, and next, we study the wake flow of a VAWT placed in the ABL.

Turbine Performance and Power Extraction
In this subsection, we are interested in how the power production of the turbine is affected by different combinations of tip-speed ratio, TSR, and chord length, c. For this purpose, 117 simulations have been performed to obtain the power coefficient, C P , of the turbine as a function of both TSR and c, i.e., C P (TSR, c). Figure 4 shows how C P varies with different values of TSR and c. Figure 4a is generated with a resolution of 0.5 m for chord length and 0.5 for TSR. It can be seen that, as we increase the chord length, the useful TSR range (a range in which C P is higher than a certain value) decreases. Moreover, the figure shows that the maximum power coefficient of the turbine occurs for a TSR of 4.5 and a chord length of 1.5 m (which corresponds to a solidity of Nc/R = 0.18, where N is the number of blades and R is the rotor radius). This combination results in a power extraction, P, of 1.3 MW and a C P of 0.47 (C P is defined as C P = P/(0.5ρDHU eq 3 ), where ρ is the fluid density and considered equal to 1.225 kg/m 3 , D is the rotor diameter and H is the rotor height).

VAWT Wake
In this subsection, we have picked the optimum combination of TSR and chord length (TSR = 4.5 and c = 1.5 m) for the VAWT rotor and studied the wake flow behind it. Figure 5 shows the instantaneous streamwise velocity field of the flow in three different orthogonal planes. In all of the following figures in this section containing contour plots, the black circles and rectangles represent the outline of the locus of the blades. The sense of the rotation of the turbine blades is counterclockwise when seen from above. The wake of the VAWT and the highly turbulent nature of the flow are obvious in this figure and in the Videos S1 and S2 included in the Supplementary Material. It should be noted that the average thrust coefficient of the turbine (defined as C T = T/(0.5ρDHU eq 2 ), where T is the total thrust force of the turbine in the x direction) in this case is found to be 0.8.  7 show contour plots of the mean streamwise velocity in the x-y plane at the equator height of the turbine and in the x-z mid-plane of the turbine. It can be seen in these figures that it takes a long distance for the wake to recover; at a downwind distance as large as 14 rotor diameters, the wake center velocity reaches only 85% of the incoming velocity. Moreover, Figure 8 shows the mean velocity contours in six y-z planes downstream of the turbine. In all of these figures (Figures 6-8), one can observe how the wake recovers in the streamwise direction (after a certain distance) and how it expands in the spanwise direction as it advances farther downstream. To have a more quantitative and precise insight about the VAWT wake, Figures 9 and 10 can be consulted. Figure 9 shows spanwise profiles of the mean streamwise velocity in a horizontal plane at the equator height of the turbine in eight downstream positions. Besides, Figure 10 presents vertical profiles of the mean streamwise velocity in the x-z mid-plane of the turbine at different downstream positions. An interesting observation that can be made from Figures 6, 9 and 10 is that the maximum velocity deficit occurs at a downstream distance of about 2.7 rotor diameters; this distance is significantly larger than the equivalent one for the case of HAWT wakes [23]. After the point where the maximum velocity deficit (more than 65% of U eq in this case) occurs has been reached, the wake starts to recover with a relatively high recovery rate (defined here as the magnitude of the rate of change of the maximum velocity deficit with streamwise distance). As we go farther downstream, the rate of the wake recovery decreases considerably; so that in the distances as large as 17 rotor diameters, where the velocity deficit reaches values of about 90%, the recovery rate is comparably very small. Another group of crucial quantities that has a significant importance in characterizing turbine wakes is the turbulence-related statistics, such as turbulence intensity and turbulent fluxes. These quantities are especially important for the design of wind farms, due to their role in both wake recovery and mechanical loads on turbine blades. Figure 11 shows contours of turbulence intensity (TI) in two different orthogonal planes (x-y and x-z) in the wake of the turbine. Here, the turbulence intensity is defined as TI = σ u /U eq . In addition, Figure 12 shows the distribution of TI in y-z planes at different downstream locations. In Figure 11a, it can be seen that two branches of high TI regions start to develop from the two spanwise extremities of the rotor swept surface (the black circle in the figure). These two branches grow in spanwise width as we go further downstream, until the point where they meet each other (for this case, in about 3.5 rotor diameters downstream of the turbine in the horizontal mid-plane of the turbine). Starting from the turbine area, the TI in each of these branches increases, until a point where the maximum TI occurs (about 3.8 rotor diameters downstream in this case); after this maximum point, the TI starts to decrease as the flow advances downstream, while the width of the branches continues to expand. Figure 13 examines the previous figure quantitatively, by showing the spanwise profiles of the TI in the equator height of the turbine. One can readily see that at each downstream position, the horizontal TI profiles have two maxima at two spanwise positions, which correspond to the two aforementioned TI branches. Although slight asymmetries can still be seen in the TI values of the two branches, the lateral asymmetry is not significantly pronounced. It should be noted that the degree to which the VAWT wake is laterally asymmetric is influenced by parameters, such as TSR, airfoil type and the Reynolds number in which the turbine is working.   In Figure 11b, one can observe a similar behavior by noticing the two high TI regions originating from the upper and lower extremities of the turbine; however, in this case, the TI originating from the upper edge of the blades is clearly larger than the one originating from the lower edge. To further quantify this, and to have a better understanding of the vertical variation of the TI in a VAWT wake, one can study Figure 14, in which vertical profiles of TI are shown at different downstream positions. In this figure, it can also be seen that in the region below the turbine blades' lower edge, the turbulence intensity has even decreased to values lower than the inflow TI; this behavior has also been observed in HAWT wakes, as well (e.g., [23]).
Furthermore, turbulent momentum fluxes in the VAWT wake are believed to be worthy of inspection, as they quantify the rate of flow entrainment into the wake, which is responsible for the recovery and lateral expansion of the wake. Figure 15 shows the normalized lateral turbulent flux (u v ) at the horizontal mid-plane of the turbine. The positive and negative regions of u v , which are located on the two lateral edges of the wake, show an inward entrainment of momentum into the wake region. This lateral entrainment can also be seen in the y-z planes in Figure 16. Figure 17 shows the spanwise profiles of u v at the equator height of the turbine at different downstream distances from the turbine. We notice that for this case, the maximum absolute value of u v at the equator height of the turbine occurs between 4.5 and 5D (4.9D for positive values and 4.7D for negative values) downstream of the turbine, which is about 1D farther downwind compared to the maximum TI point. It should be noted that, again in this figure, only a slight lateral asymmetry (in terms of |u v |) can be observed.

Summary
Acknowledging the prospects of VAWTs as alternative wind energy extractors along with HAWTs in a future clean-energy outlook, which is likely to be marked by diversity, targeted research on VAWTs' performance is deemed to be highly useful and necessary. One of the research targets, which is especially crucial in designing potential VAWT farms, is to characterize VAWT wakes; a target which is still considerably underachieved for VAWTs, particularly with respect to HAWTs. In this view, one of the approaches that can greatly contribute to the cause is to use turbulence-resolving numerical simulation techniques, which can provide plenitude of high-resolution spatial and temporal information about the flow field and lead to valuable insight into the behavior of the turbine wake.
In this study, we used a previously-validated large-eddy simulation framework, in which an actuator line model is employed to parameterize the blade forces on the flow, to simulate the atmospheric boundary layer flow through stand-alone VAWTs placed on a flat terrain. For a typical straight-bladed 1-MW VAWT rotor design, first, the variation of the power coefficient with the tip-speed ratio and the chord length of the blades was studied. In doing so, the optimum combination of TSR and solidity (Nc/R), which yielded the maximum power coefficient of 0.47, was found to be 4.5 and 0.18, respectively. Second, for a VAWT with this optimum combination, a detailed study on the characteristics of its wake was performed, in which different mean and turbulence statistics were inspected. The mean velocity in the wake was found to need a long distance to recover; for example, the wake requires a distance of 14 rotor diameters to recover its center velocity to 85% of the incoming velocity. It was also seen that for this case, the point with the maximum velocity deficit is located 2.7 rotor diameters downstream of the center of the turbine (at the equator height of the turbine), and only after this point, the wake recovery starts with a rate (based on the change of the maximum velocity deficit) that is decreasing with streamwise distance. The turbulence intensity was observed to reach its maximum value (at the equator height of the turbine) 3.8 rotor diameters downstream of the VAWT. As we go towards the upper and lower extremities of the rotor, the height-specific maximum of the TI moves closer to the turbine and its value also increases. Turbulent momentum fluxes, which are a gauge for flow entrainment and, as a consequence, are responsible for the recovery of the wake, were also quantified, and it was shown that in the equator height of the turbine, the magnitude of the lateral flux peaks about 1D farther downwind of the maximum TI point. The above-mentioned mean and turbulence statistics corresponding to the optimum tip-speed ratio show only slight lateral asymmetries in the wake. However, significant vertical asymmetries were observed in terms of both the TI and magnitude of momentum fluxes, with higher values at the upper edge of the blades compared to the ones at the lower edge.
This study paves the way to further explore VAWT wakes and to discover the effects of different relevant parameters on the wake behavior. Moreover, it can serve as a solid foundation for future studies on performance, characteristics and optimization of VAWT farms.
Supplementary Materials: Zenodo DOI:10.5281/zenodo.51387 (https://zenodo.org/record/51387). Video S1: Normalized instantaneous streamwise velocity field both on a vertical plane (x-z) going through the center of the turbine and on a horizontal plane at the equator height of the turbine (Note: the physical time corresponding to this video is 1 minute and 17 seconds, and the size of the blades is magnified for illustration purposes). Video S2: Normalized instantaneous streamwise velocity field on a horizontal plane at the equator height of the turbine for two cases: when the turbine starts to operate (top) and when the flow has reached statistically steady condition (bottom) (Note: the physical time corresponding to both videos is 1 minute and 17 seconds, and the size of the blades is magnified for illustration purposes).

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this Appendix, the procedure of the method with which the dynamic stall phenomenon is modeled is described in detail. The dynamic stall model is based on the modified MIT model developed by Noll and Ham [27], which is a practical modification of the original MIT model [30]. This model has the advantage of being simple and easy to use and also has been found to work better for VAWTs compared to other available models [1]. It is noteworthy that the following procedure can be implemented for both VAWTs and HAWTs.
Dynamic stall is a phenomenon that occurs for an airfoil when the angle of attack of the incident flow keeps changing with time and its rate of change (i.e.,α = dα dt ) is sufficiently large. For a blade element of a turbine (either VAWT or HAWT) (placed in a turbulent flow), the change of α with time can be originated by three main sources: (1) the turbulent fluctuations of the incident flow; (2) the changes (spatial or temporal) in the mean incident flow; and (3) the rotation of the blades. Of these three reasons, the second one is normally specific to HAWTs, since an HAWT blade element experiences the variation of the boundary layer mean velocity profile at different heights; which is not the case for a VAWT blade element, as it moves at a constant height. However, the third reason is specific to VAWTs, because the geometry of a VAWT rotor is such that α (for a given blade element) oscillates between a maximum positive value and a minimum negative value in each revolution (even with a uniform inflow); however, for an HAWT blade element, assuming a uniform mean inflow, α remains constant during one revolution. Since the MIT model (and other similar practical models) is (are) only appropriate for the large-scale behavior of α in time, in our implementation of this model, the dynamic stall effects arising from the above-mentioned second and third sources, as well as the relatively large-scale turbulent fluctuations (from the first source) are modeled, while the changes of α arising from the relatively small-scale turbulent fluctuations of the incident flow are filtered out.
In order to implement the above-mentioned procedure,α is calculated from a time-averaged and smoothed curve of α f = α f it (θ) during one revolution. For this purpose, the angle of attack at each azimuthal angle is time-averaged during each N rev revolutions of the blades, and then, a polynomial curve, α f it (θ), is fitted on the time-averaged curve, α avg (θ). For the rest of the dynamic stall calculations, it is the α f it (θ) curve that is used. Figure A1 shows an example for this procedure for TSR = 2 and c = 2 m. The azimuthal angle, θ, is considered to increase counterclockwise (when seen from above) from −90 • to 270 • , in a way that θ = 0 • and θ = 180 • correspond to the most downstream and the most upstream points of the rotor, respectively. It is also noteworthy to mention again that the sense of the rotation of the turbine blades is counterclockwise when seen from above. Here, for the curve fitting, an eighth order polynomial is used to detect the two extrema accurately. Subsequently, we implement the modified MIT model on the α f it (θ) curve and construct C L,DS (α) and C D,DS (α) curves, which are lift and drag coefficients as a function of the angle of attack considering dynamic stall. In the modified MIT model, we use the tabulated airfoil data for lift and drag coefficients, and based on that, C L,DS (α) and C D,DS (α) are constructed. Based on the tabulated airfoil data, we can determine the static stall angle, α SS > 0, and the lift coefficient at static stall, C L,SS > 0. The static lift and drag coefficient functions derived from the tabulated airfoil data are designated as C L,table (α) and C D,table (α) hereafter. Moreover, the slope of the C L,table (α) curve before the static stall can be calculated as a s = C L,SS /α SS , considering that in this region, normally, C L,table (α) is linear.
As can be seen in Figure A1, the global (i.e., the curve-fitted) behavior of |α| in one revolution of a blade element is such that |α| twice (once for positive α values and once for negative α values) increases from zero to a maximum value and then decreases to zero again. In each of these increase-decrease cycles of |α|, the MIT dynamic stall model casts the flow in one of the four below dynamic stall states: State 1 occurs when |α| ≤ α SS . In this state, both lift and drag coefficients are extracted directly from the static tabulated airfoil data: State 2 occurs when α SS < |α| < α DS and α fα > 0 (i.e., |α| is increasing in time). α DS is calculated with the following formula: where c is the blade chord length, V rel is the magnitude of the relative velocity (which is also a function of the azimuthal angle),α = Ωdα f it /dθ, Ω is the angular velocity of the blade and γ is a constant that has a dimension of an angle and is weakly a function of the airfoil type and is determined experimentally [27]. If an experimental value for γ is not available, a value of one radian is recommended [31]. We keep calculating α DS in this state, until the point at which |α| is on the verge of becoming larger than α DS (i.e., the point at which the model goes to State 3). We designate this last value of α DS as α DS, f inal , and with this value, we calculate the maximum value of C L,DS (i.e., C L,max ): and we apply the following clipping conditions on C L,max : Throughout this state, the lift coefficient is extrapolated from static values, and the drag coefficient is still directly extracted from the static tabulated data: where (as in Noll and Ham [27]) a sine function is used for extrapolation (noting that in the range of angles of attack, on which we normally apply the model, |α| is small, and we have sin(|α|) ≈ |α|).
State 3 occurs when α DS, f inal < |α| and α fα > 0 (i.e., |α| is still increasing in time). As soon as the model enters State 3, we start to calculate the elapsed time from the moment in which State 3 is triggered; in other words, we start to calculate the time elapsed after the α DS, f inal value has been reached; we call this time t DS . In this state, the lift and drag coefficients are calculated as: However, in this state, we only keep using Equations (A8) and (A9) as long as these conditions are both satisfied: C L,DS ≤ C L,max and t DS V rel /c < 1; otherwise, we set the lift coefficient to the C L,max value and calculate the drag coefficient accordingly (as shown below). We designate the value of |α| of the moment in which either of the aforesaid conditions is on the verge of being violated as α C L,max .
If C L,DS > C L,max Or t DS V rel c ≥ 1 : State 4 occurs when |α| > α SS and α fα ≤ 0 (i.e., when |α| starts to decrease with time). We designate the azimuthal angle of the moment in which |α| starts to decrease as θ α max . At this stage, C L,DS is lowered exponentially (in time) from C L,max to C L,SS .
where R is the radius of the blade element about the axis of rotation (in the case of a VAWT, R is simply the radius of the VAWT rotor).
As can be noticed in the above procedure, α(t) (i.e., α(θ)) needs to be a smooth function for the above model to work. Because of this, we use α f = α f it (θ) (i.e., the time-averaged and curve-fitted value of α) in the above procedure instead of α. Thus, at the end of each N rev revolution and after getting the α f it (θ) function, we apply the MIT model on this curve, and we construct the C L,DS (α) and C D,DS (α) functions, which will be used in the next N rev revolutions. For the first N rev revolutions (for which we still do not have α f it (θ)), one can preliminarily just use the static tabulated airfoil data. Figure A2 shows an example of a constructed C L,DS (α) curve under dynamic stall, which corresponds to the α f it (θ) shown in Figure A1. The aforesaid states of the model are shown in the figure. As can be seen in this figure, to construct this curve, we need both the tabulated airfoil data and some parameters, which we should obtain from the MIT model. As a summary, all of the necessary data and parameters required to construct the C L,DS (α) and C D,DS (α) curves are listed below: (1) Tabulated airfoil data: C L,table (α) and C D,table (α); (2) α SS ; (3) C L,SS ; (4) a s ; (5) γ; (6) α DS, f inal ; (7) α C L,max ; (8) C L,max ; (9) θ α max It should be noted that for the last four items, two values are obtained for each revolution: one for the α ≥ 0 cycle and one for the α < 0 cycle. In our simulations, we have used γ = 1 radian and N rev = 15.
A comprehensive and step-by-step procedure to implement the MIT dynamic stall model is given in the flowchart of Figure A3; the flowchart is deliberately given in a way that can be followed for coding purposes in any common programming language.