Simulating Dynamic Stall Effects for Vertical Axis Wind Turbines Applying a Double Multiple Streamtube Model

The complex unsteady aerodynamics of vertical axis wind turbines (VAWT) poses significant challenges to the simulation tools. Dynamic stall is one of the phenomena associated with the unsteady conditions for VAWTs, and it is in the focus of the study. Two dynamic stall models are compared: the widely-used Gormont model and a Leishman–Beddoes-type model. The models are included in a double multiple streamtube model. The effects of flow curvature and flow expansion are also considered. The model results are assessed against the measured data on a Darrieus turbine with curved blades. To study the dynamic stall effects, the comparison of force coefficients between the simulations and experiments is done at low tip speed ratios. Simulations show that the Leishman–Beddoes model outperforms the Gormont model for all tested conditions.


Introduction
The need for sustainable energy sources with no carbon dioxide emissions has caused an increase of interest in wind power.Today, the vast majority of wind power plants are horizontal axis wind turbines (HAWT).However, during the last few years, interest in vertical axis wind turbines (VAWT) has increased, due to their potential to reduce the energy cost [1,2].For large turbines, the concept of VAWT can be beneficial compared to HAWT [1,3].The VAWT size can be adjusted to optimize the absorbed power, i.e., both the diameter and height of rotor can be changed.Moreover, the lower center of gravity (compared to HAWT) makes the VAWT design favorable for floating platforms [4].In addition, VAWTs are omni-directional and therefore do not require a yawing system.The generator for a VAWT does not need to be optimized for weight and size, since it can be located at the ground level, which simplifies installation and maintenance.Several projects on multi-megawatt VAWTs for offshore use are being carried out [5][6][7].
A design with less moving parts is the principal advantage of VAWTs.However, cyclic loads on the blades, tower and drive train together with complex unsteady aerodynamics are of a major concern for VAWTs [1].This study focuses on modeling the unsteady aerodynamics.The simulation models have to account for the continuous changes of the angles of attack, which are natural during VAWT operation.The amplitude of the angles of attack increases with decreased tip speed ratio (TSR), and it is at the low TSRs where the dynamic stall is present [8].
VAWTs are usually stall regulated due to the complexities with pitch regulation [9,10].Therefore, the effect of dynamic stall is very important for the operation of VAWTs.For high wind speeds, where the rotational velocity is limited, the TSR decreases, which causes an increase in the amplitude of the angles of attack (for VAWTs with fixed blades).Consequently, the blades fall into a stall, causing a drop in lift force.Thus, VAWT blades are passively controlled by the dynamic stall, and therefore, it is important to account for it when modeling the turbine.Moreover, the amplitude of the oscillations of the blade forces changes with dynamic stall, which is related to fatigue problems.
Several approaches exist to model the effects of dynamic stall.An overview of the dynamic stall models adapted for the VAWTs applications can be found in [11,12].The Gormont model [13], modified by Massé and Berg [14,15], is widely used for VAWTs simulations due to its simplicity of implementation [16][17][18][19][20].A more complex method was developed by Beddoes [21], which models both attached and separated flows.The model was further modified by Leishman [22,23], and therefore, it is commonly called the Leishman-Beddoes model.This model has shown better agreement (compared to the Gormont model) for the experiments on pitched airfoils [24].Both the Leishman-Beddoes and Gormont models were originally developed for helicopter rotors and were later adopted for the simulations of VAWTs.
Studies on the VAWT blade loads during the dynamic stall are rare, due to the lack of existing experimental data.Scheurich et al. [25] have conducted a study in which the previously-measured data on a VAWT at the Reynolds number Re = 40, 000 (the data were taken from [26]) was compared against the vortex model coupled with the the Leishman-Beddoes dynamic stall model.Ashuri et al. [27] have used a double multiple streamtube model (DMS) with the Gormont dynamic stall model to assess their own experimental data on a vertical axis marine current turbine with three straight blades.
The present study focuses on VAWTs operating at low TSRs, where the dynamic stall is present.The simulation results are performed with the DMS model and assessed against the experimental data on the Sandia 17-m turbine with two curved blades [28].The experimental data on this turbine were chosen, because of its large scale size, which gives the Reynolds number Re > 1 × 10 6 .Two dynamic stall models are compared: the Gormont model, with the modifications of Massé and Berg, and the Leishman-Beddoes model, modified by Sheng et al. [29].
The aim of the work is to compare two dynamic stalls combined with the DMS model against measured data on a VAWT.Normal and tangential force coefficients at the mid-span of the rotor are studied, and the accuracy of the models is assessed.

Method
The DMS model, which is built on blade element momentum theory, is chosen for the turbine modeling.The turbine is divided into a number of streamtubes, and the principle of momentum conservation is used for each streamtube.The DMS model is a fast method to calculate the aerodynamic performance of a VAWT, and it has shown reasonable agreement with experiments on Darrieus turbines with curved blades [12].Two versions of the DMS model were developed around the same time, one by Paraschivoiu [30] and one by Read and Sharpe [31].The difference is that the Paraschivoiu model assumes that the streamlines are parallel to each other, and the Read and Sharpe model includes the streamtubes expansion in the horizontal plane.In the present study, the DMS model is formulated as in [32], and the flow expansion model is included, as by [31].A database on the static lift and drag coefficients is required for the DMS model.The widely-used database for symmetrical airfoils, formed by Sandia National Laboratories [33], is adapted for the present model.

Leishman-Beddoes Model
The Leishman-Beddoes model (LB) [22,23], which was further developed by Sheng et al. [29], is used in this work.It is described in detail in [24], and the main principles of the model are reviewed here.
The LB model consists of three parts: unsteady attached flow, stall onset and separated flow.The unsteady attached flow solution is comprised of circulatory and impulsive loadings, which are caused by the time varying bound vortex.The circulatory normal force coefficient is expressed as: where C Nα is the slope of the static normal force coefficient at certain a Reynolds number, and α En is an equivalent angle of attack: where α is the geometrical angle of attack and X n , Y n and Z n are the deficiency functions, which are empirically derived based on the flow velocity and the pitch rate, and they can be found in [24].As seen from Equations ( 1) and (2), the LB is the iterative model with indexes n and n − 1 representing the current and previous time steps.A delayed angle of attack due to the lag in pressure response is calculated in the unsteady attached flow part, as well: where D α is the deficiency function: where T α is an empirically-derived time constant, which is T α = 5.78 for the NACA (National Advisory Committee for Aeronautics) 0015 airfoil, and ∆s is the non-dimensional time-step: where V is the incoming flow velocity and c is the chord length.Due to the flow reversal within the boundary layer, a leading edge vortex forms at the airfoil surface.The critical angle of attack α cr is used to define the condition at which the dynamic stall may begin: where the reduced pitch rate r n is: where α is the pitch rate and r 0 is the reduced pitch rate, which delimits the quasi-steady stall and the dynamic stall, which is r 0 = 0.01.α ss and α ds0 are the static stall-onset angle and critical stall onset angle correspondingly, which are α ss = 14.67 • and α ds0 = 17.81 • for NACA0015.The following dynamic stall condition is used: The effects of separated flow are divided into two groups: trailing edge separation and leading edge vortex convection.The trailing edge separation is associated with the time delay in the movement of the boundary separation point, and it is obtained via Kirchhoff's flow approximation: where f is the delayed separation point and α 1 , S 1 and S 2 are the constants based on the airfoil profile and the local Reynolds number, found in [24].The boundary layer around the blade itself is time dependent, and this effect is superimposed on the pressure response delay, which is represented by α ; Equation (3).This additional delay is represented by a dynamic separation point, f : where the deficiency function D fn is: where T f is an empirically-derived time constant, which is T f = 3.The normal force coefficient for the unsteady conditions before the dynamic stall onset is obtained as follows: After the stall condition is met, the leading edge vortex convects over the surface of the airfoil towards the trailing edge.During this process, a significant increase in normal force is present: where C v N is the normal force coefficient during the vortex convection (so-called "vortex lift"), which depends on the pitch rate, and V x and B 1 are the parameters based on the local Reynolds number and airfoil profile, found in [24].The normal force decreases rapidly when the vortex passes the trailing edge.The total normal coefficient is estimated as: The tangential force coefficient is obtained via Kirchhoff's flow relation using the dynamic separation point: where η and E 0 are the empirical constants, and they are η = 1 and E 0 = 0.25 for the NACA0015 airfoil.Figure 1 shows an example of the force response of a pitching blade, simulated by the the LB model, and the features described above are present.Previous studies has shown that two iterations are sufficient to approach convergence for the LB model [11,24].Therefore, the LB model runs two loops during each iteration of the DMS model.

Gormont Model
The Gormont model [13] was developed earlier than the LB model, and several modifications were applied to adapt the model for VAWTs.The version of Massé and Berg [14,15] is used in this study, and only the fundamental principles of the model are presented.For the detailed description of the Gormont model, the reader is referred to [30].
The Gormont model empirically mimics a pitching airfoil by modifying the geometrical angle of attack into the reference angle of attack, α ref : where α is the geometrical angle of attack and K 1 and ∆α are the empirically derived parameters, which are found in [30].
The reference angle of attack is used to calculate the dynamic lift and drag coefficients: where α 0 is the zero-lift angle of attack, and parameter m is: where α ss is the static stall angle and C L and C D are the static lift and drag coefficients (at a certain angle of attack).
In the version by Massé and Berg, the lift and drag coefficients are further modified as follows: where A M = 6 is the empirical constant.Equations ( 16)- (19) show that the Gormont model is the empirical model, and the main features of the dynamic stall are not implemented.Namely, the trailing edge separation point, the dynamic stall condition and the leading edge vortex convection are not defined in the Gormont model.However, its simplicity makes the Gormont dynamic stall model common for simulations on vertical axis turbines.

Vortex Shedding
The effects of dynamic stall are more complex during the turbine operation than during the pitching of a single blade.Figure 2 shows the vortex shedding structure of a Darrieus turbine at low TSR.This figure is based on velocity measurements of a straight-bladed Darrieus turbine operating at a TSR of 2.14 in a water channel, obtained by Brochier et al. [34].It is shown that both leading and trailing edge vortices are released and swept away when the blade passes Quadrant III; thus, the flow is fully separated.This can be both due to the circulatory motion of the blade and due to the highly turbulent flow, as noted in [34].An additional modification to the LB model is applied to account for the vortex shedding.This is done by modeling the fast release of trailing and leading edge vortices: the delay in the angle of attack, leading edge separation point and vortex lift are set to zero when the blade passes Quadrant III: The effect of vortex shedding is not implemented in the simpler Gormont model, as it does not track the leading edge vortex.

Flow Expansion
Since the flow slows down at the downwind part of the turbine, there will be an expansion of streamtubes due to the continuity principle, and the streamlines will no longer be parallel to each other.This is illustrated in Figure 3.The flow expansion model by Read and Sharpe [31] is used in the present work, and the change in streamlines direction is taken into account.The model assumes that the flow expands in the horizontal plane.In Figure 4, the cross-section of a streamtube is shown.θ is the azimuth angle of the center of streamtube, w is the width of a streamtube, γ is the angle between the vector of free-stream velocity, V ∞ , and the streamline V is the flow velocity at the turbine.Subscripts up and dw stand for the up-and down-stream parts of the rotor.V up and V dw are calculated within the DMS model.The main principle behind the flow expansion model is the use of the continuity equation: An iterative process is required to find the size of each streamtube for known V up and V dw .θ − γ instead of θ was used in the DMS model in the expressions for relative wind velocity, angle of attack and the velocity slowdown to account for the change in the streamlines direction.
Cross-section of a streamtube.

Flow Curvature
The original DMS model assumes a blade cross-section as a point, and all of the equations are derived using this assumption.However, there is a variation of the instantaneous flow velocity along the chord, indicating that the flow is of a curvilinear nature [35].Consequently, the angle of attack is not constant along the chord length.To account for the curvature effects, the angle of attack is modified further.The assumption of the two-dimensional potential flow around the flat plate and the Joukowski transformation together with the Kutta condition were used.The full derivation is documented in [32], and the expression for the modified angle of attack is: where δ is the blade pitch angle, Ω is the turbine rotational speed, r is the turbine local radius, x 0 is the normalized blade attachment point, c is the blade chord, V rel is the relative flow velocity at the blade and i denotes the streamtube index.The last two terms in Equation ( 22) account for the chord-to-radius ratio and for the blade mounting position.

Simulation Parameters
The Sandia 17-m turbine is studied in this work.It is a two-bladed Darrieus turbine with curved blades attached at both ends.The blade profile is NACA0015 with a chord length of c = 0.61 m.An illustration of the turbine is shown in Figure 5, and the reader is referred to [36] for a detailed description of the rotor design.The surface pressure measurements for this turbine were performed by Sandia National Laboratories and are documented in [28].The pressure transducers were installed at the mid-span of the blade.Tangential and normal force coefficients, C T and C N , were obtained by using the pressure distributions and the measured incident wind speed at the blade, V rel .C T and C N are defined as: where F T and F N are tangential and normal forces per unit span and ρ is the air density.The sign convention of the forces and the definition of the rotor azimuth angle are shown in Figure 6.The turbine rotational speed was kept constant at Ω = 38.7 rpm during the experiments.The current study is in 2D, as the measurements were done at the mid-disk of the turbine, and the average Reynolds number is 1.4 × 10 6 .The mid-disk is divided into 120 streamtubes, 60 on both the upwind and downwind sides of the disk.The simulations are performed for the range of TSRs from 2.20 to 3.09 to cover the dynamic stall region.The flow expansion and the flow curvature effects are taken into account as described above.

Results and Discussion
This section presents the assessment of the simulated results compared to the measured data, both without dynamic stall modeling and with the LB and Gormont models.The simulation results without dynamic stall modeling are denoted as "no DS".The normal and tangential force coefficients at the blades' mid-span are studied.General comments on the simulation model are given at the end of this section.

Assessment of the Model
In Figure 7, the models are tested for the deep stall conditions at a TSR of 2.20.From the C N values in Figure 7a, it is noted that the blade is stalling at the azimuthal position of θ = 70 • .However, the results of both the LB and Gormont models predict stall earlier, at θ = 60 • , and the slope of the simulated C N -curve at θ < 45 • is higher than in the experimental data.The higher C N -slope is also present at TSRs of 2.33 and 2.49 for the simulations with the LB and Gormont models and without dynamic stall modeling.Therefore, the deviation in the C N -curve at θ < 45 • may come from the static lift and drag data, which are used by both the dynamic stall models and by the DMS model.However, the further overshoot of the C N response at 45 • < θ < 70 • cannot be the issue of the dynamic stall models, as it also takes place at the higher TSR.Thus, the authors presume that this overshoot originates from the DMS model.
There is a peak in the experimental data of C N at 213 • < θ < 250 • , centered at θ = 225 • , and both dynamic stall models do not reproduce it.This peak is likely to be the collision of the blade at the downwind side with its released vortex from the upwind side.This wake interaction has been noted by Akins [28], when the experiment was conducted, and it can also be seen in Figure 2. To hit the blade at θ = 225 • , the impacting vortex would originate at θ = 135 • ; Figure 2. Indeed, there is a drop in C N response at θ = 135 • , which should correspond to the vortex release.However, since the vortices are not modeled by the DMS, the wake interaction cannot be reproduced.
The Gormont and LB models show a similar response in C N in the second and fourth quadrants, while the C N response without dynamic stall modeling has peaks at 150 • < θ < 165 • and 290 • < θ < 310 • .Those peaks are present, because the static lift force coefficient increases instantaneously when the blade returns from the stall, and since the dynamic stall model is not applied, the delay in the lift force is not modeled.The difference between the LB and Gormont models' response in the third quadrant is due to their ability to represent the circulatory motion of the blade.It was noted in Section 2.1.3that at the third quadrant, both leading and trailing edge vortices are swept away, and the flow around the blade is fully separated.This is taken into account by the LB model, but it cannot be taken into account by the Gormont model, where the the vortex separation point is not modeled.
The tangential force coefficient C T is shown in Figure 7b.The experimental values of C T do not cross or even approach zero values at the azimuthal angle of θ = 180 • , when the blade moves directly downwind, and the angle of attack is zero.This is due to the error in the response of the pressure sensors, which were used in the experiments, and it is reported by Akins [28] as the measurements error.Such a behavior of the measured C T at θ = 180 • is present at all tested TSRs; see Figures 7-10.Therefore, the simulated results should not be considered erroneous in this region.The slope of the simulated C T -curve is shifted relatively to the measured C T in the upwind region.This displacement is present for all tested TSRs and is predicted by both dynamic stall models.Thus, it is presumed that it is an issue of the DMS model.The effect of the wake interaction observed in the C N data at 213 • < θ < 250 • is represented by the drop in the C T in the same region.As with the C N response, the wake interaction cannot be represented for the C T with the present simulation model.The ability of the LB model to trigger the release of trailing and leading vortices is beneficial when comparing to the result of the Gormont model at 180 • < θ < 270 • .An inconsistent behavior of the C T -curve by the Gormont model is observed in the region of 150 • < θ < 185 • , and it is present at the second quadrant for all TSRs.Such a response is due to the high dependence of the Gormont model on the static data; see Section 2.1.2.
The predicted C T values without dynamic stall modeling have two peaks at the upwind (at 55 • < θ < 60 • and 155 • < θ < 160 • ) and two peaks at the downwind (at 200 • < θ < 215 • and 290 • < θ < 310 • ) at TSRs of 2.20, 2.33 and 2.49; Figures 7-9.These peaks are present, because both the stall delay and flow reattachment are excluded when predicting the results without dynamic stall modeling.
The results at TSRs of 2.33 and 2.49 are shown in Figures 8 and 9, respectively.The wake interaction at TSRs of 2.33 and 2.49 is not as noticeable as at a TSR of 2.20, where it is observed on the measured C N data at the third quadrant.The agreement between the C N predictions and the measured C N is better with increasing TSR.The models error is largest at a TSR of 2.20; see Table 1, corresponding to the influence of the wake interaction, which cannot be modeled.Similar to the C N results, the agreement between the measured and predicted C T data is better with increasing TSR.The method of the comparison of the models' accuracy is discussed further in Section 4.2.In Figure 10, the results for a TSR of 3.09 are shown.In this condition, the stall at the upwind is less pronounced, and the wake interaction at the downwind is no longer evident; Figure 10a.However, the overshoot in the predicted C N is still present, similar to the one at the lower TSRs, and it is discussed earlier in this section.At the third quadrant, the C N values by the LB model are closer to the experimental data than the Gormont model due to the separated flow in this region, though it is not as pronounced as for the lower TSRs.In Figure 10b, the shift in the C T curve at the first quadrant is more evident than at the lower TSRs, and this is due to the streamtube expansion model, as mentioned earlier.The inconsistent behavior of the C T response at the second quadrant, noted at the lower TSRs, is still distinct.At a TSR of 3.09, the predictions without dynamic stall modeling show peaks that are observed neither on the measured data nor on the simulations with the LB and Gormont dynamic stall models.These are the peaks on the C N and C T data in the second quadrant, and they are due to the excluded delay in the stall onset and flow reattachment, which is also noted at the lower TSRs.A similar behavior of the predicted force response without dynamic stall modeling is noted at the lower TSRs, and it is discussed earlier in the section.

Final Discussions
The LB model shows better agreement with the measured data compared to the results by both the Gormont model and the modeling without dynamic stall.While the C N results by all three models are similar (except the peak at the second quadrant without dynamic stall modeling), the C T results show the significant benefit of the LB model.Moreover, the Gormont model does not show the improvement of the C T response compared to the modeling without dynamic stall.
The quantitative comparison of the results for all of the TSRs is shown in Table 1.The RMS of the difference between the simulated results and the experimental data is presented, which is a simple method to evaluate the models' accuracy.The disadvantage of using the RMS is that only the deviation in the force response values is taken into account, i.e., only the difference along the y-axis is considered.If the curve of predicted values is shifted along the x-axis, but has the same shape as the curve of the measured data, the RMS can still be as high as for an oscillating curve of a different shape.When comparing the C N results by the LB model and without dynamic stall modeling, it is observed that the RMS error of the C N without the dynamic stall model is lower than with the LB model, in the case of T SR = 3.09.The peak of C N at θ = 150 • without dynamic stall modeling has no physical meaning (as there is a delay in the flow reattachment); however, since this peak is closer to the measured data, it happens that the error is smaller.Thus, the inconsistent behavior cannot be quantified with the RMS, and more advanced methods for the global comparison of the models should be used, which is not studied in this work.
When designing the turbine, the maximum value of the normal force over the revolution is used to represent the loads and to dimension the rotor and tower.The mean of the tangential force corresponds to the power output.However, the maximum values of the C N and the mean of the C T do not reflect the oscillations during revolution, and the unsteady effects (including dynamic stall and vortex shedding) can be omitted, while the maximum C N and the average C T can be close to those at steady conditions.Considering the problems associated with fatigue, it is utterly important that the model is stable and does not create inconsistent oscillations.In this light, the LB model is beneficial over the Gormont model and over the case without dynamic stall modeling, though all of the models produce deviations between the predicted and measured data.

Conclusions
The double multiple streamtube model was used to simulate a vertical axis wind turbine at low TSRs, and the results are assessed against the experimental data in two-dimensional space.Two dynamic stall models were used: the Leishman-Beddoes-based model with improvements for the conditions of VAWTs and a Gormont model.Additionally, simulations without dynamic stall modeling were presented.For all test conditions, the Leishman-Beddoes-type model shows better agreement with the experimental data.Although the model does not reproduce all of the unsteady effects observed in the experimental data, it shows higher accuracy and is more stable compared to the Gormont model, which is widely used for simulations of vertical axis turbines.

Table 1 .
The root mean square (RMS) of the difference between the predicted and measured results.LB, Leishman-Beddoes; DS, dynamic stall.