Development of a Second Order Dynamic Stall Model

Dynamic stall phenomena bring risk for negative damping and instability in wind turbine blades. It is crucial to model these phenomena accurately to reduce inaccuracies in predicting design driving (fatigue) loads. Inaccuracies in current dynamic stall models may be due to the facts that they are not properly designed for high angles of attack, and that they do not specifically describe vortex shedding behaviour. The Snel second order dynamic stall model attempts to explicitly model 10 unsteady vortex shedding. This model could therefore be a valuable addition to DNV GL’s turbine design software Bladed. In this thesis the model has been validated with oscillating airfoil experiments and improvements have been proposed for reducing inaccuracies. The proposed changes led to an overall reduction in error between the model and experimental data. Furthermore the vibration frequency prediction improved significantly. The improved model has been implemented in Bladed and tested against small scale turbine experiments at parked conditions. At high angles of attack the model looks promising for reducing 15 mismatches between predicated and measured (fatigue) loading. Leading to possible lower safety factors for design and more cost efficient designs for future wind turbines.


Introduction
Wind turbines operate in highly unsteady aerodynamic environments (Leishman, 2002). For design and certification, design load cases (DLCs) have been set which describe the conditions that wind turbine designs have to withstand (DNV GL, 2016). Some of the design driving DLCs are those for parked and idling conditions where wind turbine blades will experience high angles of attack (AoAs), leading to (dynamic) stall behaviour (Schreck et al., 2000). The wind turbine yaw angle is defined as the angle in the horizontal plane between the free-stream wind direction and the wind turbine rotor shaft. It can be noted that when the turbine is parked, and the blades are pitched to 90 • , the yaw angle effectively becomes the inflow angle. When the yaw system is not operating during parked conditions due to a failure, the blades will experience flow from all directions. For particular wind directions the flow on the blades is separated, leading to dynamic stall effects. These effects may already appear at inflow angles of below 30 and 20 • (Gonzalez and Mund-uate, 2007). Therefore, accurate modelling of dynamic stall is therefore crucial in wind turbine design (Choudry et al., 2014).
Dynamic stall is a phenomenon leading to larger variations in lift, drag, and pitching moments on the aerofoil than would be observed during steady operation (Choudry et al., 2014). This then creates larger aerodynamic forces on the blades than expected during steady conditions (Leishman, 2002). Dynamic stall happens with dynamic variation in the inflow and/or the effective angle of attack and can be viewed as a delay in the onset of stall. Recirculation of flow after the static stall angle starts near the trailing edge and rapidly moves towards the leading edge, leading to the formation of a large dynamic stall clockwise vortex at the leading edge at increasing angles of attack. The dynamic stall vortex will travel along the suction side, leading towards the trailing edge before detaching completely. Full separation will occur when the dynamic stall vortex is completely detached. This moment is called the "break" or "dynamic stall onset". As a re-sult, low lift remains until reattachment of the flow. However, a time delay for reattachment of the flow is present as well. After reattachment the process repeats, creating a hysteresis loop. Dynamic stall phenomena carry the risk of negative damping and instability, especially if the aerofoil is oscillating in and out of stall (McCroskey, 1981). A visual description for dynamic stall is presented in Fig. 1.
When keeping the aerofoil pitched in (deep) stall for longer periods of time, periodic vortex shedding will occur. A single large dynamic stall vortex will no longer be shed, but rather multiple periodic vortices from both the leading and trailing edge will be shed. This will induce time-varying loads on the blades (Riziotis et al., 2010). The periodic vortex shedding is characterized by the Strouhal number representing the dimensionless frequency of shedding (Pellegrino and Meskell, 2013). The Strouhal number is defined following Eq. (1): in which f notes the characteristic vortex shedding frequency, c the aerofoil chord (sometimes the projected chord length perpendicular to the incoming flow), and U the wind velocity at the wind turbine blade section. Synchronization of the natural and Strouhal frequencies (a "lock-in") will lead to resonance (Pellegrino and Meskell, 2013). Locked-in vortex-induced vibrations are a potential threat in standstill conditions as the turbine size is increasing. Modern aeroelastic tools with dynamic stall models are only able to provide accurate deep stall loads at conditions close to maximal lift, so relatively small angles of attack (Riziotis et al., 2010). The same study noted that today's aeroelastic tools are not properly tuned for high angles of attack. Mismatches between load predictions between measurements and engineering tools have been found to be as high as 20 % for highwind-speed dynamic stall conditions (Schreck, 2002). Stateof-the-art aerodynamic models overestimate fatigue loading by 15 % (Schepers and Snel, 2007). Madsen et al. (2019) show promise for using computational fluid dynamics (CFD) in early stages of wind turbine design and shape optimization, and Sørensen et al. (2016) showed that CFD models had good agreement with instrumented rotor experiments. So although CFD is becoming more available and useful in wind turbine design, it requires large computational power. Therefore these tools are not yet fit for practical design calculations. In the industry there is a need for relatively fast and accurate engineering models predicting key loading. Hollierhoek et al. (2013) studied different dynamic stall models, namely Beddoes-Leishman, ONERA, and Snel models. However, a clear single best model was not found. Gonzalez and Munduate (2007) and Wala et al. (2018) both showed promising results using modified and optimized Beddoes-Leishman models compared with experimental data. Inaccuracies in dynamic stall models may be due to the abovedescribed fact that they are not properly designed for high angles of attack and that there are some of them that do not specifically describe vortex shedding behaviour. The Snel second-order dynamic stall model does attempt to explicitly model unsteady vortex shedding. However, there is still the need for further tuning and validation of the model (Snel, 1997). Snel takes the model from Truong (1993) as a starting point. Truong (2017) already showed promise by modifying his model for unsteady vortex shedding, and Snel takes a similar approach. However, as described in Sect. 2, the Snel model differs from the Truong model by incorporating the steady lift coefficient and is therefore interesting to study and modify. This paper will provide a detailed analysis of the Snel second-order model and will try to answer the following main research question: what are possible ways to improve predictions of blade vibrations during dynamic stall in parked conditions using the Snel second-order model? This paper will have the following outline.
-The Snel model is validated against experimental data.
-Proposed changes are presented based on the validation results to improve the model predictions. These include a dimensional analysis, calculation of the slope for the potential lift coefficient, application of the normal force coefficient, an investigation into the downstroke and vortex shedding predictions of the model, and finally an optimization for empirical constants.
-Attention will be paid to vibration prediction as this influences turbine fatigue loading, which has a large impact on the design of the wind turbine.
-An absolute error analysis is carried out before and after the improvements to quantify the increase in performance. - The improved model will be tested with actual smallscale turbine experimental data to assess performance in combination with Bladed.

Analysis and validation of the Snel second-order model
This section will validate the Snel model and propose adaptations to the model to improve the performance. The description of the model is based on the description in both Snel (1997) and Hollierhoek et al. (2013). Snel (1997) derived a dynamic stall model based on the work of Truong (1993), who proposed that the dynamic lift coefficient can be distinguished into two parts, namely c l 1 and c l 2 : c l,dyn = c l,steady + c l 1 + c l 2 .
The first describes the forcing frequency response and the second term the self-exited higher-frequency dynamics. Snel follows this approach but also expresses the first part as the difference from the steady-state (time-averaged) lift coefficient (Montgomerie, 1996). The dynamic lift of the Snel model will be as in Eq.
(2). The first part, c l 1 , must decay to zero when no excitation is present, while the second part will decay to zero for small angles of attack, but nearing stall the second part will show periodic oscillations related to vortex shedding. In the original model of Truong the first part is based on the Beddoes-Leishman dynamic stall model. The Snel model uses the SIMPLE model from Montgomerie (1996) as a departure point for the first-order correction while Truong (1993) uses the Beddoes-Leishman (B-L) model to calculate c l 1 . The modelling of the first part will therefore follow τ c l 1 + cf 10 c l 1 = ft 1 . (3) The forcing term ft 1 will be based on the time derivative of the difference between the steady and potential lift coefficients. This is shown in Eqs. (4) and (5). The function is made non-dimensional by taking the coefficient of the derivative term as the time constant usually used in dynamic stall. This constant is described in Eq. (11).
The coefficient of c l 1 can be seen as a spring trying to pull the term back to the steady state. The stiffness of the spring of this equation is given by Eq. (6). In downwards pitching motion the stiffness is higher than in upwards pitching motion.
The second part of the dynamic lift coefficient is of the second order to create the higher-frequency dynamics. This will be a non-linear mass-damper-spring system following τ 2 c l 2 + cf 21 ċ l 2 + cf 20 c l 2 = ft 2 . The spring and damping coefficients are taken from Snel (1997) and are defined as in Eqs. (8) and (9) respectively.
cf 20 = k 2 s 1 + 3 c l 2 2 1 + 3α 2 (8) Also, here it can be noted that the response will be different in pitching upwards or downwards. The forcing term is defined as In the second-order part of the model, k s is taken as the reduced vortex shedding frequency or Strouhal number. This is given a value of 0.2 as in the original model (Snel, 1997). In the above equations the time constant is given below and can be seen as the time it takes for the wind to travel half a chord distance:

Initial model implementation and validation
The Snel model as described in the section above is implemented in a MATLAB environment. The numerical implementation follows the steps described in De Vaal (2009). The time step used is 0.001 s to capture higher-frequency events from the Ohio State University experiments. The Ohio State University (OSU) experiments (Hoffmann et al., 1996) will be used for validation of the initial implementation of the Snel model. In the OSU experiments an extensive set of aerofoils have been tested for both unsteady and steady data. The measurements recorded responses to forced sinusoidal pitching oscillations. Different amplitudes, mean angles of attack, oscillation frequencies, and Reynolds numbers were tested. The focus in this paper will be on the NACA4415 and S809 aerofoils. The model parameters obtained from the OSU database, as displayed in Table 1, will be analysed. The oscillation frequencies of the cases are set such that the forcing angle of attack matches the OSU experiments. Figures 2  and 3 show the time series of the lift coefficient for cases with both low and high reduced frequencies and low and high mean angles of attack. The following observations can been made for the current model: -The current model overpredicts the loss of lift during the downstroke.
-The predicted vortex shedding does not always happen at the correct time.
-There is currently no unsteady vortex shedding at higher angles of attack. -The shedding frequency is not dependent on the reduced frequency while the experiments do show a dependency.
These points will be analysed as part of the possible improvements described in Sect. 3.3 up to Sect. 3.9.

Error analysis
To quantify the accuracy of the Snel model, the total absolute error between model predictions and experimental data is calculated. This will give an objective measure to assess proposed improvements. The Snel model is interpolated along the angle of attack to obtain the dynamic lift coefficient output at the precise angles of attack of the OSU experiment. The errors are then evaluated according to and the total absolute error of all data points is calculated using in which n denotes the number of data points of the OSU experimental data.

Proposed changes to the Snel model
This section will outline proposed modifications to the Snel model. The following areas will be investigated for modification:  a dimensional analysis of the model, the calculation method of the slope for the potential lift coefficient, application of the normal force coefficient instead of the lift force coefficient, the downstroke prediction of the model, prediction of vortex shedding, and lastly an optimization for empirical (aerofoilspecific) constants.

Dimensional analysis
It can be seen in the formulation of the model that Eqs. (8) and (10) are cast in a dimensional form. For different values for chords, wind speed, and pitching frequency, the current model will not produce identical results. In order to make them dimensionless, the time constant used in dynamic stall (Eq. 11) is added. The new constants will be such that the initial value of the constant is kept. The equations will now be cf 20 = k 2 s 1 + 3 c l 2 2 1 + 280 2 τ 2α2 , ft 2 = 0.1k s −0.15 c l,pot + 8τ ċ l,pot .
3.2 Correct slope for Cl potential The Snel model uses 2π as theoretical slope for lift coefficient at low angles of attack. This theoretical value might not be applicable to real aerofoils. Calculating the precise slope improves the accuracy of the model. Therefore the slope calculated from the aerofoil polar is used in the model. The slope is calculated between the intercept at angle of attack = 0 and the intercept at lift coefficient = 0.

Application of the normal force coefficient
The Snel model uses the lift coefficient of the steady aerofoil data. However, the lift coefficient tends to zero when angles of attack reach 90 • . At a 90 • angle of attack there is still unsteady vortex shedding present. Therefore, it would make sense to model vortex shedding to the normal force on the aerofoil instead of the lift force. The normal force coefficient together with the lift and drag coefficients are presented in Fig. 4. For implementation all the C L terms are changed to C n terms. A couple of additions must then be made to the model to obtain the dynamic lift and drag coefficients. First, the steady normal force coefficient (c n,steady ) is the sum of the steady lift and drag coefficient and the angle of attack following c n,steady = c l,steady cos(α) + c d,steady sin(α), and the steady chordwise coefficient (c c,steady ) will follow c c,steady = −c l,steady sin(α) + c d,steady cos(α).
When the dynamic normal force coefficient has been obtained, an inverse calculation yields the lift and drag coefficients: c d,dyn = c n,dyn sin(α) + C c,steady cos(α).

Downstroke of the model
The consistent differences between the implementation of the Snel second-order model and earlier implementations are lower values in the downstroke. Figure 5 shows the first-and second-order part of the model as a function of time. The second-order part ( c l 2 ) contributes highly to the lower values at the start of the downstroke, which causes a large part of the observed error. Equation (10) uses −0.15 c l,pot . Together with a negative ċ l,pot in the downstroke, ft 2 will be highly negative in the downstroke. To improve this behaviour the forcing term is set to zero for the downstroke. In addition, Figs. 2 and 3 show that the predicted shedding in the upstroke is larger than in the measurements. Therefore, the forcing term is lowered and will be changed to ft 2 = 0.1k s −0.08 c l,pot + 1.5τ ċ l,pot .
From Figs. 1 and 2 it is visible that there are higher frequencies in the downstroke. This is not modelled as cf 21 is a constant value in the downstroke; see Eq. (9). To allow shedding, the coefficient is changed to The damping in the downstroke is set higher than in the upstroke.

Prediction of vortex shedding
The shedding frequency will depend on the angle of attack. The current model does not predict a dependency and so it has to be improved in this aspect. The Strouhal number uses the chord or the projected chord perpendicular to the incoming flow. Because the projected chord length is driven by the angle of attack, it is proposed to add the projected chord length to the "spring" term (cf 20 ) of the second-order part. This allows for the desired dependency of stiffness. From Eq.
(1) it can be inferred that when projecting the chord perpendicular to the incoming flow this is effectively the same as projecting the Strouhal number. The new equation for cf 20 will now be cf 20 = 10 · k 2 s,pr 1 + 3 c l 2 2 1 + 280 2 τ 2α2 , with k s,pr = k s · sin(α).

Optimization for (aerofoil-specific) constants
An unconstrained minimization algorithm in MATLAB is used to optimize empirical constants. The algorithm searches for the lowest summation of the absolute errors, from Eq. (13), of all cases considered. The constants selected for this analysis are shown in the initial row of Table 2. The constants are selected since they are included in equations affected by the modifications from this section and also because they have a high influence on the output of the Snel model. Three optimization analyses have been carried out: a global optimization which covers all cases, an optimization which focussed only on the cases with a mean angle of  attack of 14 or 20 • , and an optimization on the low-reducedfrequency cases. The initial constants as in Table 2 are the start values for the optimization. The results for all individual cases, using the error analysis described in Sect. 2.2, are shown in Fig. 6. It can be seen that the global optimization is the most optimal. The global optimization gives the most consistently lower error compared to the initial model. In Bladed the Beddoes-Leishman dynamic stall model has three sensitive constants which are allowed to be changed by the user. For the Snel model, the optimization shows that two constants are the most sensitive. The optimization output showed significantly different values for these constants for the NACA 4415 and the S809 aerofoils. They are shown in the improved row of Table 2. These constants must be allowed to be changed by users of a turbine design code. The same goes for the Strouhal number. The constants are -C1 from Table 2, which will be called the first-order coefficient and Table 2, which will be called the second-order forcing coefficient.
The first-order coefficient C 1 will be 0.5 for the NACA 4415 aerofoil and 0.2 for the S809. For the second-order forcing coefficient C 2 will be 1 and 1.5 respectively. Table 2 shows that for Eqs. (18) and (19) the constants remain the same. Therefore these remain unchanged. Finally Eq. (21) will become as follows:

Results of the proposed modifications
The modified Snel model is implemented as described in Sect. 2.1 and tested against the same set of experiments as in Table 1. The results of the modified Snel model are shown in Figs. 7 and 8. In comparison to Figs. 2 and 3 it is noted that the shedding prediction at low reduced frequency is improved. For the higher reduced frequency the model also cap-tures the shedding slightly better, even though there is less shedding present. Furthermore, it can be seen that the frequency changes between both cases as desired. It is important to investigate the impact changes on different situations and aerofoils. Figure 9 displays the updated model in combination with the S809 aerofoil. The improved model still predicts slightly lower values than the experiments but the  overall trends are followed nicely, and shedding frequency matches the experiment well. Another goal of the proposed changes to the Snel model was to capture, predict, and match the vortex shedding and the shedding frequency of aerofoils in dynamic stall conditions. In order to check the validity of these changes in a quantitative way, a frequency domain analysis has been performed. The power spectral density (PSD) estimate is calculated using Welch's method. The Hamming window is set equal to the number of data points and the number of overlapped values to 50 % of the window length. The forcing frequency is removed from the plot as the shedding frequen-cies are higher and the forcing frequency will take up a large proportion of the PSD. The results for both aerofoils are shown in Figs. 10 and 11. From the figures it becomes clear that the Snel model captures, for both aerofoils at different mean angles of attack and different reduced frequencies, the self-induced shedding frequencies fairly correct. However, as shown in Figs. 7-9, there is still room for improvement here. All predicted shedding frequencies match frequencies observed in the measurements, whereas the intensity is not always correct. Care must be taken with the higher frequencies as the OSU measurement has a relatively low sampling  frequency and might therefore not fully capture some higherfrequency dynamics.
To quantify the improvements, the effects of all previous changes on the overall absolute error of the new model are shown in Fig. 12 together with the overall error of the initial model of Sect. 3. It is seen that the improved model outperforms the initial model in almost every case with a single exception. The initial model already gave very accurate results for that case, and the increase in error is very small compared to the reduction achieved in all other cases. It is also noted that the overall prediction of the shedding phenomena has been improved. Hence it can be concluded that the model changes developed and presented in this paper improve the performance of the Snel second-order model.

The updated Snel model performance with New MEXICO data
Model Experiments In Controlled Conditions (MEXICO) was a project in which an instrumented, three-bladed turbine of 4.5 m rotor diameter was tested (Schepers et al., 2012). MEXICO was carried out in the Large Low-Speed Facility (LLF) of the German-Dutch Wind Tunnels (DNW) (Schepers et al., 2012). The blades were fitted with pres-  . This data set is particularly valuable for the testing of dynamic stall models as it represents standard load cases set out by the IEC (DNV GL, 2016). For a complete overview of the MEXICO and New MEXICO projects, reference is made to Schepers et al. (2012). Khan (2018) conducted similar research with different dynamic stall models, and their study will be used as a baseline in this paper. The improved Snel model from this paper will be tested against the New MEX-ICO data sets of Table 3 ( Figure 15 shows the normal force distribution for the blade at azimuth 120 • for the yawed case with 30 • wind direction and 90 • pitch. Interesting to notice for this case is that the angle of attack will be negative. The Snel model does not show any unsteadiness in the normal force distribution. Further research is needed to explain and improve the Snel model for negative angles of attack. A single reason for the higher-frequency prediction has not been found. Several possibilities are suggested. First, the modal damping is not specified in the New MEXICO turbine data and was therefore assumed to be 0.5 %. Second, the im-  (Khan, 2018). Research by Skrypinski et al. (2014) and Zou et al. (2015) shows that for aerofoils with an effective angle of attack of around 90 • the Strouhal number should be around 0.10-0.13 Hz. This ought to be implemented in the model and tested in further research. Third, wind tunnel effects are not modelled in Bladed. The wind field in Bladed has zero turbulence and wind tunnel effects. More research on the Snel model in Bladed with actual turbine data is advised.

Conclusion and recommendations
The Snel model has been validated with OSU experimental data, and following this validation propositions for improvements have been made. The improvements to the model have been tested and led to a reduction in the overall error between the Snel model and the OSU experimental data. Furthermore, an improvement in the prediction of both the amplitude and frequency of vibrations in the measurements has been accomplished. The improved model is implemented in the Bladed turbine design software and tested against the New MEXICO experimental data. Prediction of normal force distributions along the blade seems to match earlier implementations in other turbine design codes while the mean value of the normal force is not correct. This is a major area for further research and improvement. The Snel model predicts the amplitude of the normal force vibrations well while the predicted frequency is higher than in the experiment. A single reason for this has not been found, and therefore further research into the Snel model is advised. Truong (2017) proposed a similar modification of his original model as Snel (1997) and therefore also similar to the model adapta-tions presented in this paper. The main difference lies in the incorporation of the steady lift data as described in Sect. 2.
The authors highly recommend a study comparing both the model described in Truong (2017) and the model presented in this paper. The model as proposed is formulated on the basis of variations in angles of attack. It is recommended for further research to delve into the possibility to adept this model to dynamic variation in inflow velocity for the case of rotating and vibrating blades. The proposed Snel second-order dynamic stall model might become a valuable addition to the modelling of dynamic behaviour in stall conditions. As the conditions tested in this paper are often design driving, the Snel model looks promising for more accurate prediction of design-driving (fatigue and extreme) loads and more cost-efficient wind turbine designs.
Code and data availability. The paper uses the publicly available OSU oscillating aerofoil experiment data. The implementation of the Snel model in this paper (and the final improved model) in the MATLAB environment is available and can be requested from the corresponding author.
Author contributions. GS has been the thesis supervisor within the MSc programme. He has been involved in writing (reviewing and editing) of the final thesis and the corresponding paper. MK was thesis supervisor from DNV GL and contributed to the implementation of the model in the MATLAB environment. Furthermore, MK has been a part of the investigation, validation, and improvement of the model as well as implementation of the model in the Bladed turbine design software. Finally, he has reviewed and edited the final thesis and corresponding paper.
Competing interests. The authors declare that they have no conflict of interest.

Special issue statement.
This article is part of the special issue "Wind Energy Science Conference 2019". It is a result of the Wind Energy Science Conference 2019, Cork, Ireland, 17-20 June 2019.