Articles | Volume 5, issue 3
Wind Energ. Sci., 5, 1037–1058, 2020
https://doi.org/10.5194/wes-5-1037-2020
Wind Energ. Sci., 5, 1037–1058, 2020
https://doi.org/10.5194/wes-5-1037-2020

Research article 20 Aug 2020

Research article | 20 Aug 2020

An improved second-order dynamic stall model for wind turbine airfoils

An improved second-order dynamic stall model for wind turbine airfoils
Galih Bangga1, Thorsten Lutz1, and Matthias Arnold2 Galih Bangga et al.
  • 1Institute of Aerodynamics and Gas Dynamics (IAG), University of Stuttgart, 70569 Stuttgart, Germany
  • 2Wobben Research and Development GmbH, 26607 Aurich, Germany

Correspondence: Galih Bangga (bangga@iag.uni-stuttgart.de)

Abstract

Robust and accurate dynamic stall modeling remains one of the most difficult tasks in wind turbine load calculations despite its long research effort in the past. In the present paper, a new second-order dynamic stall model is developed with the main aim to model the higher harmonics of the vortex shedding while retaining its robustness for various flow conditions and airfoils. Comprehensive investigations and tests are performed at various flow conditions. The occurring physical characteristics for each case are discussed and evaluated in the present studies. The improved model is also tested on four different airfoils with different relative thicknesses. The validation against measurement data demonstrates that the improved model is able to reproduce the dynamic polar accurately without airfoil-specific parameter calibration for each investigated flow condition and airfoil. This can deliver further benefits to industrial applications where experimental/reference data for calibrating the model are not always available.

1 Introduction

An accurate prediction of wind turbine blade loads is influenced by many parameters including 3D and unsteady effects. The first mainly occurs in the root and tip areas of the blade due to radial flow and induced velocity influences, respectively (Bangga2018). The latter can occur due to variation in the inflow conditions caused by yaw misalignment, wind turbulence, shear and gusts, tower shadow, and aeroelastic effects of the blade. The abovementioned phenomena may result in dynamic stall (DS). Experimental studies (Martin et al.1974; Carr et al.1977; McAlister et al.1978) showed that the aerodynamic forces can differ significantly in comparison to the static condition. DS is often initiated by the generation of a leading edge vortex (LEV), which increases the positive circulation effect on the airfoil suction side, causing delayed stall (Bangga2019). This intense leading edge vortex is convected downstream along the airfoil towards the trailing edge. At the same time, the lift force increases significantly, and the pitching moment becomes more negative compared to the static values. A significant drag increase is observed at large angles of attack. An example is shown in Fig. 1. Afterwards, a trailing edge vortex (TEV) with opposite rotational direction than LEV is formed, which pushes the leading edge vortex towards the wake area. This onset may result in a significant drop of the lift coefficient (CL) and can be dangerous for the blade structure itself, although dynamic stall also generally enhances aerodynamic damping.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f01

Figure 1Typical dynamic stall behavior of S801 airfoil. Data obtained from Ramsay et al. (1996).

Download

To model the behavior of the airfoil under these situations, semiempirical models can be used. The models are known to produce reasonable results without any notable increase in computational effort. Despite that, these models usually cannot reproduce higher harmonics of the load fluctuations. Furthermore, the applied constants shall be adjusted according to the flow conditions and airfoils. Leishman and Beddoes (LB) (Leishman and Beddoes1989) have developed a model for dynamic stall combining the flow delay effects of attached flow with an approximate representation of the development and effect of separation (Larsen et al.2007). This model was developed for helicopter applications and therefore includes a fairly elaborate representation of the nonstationary attached flow depending on the Mach number and a rather complex structure of the equations representing the time delays (Larsen et al.2007). Hansen et al. (2004) simplified the model for wind turbine applications by removing the consideration of compressibility effects and the leading edge separation. The latter was argued because the relative thickness of wind turbine airfoil is typically no less than 15 %. This model was called the Risø model in Larsen et al. (2007). Examples of the other models are given by Øye (Øye1991), Tran and Petot (Tran and Petot1980) (ONERA model), and Tarzanin (Boeing Vertol model) (Tarzanin1972). To better model the vortex shedding characteristics at large angles of attack, second-order dynamic stall models were introduced. An example of this model was given by Snel (1997) which makes use of the difference between the inviscid and viscous static polar data as a main forcing term for the dynamic polar reconstruction, in contrast to the LB model that uses the changes of the angle of attack over the time. An improved version of the Snel model was proposed recently by Adema et al. (2020) to cover for the increased shedding effects in the downstroke phase. All abovementioned models employ the static polar data and dynamic flow parameters as the input needed for the dynamic polar reconstruction. Then, the models compute the dynamic force difference required for the reconstruction process.

Although many studies have been dedicated to dynamic stall modeling (Gupta and Leishman2006; Larsen et al.2007; Adema et al.2020; Elgammi and Sant2016; Wang and Zhao2015; Sheng et al.2006; Galbraith2007; Sheng et al.2008), engineering calculations in the industry are still relying on the very basic classical dynamic stall models such as the Leishman–Beddoes and Snel models. The reason is the simplicity to tune in the models for different airfoils and for different flow conditions. Therefore, one major key for a model to be used in industrial applications is robustness of the model itself; i.e., the model is easy to apply with a small number of well-defined user parameters. One of the purposes of this paper is to document widely used dynamic stall models in research and industries. These include the first-order LB model and the second-order Snel model. A very recently improved Snel model according to Adema et al. (2020) will also be evaluated. The mathematical formulations of these models will be presented in this report. Weaknesses of existing dynamic stall modeling shall be identified, and possible corrections to those limitations will be described. Finally, a new second-order dynamic stall modeling will be proposed that is able to model not only the second-order lift and drag forces, but also the pitching moment along with calculation examples in comparison to experimental data for different airfoils and flow conditions.

The paper is organized as follows. Section 2 describes the mathematical formulation of existing dynamic stall models and the new model developed in this work. Then, in Sect. 3 assessments are carried out on the performance of the second-order dynamic stall models in comparison with measurement data. The new model is further tested at various flow conditions and to examine its robustness on four different airfoils without further calibrating the constants. Finally, all results will be concluded in Sect. 4.

2 Mathematical formulations

In this section the mathematical formulations of each model are described in detail. The reasons are mainly to provide information on how each model was employed and to gain deeper insights for further developing the new model. Note that each existing model was developed by different authors; thus different symbols and formulation methods were adopted in those publications (Beddoes1982; Leishman1988; Leishman and Beddoes1989; Snel1997; Adema et al.2020). In this paper, all models are described in a consistent way for clarity and for an easier interpretation and implementation process.

2.1 Leishman–Beddoes model

The original Leishman–Beddoes model is composed of three main contributions representing various flow regimes: (1) unsteady attached flow, (2) unsteady separated flow and (3) dynamic stall. The present section will elaborate the mathematical description and its physical interpretation of each module. Figure 2 illustrates several main parameters needed for modeling the dynamic stall characteristics.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f02

Figure 2Illustration of the main aerodynamic parameters needed for modeling the dynamic stall characteristics.

Download

2.1.1 Unsteady attached flow

In this module, the unsteady aerodynamic response of the loads is represented by the time delay effects. The indical formulas were constructed based on the work of Beddoes (1982) and have been refined by Leishman (1988). The loads are assumed to originate from two main sources: one for the initial noncirculatory loading from the piston theory and another for the circulatory loading which builds up quickly to the steady-state value (Leishman and Beddoes1989). In the formulation, the relative distance traveled by the airfoil in terms of semi-chords is represented by s=2Vt/c, which can also be used to describe the nondimensional time. Note that V, t and c are free-stream wind speed, time and chord length, respectively. For a continuously changing angle of attack αn, the effective angle of attack (αen) can be represented as

(1) α e n = α n - X n - Y n ,

where n is the current sample time. The last two terms describe the deficiency functions that are given by

(2)Xn=Xn-1exp-b1β2Δs+A1Δαnexp-b1β2Δs/2,(3)Yn=Yn-1exp-b2β2Δs+A2Δαnexp-b2β2Δs/2,

where

(4)Δαn=αn+1-αn,(5)Δs=sn-sn-1.

In these equations, b1, b2, A1 and A2 are constants. The variable β represents the compressibility effects and is formulated as β=1-M2. Because information about the previous cycle is needed in the formulations, initializations are required. The solution needs to develop for a certain time until convergence of the resulting unsteady loads is obtained.

The circulatory normal force due to an accumulating series of step inputs in angle of attack can be obtained using

(6) C N n C = d C N d α α e n - α 0 INV .

The variable α0INV is the angle of attack for zero inviscid normal force. The original formulation of the model disregarded the use of α0INV. However, this term is important when the airfoil has a finite camber. This has been pointed out as well by Hansen et al. (2004).

The noncirculatory (impulsive) normal force is obtained by

(7) C N n I = 4 K α T I M Δ α n Δ t - D n ,

where TI is given by TI=Mc/V. The deficiency function Dn is given by

(8) D n = D n - 1 exp - Δ t K α T I + Δ α n - Δ α n - 1 Δ t exp - Δ t 2 K α T I

and Δt=tn-tn-1.

The total normal force coefficient under attached flow conditions is given by the sum of circulatory and noncirculatory components as

(9) C N n P = C N n C + C N n I .

2.1.2 Unsteady separated flow

Leishman and Beddoes (1989) stated that the onset of leading edge separation is the most important aspect in dynamic stall modeling. The condition under which leading edge stall occurs is controlled by a critical leading edge pressure coefficient that is linked into the formulation by defining a lagged normal force coefficient CNnP1 as

(10) C N n P 1 = C N n P - D p n ,

where Dpn is given by

(11) D p n = D p n - 1 exp - Δ s T p + C N n P - C N n - 1 P exp - Δ s 2 T p .

It has been investigated by Leishman and Beddoes (1989) that the calibration time constant Tp is largely independent of the airfoil shape. The substitute value of the effective angle of attack incorporating the leading edge pressure lag response may be obtained using

(12) α f n = α 0 INV + C N n P 1 d C N / d α .

In most airfoil shapes, the progressive trailing edge separation causes loss of circulation and introduces nonlinear effects on the lift, drag and pitching moment, especially on cambered airfoils. This is even more important for wind turbine airfoils because the relative thickness is large. To derive a correlation between the normal force coefficient and the separation location (fn), the relation based on the flat plate from Kirchhoff and Helmholtz can be used, which reads

(13) C N n VISC = d C N d α 1 + f n 2 2 α n - α 0 VISC .

The location of the separation point is usually obtained by a curve-fitting procedure in literature. For example, Leishman and Beddoes (1989) proposed the following correlation:

(14) f n = 1 - 0.3 exp α n - α 1 S 1 ; α f n α 1 0.04 + 0.66 exp α 1 - α n S 2 ; α f n > α 1 .

The coefficients S1 and S2 define the static stall characteristic while α1 defines the static stall angle. The derivation was based on the NACA 0012, HH-02 and SC-1095 airfoils that have a single break point of the static lift force coefficient. Gupta and Leishman (2006) proposed the formulation for the S809 airfoil as

(15) f n = c 1 + a 1 exp S 1 α n ; α f n α 1 c 2 + a 2 exp S 2 α n ; α 1 < α f n < α 2 c 3 + a 3 exp S 3 α n : α f n α 2

that has two break points (α1 and α2) of the static lift force coefficient, where c1, c2, c3, a1, a2 and a3 are constants.

The additional effects of the unsteady boundary layer response may be represented by application of a first-order lag to the value of fn to produce the final value for the unsteady trailing edge separation point f2n (Leishman and Beddoes1989). This can be represented as

(16) f 2 n = f n - D f n ,

where Dfn is given by

(17) D f n = D f n - 1 exp - Δ s T f + f n - f n - 1 exp - Δ s 2 T f ,

and Tf is a constant. Then, the unsteady viscous normal force coefficient for each sample time can be obtained using

(18) C N n f = d C N d α 1 + f 2 n 2 2 α e n - α 0 VISC + C N n I .

The tangential component of the force can be obtained by (Leishman and Beddoes1989)

(19) C T n f = - η d C N d α α e n 2 f 2 n .

Note that positive CTnf is defined in the direction of the trailing edge while η is a constant.

According to Leishman and Beddoes (1989) and Gupta and Leishman (2006), a general expression for the pitching moment behavior cannot be obtained from Kirchhoff theory, and an alternative empirical relation must be formulated. Gupta and Leishman (2006) proposed the following formulation for the S809 airfoil

(20) C M f = C M 0 + K 0 + K 1 1 - f 2 n + K 2 sin π f 2 n m ; α n α 2 C M 0 + K 0 + K 3 exp K 4 f 2 n m ; α n > α 2 .

where CM0 defines the moment coefficient at zero normal force and K0 is the mean offset of the aerodynamic center from the quarter chord position. K1, K2, K3, K4 and m are constants.

2.1.3 Dynamic stall

The third part of the model describes the post-stall characteristics where the vortical disturbances near the leading edge become stronger. The effect of vortex shedding is given by defining the vortex lift as the difference between the linearized value of the unsteady circulatory normal force and the unsteady nonlinear normal force obtained from the Kirchhoff approximation, which reads

(21) C V n = C N n C 1 - K n ,

where Kn is given by

(22) K n = 1 4 1 + f 2 n 2 .

The normal force is allowed to decay, but it is updated with a new increment in the normal force based on prior forcing conditions, which can be defined as

(23) C N n V = C N n - 1 V exp - Δ s T v + C V n - C V n - 1 exp - Δ s 2 T v ; if 0 < τ v n < T vl C N n - 1 V exp - Δ s T v ; otherwise ,

where Tv and Tvl are the vortex decay and center of pressure travel time constants, respectively. The nondimensional vortex time is given by (Pereira et al.2011; Elgammi and Sant2016)

(24) τ v n = τ v n - 1 + 0.45 Δ t c V ; if C N n P 1 > C N CRIT 0 ; if C N n P 1 < C N CRIT and Δ α n > 0 ,

with CNCRIT being the inviscid critical static normal force, usually indicated by the break of the (viscous) moment polar at the critical angle of attack αnCRIT. This can be formulated as

(25) C N CRIT = d C N d α α n CRIT - α 0 INV .

The idealized variation in the center of pressure with the convection of the leading edge vortex can be modeled by

(26) C Pv n = K v 1 - cos π τ v T vl .

The dynamic moment coefficient can be formulated as

(27) C M n V = - C Pv n C N n V .

Therefore, the total dynamic loading on the airfoil from all modules can be written as

(28)CNnD=CNnf+CNnV,(29)CTnD=CTnf,(30)CMnD=CMnf+CMnV,

and by converting these forces into lift and drag, one obtains

(31)CLnD=CNnDcosαn-CTnDsinαn,(32)CDnD=CNnDsinαn+CTnDcosαn.

2.1.4 Note to present implementation

In Eqs. (14) and (15), a curve-fitting procedure is usually adopted in literature. In this sense, the parameters or even the formulation need to be adjusted when the airfoil is different. Therefore, in the present implementation, the separation point is derived directly from the static polar data using the inversion of Eq. (13) as

(33) f n = 2 C N n VISC d C N d α α f n - α 0 VISC - 1.0 2 .

The same approach was used for example by Hansen et al. (2004). This way, the user can avoid dealing with curve fitting adjustment (which requires changes on the constants for different airfoils and flow conditions) as long as the static polar data are available.

In the original formulation, the pitching moment is also obtained by a curve-fitting procedure in Eq. (20). Again, this kind of approach is not straightforward as the user needs to perform curve fitting of the polar data. In the present implementation, the moment coefficient is easily obtained from the static viscous polar data by interpolating the value at the effective angle of attack, incorporating the leading edge pressure time lag αfn, which reads

(34) C M n f = C M VISC α f n .

In this sense, the moment coefficient can be reconstructed easily without the need to adjust the parameters in advance, minimizing the user error.

Furthermore, to avoid discontinuity in the downstroke phase for Eq. (24), an additional condition is applied in the present implementation as

(35) τ v n = τ v n - 1 + 0.45 Δ t c V ; if C N n P 1 > C N CRIT 0 ; if C N n P 1 < C N CRIT and Δ α n 0 τ v n - 1 ; otherwise .

2.2 Snel second-order model

The history of the Snel second-order model (Snel1997) dates back to 1993 based on Truong's observation on dynamic lift coefficient characteristics (Truong1993). Truong proposed that the difference between the static and dynamic lift can be divided into two terms: the forcing frequency response and the higher-frequency dynamics of a self-excited nature. The total dynamic response of the airfoil is formulated as

(36)CLnD=CLnVISC+ΔCLnD1+ΔCLnD2,(37)CDnD=CDnVISC+,(38)CMnD=CMnVISC+,

with D1 and D2 being the first- and second-order corrections, respectively. The first correction is modeled using an ordinary differential equation (ODE) by applying a spring-damping like function as

(39) τ Δ C ˙ L n D 1 + K f 10 n Δ C L n D 1 = F 1 n .

The frequency of the first-order corrected lift follows the frequency of the forcing term F1. This term is based on the time derivative of the difference between the steady inviscid coefficient CLnINV and viscous lift coefficient CLnVISC of an airfoil (ΔCLnINV) as

(40)F1n=τΔC˙LnINV,(41)ΔCLnINV=CLnINV-CLnVISC=dCLdααn-α0INV-CLnVISC,

with n and dCL∕dα as the current sample time and inviscid lift gradient, respectively. The time constant τ in the above equation represents the time required for the flow to travel half a chord distance as

(42) τ = c 2 V .

The “stiffness” coefficient of the first-order term Kf10n can be expressed as

(43) K f 10 n = 1 + 0.5 Δ C L n INV 8 1 + 60 τ α n ˙ ; if α n ˙ C L n INV 0 1 + 0.5 Δ C L n INV 8 1 + 80 τ α n ˙ ; if α n ˙ C L n INV > 0 .

As shown in Faber (2018), the above equation becomes numerically unstable if αn˙ is large (increasing reduced frequency above 0.1) for αn˙CLnINV0. The reason is that the denominator goes to zero and then negative, causing numerical integration instability. Thus, based on pure intuition the denominator value was set to a minimum of 2.0 in Faber (2018). In the present implementation, a similar approach is adopted but the limit differs. Instead, the minimum denominator value is limited to 1×10-5, because it yields more physical results for several cases tested by the authors.

To incorporate the higher-order frequency dynamics, a second-order ODE is used to describe the second-order correction term. The general form may be written as

(44) τ 2 Δ C ¨ L n D 2 + K f 21 n Δ C ˙ L n D 2 + K f 20 n Δ C L n D 2 = F 2 n

similar to the first-order correction. The frequency of the higher-order dynamics is determined by the forcing term F2n, defined as

(45) F 2 n = 0.1 k s - 0.15 Δ C L n INV + 0.05 Δ C ˙ L n INV .

It is noted that the value 0.1 as a constant was chosen according to Adema et al. (2020). This is not a fixed value and can be adjusted based on the evaluated cases as seen in literature (Adema et al.2020; Snel1997; Holierhoek et al.2013; Faber2018; Khan2018). Variable ks is a constant with a typical value of 0.2. The spring coefficient is given by

(46) K f 20 n = k s 2 1 + 3 Δ C L n D 2 2 1 + 3 α n ˙ 2

and the damping coefficient as

(47) K f 21 n = 60 τ k s - 0.01 Δ C L n INV - 0.5 + 2 Δ C L n D 2 2 ; if α n ˙ > 0 2 τ k s ; if α n ˙ 0 .

2.3 Adema–Snel second-order model

The recently developed model of Adema et al. (2020) improves the original Snel model (Snel1997) in several aspects. Instead of using the lift coefficient (CL), the normal force coefficient (CN) is used, similar to the LB model (Leishman and Beddoes1989). The total dynamic response of the airfoil is formulated as follows.

(48)CNnD=CNnVISC+ΔCNnD1+ΔCNnD2(49)CTnD=CTnVISC+(50)CMnD=CMnVISC+

The model introduces some modifications of the original model in terms of (1) projected ks, (2) the first-order coefficient and (3) the second-order coefficient. The mathematical formulation of the first-order term of the model is listed as

(51)τΔC˙NnD1+Kf10nΔCNnD1=F1n,(52)F1n=τΔC˙NnINV,(53)ΔCNnINV=CNnINV-CNnVISC=dCNdααn-α0INV-CNnVISC,(54)Kf10n=1+0.2ΔCNnINV81+60ταn˙;ifαn˙CNnINV01+0.2ΔCNnINV81+80ταn˙;ifαn˙CNnINV>0,

and for the second-order correction term as

(55)τ2ΔC¨NnD2+Kf21nΔC˙NnD2+Kf20nΔCNnD2=F2n,(56)F2n=0.01ks-0.04ΔCNnINV+1.5τΔC˙NnINV,(57)Kf20n=10kssinαn21+3ΔCNnD221+2802τ2αn˙2,(58)Kf21n=60τks-0.01ΔCNnINV-0.5+2ΔCNnD22;ifαn˙>060τks-0.01ΔCNnINV-0.5+14ΔCNnD22;ifαn˙0.

One may notice that Eq. (56) contains τ in the second term of the right-hand side (RHS). This is intended to remove the dependency of the model on the velocity as the input parameter. The other main difference with the original model is also observed in Eq. (57), where ks is projected by sin αn. At last, the downstroke motion of the second-order term of Eq. (58) is modified to enable vortex shedding effects.

To sum up the characteristics of the above-discussed state-of-the-art dynamic stall models, Table 1 lists the properties of each model and in which aspects the model can be improved further.

Table 1Properties of the discussed dynamic stall models.

Download Print Version | Download XLSX

2.4 New second-order IAG model

The proposed IAG model is developed based on knowledge gained from four different models: Leishman–Beddoes, Snel, Adema–Snel and ONERA (Tran and Petot1980; Dat and Tran1981; Petot1989) models with modifications. Similar to the modern models like those from Snel (and ONERA) and its derivatives, the present model is constructed by two main terms: the first-order and second-order corrections. The total dynamic response of the airfoil is formulated as

(59)CLnD=CLnD1+ΔCLnD2,(60)CDnD=CDnD1+ΔCDnD2,(61)CMnD=CMnD1+ΔCMnD2,

with D1 and D2 being the first- and second-order corrections, respectively. Below the description of the modifications made for the new model will be discussed in detail.

2.4.1 First-order correction

Based on the Hopf bifurcation model of Truong (1993) that used the LB model as the starting point of the first-order correction, the present model operates similarly. Despite that, the LB model is not transferred into the state-space formulation, but it is retained as the indical formulation. The model applies the superposition of the solution using a finite-difference approximation to Duhamel's integral to construct the cumulative effect on an arbitrary time history of angle of attack. The LB model described in Sect. 2.1.1 to 2.1.3 will be used with the following modifications.

In the above LB model, predictions for drag are not accurate as will be shown in Sect. 3.1. This inaccuracy lies in the determination of η in Eq. (19) for the tangential force component because drag is more sensitive to tangential force than the lift force is. Also, to maintain simplicity, parameter η is removed, and the tangential force is obtained from the static data at the time-lagged angle of attack αfn by

(62) C T n f = C T VISC α f n .

If one uses this formulation directly, at some point drag still becomes lower than the static drag value by a significant amount. By evaluating the experimental data for several airfoils and various flow conditions, this is not physical at small angles of attack, especially in the downstroke regime, where it usually just returns to the static value. In fact, those experimental data infer that strong drag hysteresis occurs only at high angles of attack beyond stall. Similarly, in the upstroke regime the drag value increases only slightly (approximately only 20 %). In Fig. 3, one can see that drag hysteresis occurs when ζnζv, with

(63) ζ n = 1 π d C N d α 1 + f n 2 2 ,

and ζv is a constant with a value of 0.76. Based on these observations, a simple drag limiting factor is adopted when ζn≥0.76 as

(64) C D n D = 1.2 C D n VISC ; if C D n D > 1.2 C D n VISC and C N n P - C N n - 1 P 0.0 C D n VISC ; if C N n P - C N n - 1 P < 0.0 C D n D ; otherwise .

Note that for the purpose of numerical implementation, it is always recommended in practice to adopt relaxation to avoid any discontinuity which may present in the above formulation. Furthermore, the value of ζv may be chosen differently for different airfoils depending on the vortex lift influence on drag. Further studies on this aspect are encouraged. The effects of these modifications are displayed in Fig. 4.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f03

Figure 3Relation between drag hysteresis in the stall regime with weighted separation parameter ζ for four airfoils. From (a) to (d): S801 (13.5 %), NACA4415 (15 %), S809 (21 %) and S814 (24 %).

Download

It will also be shown in Sect. 3.1 that predicting moment coefficient directly from the static polar data by means of the time-lagged angle of attack has its drawback in the correct damping effect calculation. One may obtain better results by using the “fitting function” as in Eq. (20) instead of using Eq. (34). However, this limits the usability for different airfoils, since the fitting has to be done again for each individual airfoil. For wind turbine simulations, this is fairly impractical because a wind turbine blade is usually constructed by several different airfoils, not to mention the interpolated shapes in between each airfoil position. To solve for this issue, a relatively simple approach is introduced by applying a time delay to the circulatory moment response as

(65) C M n C = C M n - 1 C exp - Δ s T M U - C P f n C V n - C V n - 1 exp - Δ s 2 T M U ; if τ v n < T vl and Δ α n 0 C M n - 1 C exp - Δ s T M D - C P f n C V n - C V n - 1 exp - Δ s 2 T M D ; if Δ α n < 0 C M n - 1 C ; otherwise

where

(66) C P f n = K f C d C N d α α n CRIT - α 0 INV ,

with KfC, TMU and TMD being constants relatively insensitive to airfoils. Furthermore, the second condition of Eq. (35) is modified to avoid discontinuity which occurs at a large reduced frequency (e.g., k=0.2) for any LB-based models without recalibration of the time constant as

(67) τ v n = τ v n - 1 + 0.45 Δ t c V ; if C N n P 1 > C N CRIT τ v n - 1 exp ( - Δ s ) ; if C N n P 1 < C N CRIT and Δ α n 0 τ v n - 1 ; otherwise .

The effects of these modifications are displayed in Fig. 5.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f04

Figure 4Drag reconstruction in comparison with the experimental data for the S801 airfoil (Ramsay et al.1996) applying (a) Eq. (19), (b) Eq. (62), and (c) Eqs. (62) and (64).

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f05

Figure 5Moment reconstruction in comparison with the experimental data for the S801 airfoil (Ramsay et al.1996) applying (a) Eq. (34) and (b) Eq. (65).

Download

The total first-order dynamic response of the airfoil is formulated as

(68)CNnD1=CNnf+CNnV,(69)CTnD1=CTnf,(70)CMnD1=CMnf+CMnV+CMnC,

where the definition and description of each variable were given in Sect. 2.1 for the LB model. Thus the first-order lift and drag responses can be obtained by

(71)CLnD1=CNnD1cosαn-CTnD1sinαn,(72)CDnD1=CNnD1sinαn+CTnD1cosαn.

2.4.2 Second-order correction

The second-order correction takes the form of the non-linear ordinary differential equation according to the second-order correction of the Snel (Snel1997) and Adema–Snel (Adema et al.2020) models. Generally, the basis of implementation of the present model mostly uses the Adema–Snel (Adema et al.2020) model where the normal force is used instead of the lift force as for the original Snel model (Snel1997) as

(73) Δ C ¨ N n D 2 + K f 21 n Δ C ˙ N n D 2 + K f 20 n Δ C N n D 2 = F 2 n .

This way, shedding effects apply not only to the lift force but also to the drag force. Note that τ is not present in Eq. (73) in contrast to the original formulation in Eqs. (44) and (55). The equation is changed because the time derivatives in the above equation are no longer with respect to time but to s=2Vt/c, similar to the ONERA model (Tran and Petot1980; Dat and Tran1981; Petot1989). This is done to nondimensionalize the equations.

In Eq. (57), the ks was projected as a function of the angle of attack by sin αn. This modification causes problems when the hysteresis effect takes place in both positive and negative angles because Eq. (57) will be zero and then negative, causing instability of the ODE. Thus, the original form of the Snel model (Snel1997) is retained in the present model, but the constant is modified as

(74) K f 20 n = 20 k s 2 1 + 3 Δ C N n D 2 2 1 + 3 α n ˙ 2 .

The idea for the downstroke damping as in Eq. (58) is adopted in the present model; the following form and constants are used:

(75) K f 21 n = 150 k s - 0.01 Δ C N n INV - 0.5 + 2 Δ C N n D 2 2 ; if α n ˙ > 0 30 k s - 0.01 Δ C N n INV - 0.5 + 14 Δ C N n D 2 2 ; if α n ˙ 0 and α n α n CRIT 0.2 k s ; if α n ˙ 0 and α n < α n CRIT .

Note again that τ is not present in the above equation. The original formulation in Eq. (58) replaces the damped oscillator when αn˙0 for a self-excited oscillator of Van der Pol type with more damping. This is in contrast with the implementation done in Truong (2016), where the self-excited oscillator is only replaced by the damped oscillator, when the flow is reattached on the return cycle. Under such circumstances, the oscillatory behavior still subsists in the return cycle, albeit with smaller amplitude. To accommodate this aspect, the last term of Eq. (75) is applied when the angle is smaller than αnCRIT. As for the forcing term, the original form of the Snel model (Snel1997) is adopted as

(76) F 2 n = 0.5 k s - 0.15 Δ C N n INV + 0.05 Δ C ˙ N n INV .

To facilitate the inclusion of the higher harmonic effects for the pitching moment, the idealized center of pressure obtained in the first-order correction given in Eq. (26) is passed into the second-order model. Thus, the dynamic moment reaction takes the form

(77) Δ C M n D 2 = - C Pv n Δ C N n D 2 .

Regarding the tangential force, a similar assumption is made as in Eq. (49) where the influence of ΔCTnD2 is neglected in the formulation, but ΔCTnD1 is considered in the first-order term correction described above. Finally, the second-order term of the lift (ΔCLnD2) and drag (ΔCDnD2) can be calculated easily. The effects of the second-order term are displayed in Fig. 6.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f06

Figure 6Airfoil response reconstruction in comparison with the experimental data for the S801 airfoil (Ramsay et al.1996) applying only the first-order correction and with inclusion of the second-order term. (a) Lift, (b) drag and (c) pitching moment.

Download

2.5 Constants applied for the investigated dynamic stall models

The following constants are applied in the implemented dynamic stall models. These values are kept constant throughout the paper. The constants for the Leishman–Beddoes model and for the proposed IAG model are given in Tables 2 and 3, respectively. For any model developed based on the Leishman–Beddoes type, the critical angle of attack plays a major role. This can be obtained as the angle where the viscous pitching moment breaks or when the drag increases significantly (or the stall angle). The applied critical angles are given in Table 4. The validation is done by comparing the calculations with experimental data obtained at The Ohio State University for the S801 airfoil (13.5 % relative thickness) (Ramsay et al.1996), NACA4415 airfoil (15 % relative thickness) (Hoffman et al.1996), S809 airfoil (21 % relative thickness) (Ramsay et al.1995) and S814 airfoil (24 % relative thickness) (Janiszewska et al.1996). All selected test cases for the airfoils are employed with a leading edge grit (turbulator) to enable the “soiled” effects on a wind turbine blade at a Reynolds number of around 750 K. Note that these polar data are different from those used for example by Sheng et al. (2010) where the natural transition cases were taken. Therefore, the critical angles of attack are also different. The results of the studies are presented in the following sections.

Table 2Constants applied for the Leishman–Beddoes model.

Download Print Version | Download XLSX

Table 3Constants applied for the IAG model.

Download Print Version | Download XLSX

Table 4Critical angle of attack (αnCRIT) applied for the Leishman–Beddoes and IAG models.

Download Print Version | Download XLSX

3 Results and discussion

The three state-of-the-art dynamic stall models reviewed above (Leishman–Beddoes, Snel, Adema–Snel) have been used as a basis for examining the dynamic loads of four different pitching airfoils at various flow conditions. Experience gained from those models is used to formulate a new second-order dynamic stall model, namely the IAG model, by evaluating the weakness and strength of each model. The presented second-order models need to solve a set of differential equations. For this purpose, the Euler–Heun forward integration method is used.

3.1 Comparison against experimental data

This section compares the predicted dynamic forces and the measurement data. For a fair comparison, all models are assessed with the same time step size of Δt=T/1440, with T being the pitching period. The evaluations are performed on the S801 airfoil at k=0.073. The comparison of each model is shown in Figs. 7 to 9 for the Snel, Adema and IAG models, respectively. Note that the constants of the other existing dynamic stall models are taken directly from literature without further calibration for the S801 airfoil. Therefore, it is already expected that their performance will not be optimal. The main purpose of the comparison is not to study their accuracy but to analyze the robustness of each model for a different airfoil without tuning the constants. On the other hand, the constants for the IAG model in Table 3 were obtained using the S801 airfoil. To enable a fair assessment of the model robustness, the IAG model will also be used to reconstruct the dynamic polar data of four airfoils with different relative thicknesses without changing the constants in Sect. 3.6.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f07

Figure 7Dynamic force reconstruction using the Snel model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440. S801 airfoil, k=0.073, α=20 and Δα=10.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f08

Figure 8Dynamic force reconstruction using the Adema model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440. S801 airfoil, k=0.073, α=20 and Δα=10.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f09

Figure 9Dynamic force reconstruction using the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440. S801 airfoil, k=0.073, α=20 and Δα=10.

Download

The original Snel models cannot predict the drag and moment coefficients in the original formulations. Thus, only the static polar data are shown. The Snel model actually shows an acceptable accuracy even though the constants are taken as found in literature. The higher harmonic effects are unfortunately not captured by this model. This is further refined by the Adema model which was developed as an improvement for the Snel model. The model performs fairly well for the lift and drag predictions, though the drag value at small angles of attack is a bit off. The pitching moment prediction is also not included in its formulations. These disadvantages are better treated in the proposed IAG model. Not only the prediction of the lift coefficient but also the accuracy of drag prediction are improved significantly. The modifications described in Sect. 2.4 result in the improvement at low and high angle-of-attack regimes. The model is also able to reconstruct the pitching moment polar accurately, which is important for aeroelastic calculations of wind turbine blades.

For the following sections, the proposed IAG model will be tested under various flow conditions and for several airfoils at various relative thicknesses in comparison with measurement data. Note that these calculations are performed without changing the constants to assess the robustness of the model at different flow conditions.

3.2 Effects of time signal deviation

The actual pitching motion within The Ohio State University (OSU) measurement differs slightly from the intended motion. The actual time series of the angle of attack is included in the experimental data (Ramsay et al.1996; Hoffman et al.1996; Ramsay et al.1995; Janiszewska et al.1996). To assess the effects of this time signal deviation on the aerodynamic response, the calculations using these time signal data were performed applying the IAG model. The results are compared with the experimental data and the results of the IAG model presented in Sect. 3.1. Note that these time signal data are fairly coarse and can cause problems for second-order dynamic stall models because the gradient of α change can be extremely large. To cover for this issue, the time signal is interpolated in between each available point using a third-order cubic-spline interpolation. The time signals are discretized by Δt=T/1440 over a single pitching period. The first period of oscillation is used for initialization of the time integration; thus a constant angle of attack is applied as shown in Fig. 10.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f10

Figure 10Comparison of the time series of the idealized sinusoidal angle of attack to the exact signals in the experimental campaign for the S801 airfoil, k=0.073, α=20 and Δα=10.

Download

Figure 11 presents the influence of the time signal variation on the aerodynamic performance in terms of CL, CD and CM. TS labels the exact time signals in the experimental campaign. Although the time signal difference has almost no influence on the global prediction characteristics, some deviation from the idealized sinusoidal motion can be noticed clearly. For example, one can see the increased lift buildup in the upstroke regime before stall and the location of the lift stall. Some deviations in the drag and pitching moment coefficients are observed as well. For the rest of the paper, the prediction procedure using the actual time signal from the experimental data is used for the best consistency with the experimental campaign.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f11

Figure 11Dynamic force reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign. TS labels the exact time signals in the experimental campaign. S801 airfoil, k=0.073, α=20 and Δα=10.

Download

3.3 Performance of the model for different mean angles of incidence

In this section, the effects of the mean angle of attack are evaluated. Three different angles of attack at the same inflow conditions are selected for this purpose. These are α=8, 14 and 20. Note that these mean angles of attack are only approximations since the actual time signal data from the experimental campaign are used. The selected mean angles represent the regime of attached flow, partly separated flow and fully separated flow conditions. These are helpful to assess the model performance under various flow situations.

Figure 12 presents the results for the lift coefficient under these three investigated mean angles of attack. The model performs very well for these different cases. The maximum lift is a bit overestimated in the model for the lowest α, but in general all unsteady lift characteristics in the measurement data are reproduced in a sound agreement with the experimental data. A similar behavior is shown for the drag prediction depicted in Fig. 13. The proposed model captures the increased drag effect and its shedding characteristics well. The simple modifications applied in Sect. 2.4 result in a good prediction of the drag coefficient compared with the experimental data. In Fig. 14, the prediction for pitching moment is shown. Here the predicted moment coefficient is in a good agreement with the measured values.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f12

Figure 12Lift reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign at various α values. From (a) to (c): α=8, α=14 and α=20. S801 airfoil, k=0.073 and Δα=10.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f13

Figure 13Drag reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign at various α values. From (a) to (c): α=8, α=14 and α=20. S801 airfoil, k=0.073 and Δα=10.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f14

Figure 14Pitching moment reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign at various α values. From (a) to (c): α=8, α=14 and α=20. S801 airfoil, k=0.073 and Δα=10.

Download

3.4 Performance of the model for different reduced frequencies

The effects of pitching frequency on the aerodynamic response will be discussed in this section. Three different reduced frequencies are examined, namely k=0.036, 0.073 and 0.111. The stall regime is shown here, where the prediction is the most challenging. The actual time signals as of the measurement campaign are used, following the procedure described in Sect. 3.2.

Figure 15 displays the results for the dynamic lift coefficient response. The lowest reduced frequency of 0.036 is dominated by the viscous effects. It represents the case where the “delayed” angle-of-attack response is the weakest. It can be seen that the maximum attained lift coefficient increases with increasing k. The gradient of the lift polar in the upstroke and downstroke phases also increases as well. These characteristics are present in both experimental data and predictions delivered by the IAG model. A similar behavior is also displayed in drag and pitching moments in Figs. 16 and 17, respectively. It is obvious that stall occurs much earlier for a smaller k value. One can see that the maximum amplitude of all three force components increases with increasing k. This can be dangerous for the structural stability, since the amplitude determines the fatigue loads.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f15

Figure 15Lift reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign at various k values. From (a) to (c): k=0.036, k=0.073 and k=0.111. S801 airfoil, α=20 and Δα=10.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f16

Figure 16Drag reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign at various k values. From (a) to (c): k=0.036, k=0.073 and k=0.111. S801 airfoil, α=20 and Δα=10.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f17

Figure 17Pitching moment reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign at various k values. From (a) to (c): k=0.036, k=0.073 and k=0.111. S801 airfoil, α=20 and Δα=10.

Download

To better investigate the effects of k, the IAG model is used to reconstruct the dynamic polar data at various k values by applying an idealized sinusoidal motion as presented in Fig. 18. Only the last DS cycle is shown for clarity of the observation. While the maximum amplitude of all three force components at low-frequency domains increases with increasing k (blue and green markers), the amplitudes for all three forces at high-frequency domains show different characteristics as shown in the Fourier transformation in Fig. 19, albeit with much smaller values. The higher harmonic amplitudes are attributed to flow separation effects, while for low-frequency domains are driven by the pitching motion (i.e., external unsteadiness or inflow).

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f18

Figure 18Effects of k on the aerodynamic response by the IAG model for Δt=T/1440. S801 airfoil, α=20 and Δα=10. (a–c) Polar; (d–f) time series.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f19

Figure 19Fourier transformation of the predicted forces presented in Fig. 18. fs=f/f0 with f0 being the pitching frequency.

Download

3.5 Performance of the model for different pitching amplitudes

In this section, the effects of pitching amplitude on the aerodynamic response of a pitching airfoil are investigated. The mean angle of attack is fixed at α=20. Note again that α is only an approximation because the actual time signal data from the measurement campaign are applied. This large mean angle of attack is purposely selected because the post-stall characteristic is of interest and is well known for its violent vibration, even for the static condition. The small amplitude in this case means that the whole pitch oscillation occurs within the stall regime.

Figures 20 to 22 display the dynamic force responses due to pitching motion of the airfoil predicted by the IAG model in comparison with the experimental data. The model accurately reconstructs the dynamic forces despite the predicted case being challenging within the post-stall regime. Interesting to note is that the small-pitching-amplitude case induces stronger shedding effects for lift than the larger-amplitude case. This can be explained as follows. As described by Leishman in his papers (Beddoes1982; Leishman1988; Leishman and Beddoes1989), the airfoil sees a lagged force response compared to the imposed disturbance. Therefore, in his model, a time-lagged angle of attack is introduced as the “effective” angle actually seen by the airfoil section. When the pitching motion takes place partly within the fully separated region (in the static case) and partly in the attached/partly separated flow region, the airfoil still sees the lower angle (where the flow is still attached) even though the pitching motion already reaches the post-stall regime. This effect stops/decreases when the effective angle is larger than the critical angle defined in Table 4. As the critical angle for the S801 airfoil is defined at 15.1, the lower-amplitude case is fully operating within the stall regime, where the attached flow effect is not present.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f20

Figure 20Lift reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign at various Δα values. (a) Δα=5.5; (b) Δα=10. S801 airfoil, k=0.073 and α=20.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f21

Figure 21Drag reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign at various Δα values. (a) Δα=5.5; (b) Δα=10. S801 airfoil, k=0.073 and α=20.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f22

Figure 22Pitching moment reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign at various Δα values. (a) Δα=5.5; (b) Δα=10. S801 airfoil, k=0.073 and α=20.

Download

3.6 Performance of the model for different airfoils

In this section, the performance and robustness of the proposed IAG model are assessed for airfoils with different relative thickness. All model constants in Table 3 remain the same for all calculations. The difference from one airfoil calculation to the other lies only in the critical angle of attack value as shown in Table 4. The value was obtained simply by looking at the static polar data where the viscous pitching moment breaks or when the drag increases significantly.

Despite the increased airfoil thickness from 13.5 % to 24 %, Figs. 23 to 25 demonstrate that the reconstructed dynamic forces are in a good agreement with the experimental data, not only for the general trend but also the higher harmonic effects. As is also the case for the Leishman–Beddoes model, it is important to select the appropriate value for the critical angle of attack. The simple approach used in the present paper has shown its usefulness and potentially reduces the complexity of parameter tuning for industrial applications. Elgammi and Sant (2016) for example defined two different critical angles of attack, one for CN and the other for CT that were shown to improve the prediction accuracy. Although their attempt might be beneficial, this is not followed in this work because one main aim of the studies is to reduce parameter tuning required for one case to the other cases.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f23

Figure 23Lift reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996; Hoffman et al.1996; Ramsay et al.1995; Janiszewska et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign for different airfoils. From (a) to (d): S801 (13.5 %), NACA4415 (15 %), S809 (21 %) and S814 (24 %). k=0.073, α=20 and Δα=10.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f24

Figure 24Drag reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996; Hoffman et al.1996; Ramsay et al.1995; Janiszewska et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign for different airfoils. From (a) to (d): S801 (13.5 %), NACA4415 (15 %), S809 (21 %) and S814 (24 %). k=0.073, α=20 and Δα=10.

Download

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f25

Figure 25Pitching moment reconstruction by the IAG model in comparison with the measurement data (Ramsay et al.1996; Hoffman et al.1996; Ramsay et al.1995; Janiszewska et al.1996) for Δt=T/1440 using the actual angle of attack in the experimental campaign for different airfoils. From (a) to (d): S801 (13.5 %), NACA4415 (15 %), S809 (21 %) and S814 (24 %). k=0.073, α=20 and Δα=10.

Download

3.7 Predictions of the center of pressure

To further complement the analyses conducted in Sect. 3.6, the location of the actual pressure center is calculated in this section as

(78) X p = - C M C L ,

which indicates the distance of the pressure point to the quarter chord position where CM is defined.

A correct location of the pressure point is important for determining the stability on aeroelastic simulations of wind turbine blades. The results of the calculations both for the experimental data and for the proposed IAG model are presented in Fig. 26 for all four investigated airfoils both as time series and as the polar plot. It can be seen clearly that the agreement between the experimental data and the present predictions is excellent for all investigated airfoils.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f26

Figure 26Center of pressure reconstruction in comparison with the measurement data by the IAG model for Δt=T/1440 using the actual angle of attack in the experimental campaign for different airfoils. From top to bottom: S801 (13.5 %), NACA4415 (15 %), S809 (21 %) and S814 (24 %). k=0.073, α=20 and Δα=10.

Download

3.8L2 norm of error analyses

Holierhoek et al. (2013) introduced a way of quantifying the absolute error between the experimental data and modeled lift coefficient. The general formulation reads

(79) L 2 ϕ = 1 N i N ϕ i mod - ϕ i exp 2 ,

with ϕ being the variable of interest, i the current sample and N the total number of samples. In their paper, however, only lift was considered. Here all three force components will be shown for four different airfoils. Figure 27 displays the quantified error for two different flow categories, attached and deep stall. The time series of the angle of attack was obtained from the measured data by applying a third-order cubic-spline interpolation in between each available point. One can see that generally the attached flow case is predicted very well, while the error increases as the flow condition becomes more complicated. Interestingly, especially for lift, it seems that the error decreases with increasing airfoil thickness. The reason for the larger error obtained for the thinner airfoil is attributed to the complex characteristics of the leading edge stall, causing severe load variations, especially with increasing angle of attack. Thus, it makes the prediction more challenging. Furthermore, the quantification of the error was also performed on two other dynamic stall models, Snel and Adema–Snel. The same approach for the angle of attack signal was applied. One can see that the IAG model shows its improved prediction in particular for the deep-stall case for all three force components.

https://wes.copernicus.org/articles/5/1037/2020/wes-5-1037-2020-f27

Figure 27Quantified L2 norm of error with respect to the measurement data for four airfoils. (a–c) Attached flow case (k=0.073, α=8, Δα=5.5); (d–f) deep-stall case (k=0.073, α=20, Δα=10).

Download

4 Conclusions

Comprehensive studies on the accuracy of several state-of-the-art dynamic models to predict the aerodynamic loads of a pitching airfoil have been conducted. From the studies, the strength and weaknesses of each model were highlighted. This information was then transferred to develop a new second-order dynamic stall model proposed in this paper. The new model improves the prediction for the aerodynamic forces and their higher-harmonic effects due to vortex shedding, developed for robustness to improve its usability in practical wind turbine calculations. Details on the model characteristics, modifications and treatment for numerical implementation were summarized in the present paper. The studies were conducted by examining the influence of the time step size, time signal deviation, mean angle of attack, reduced frequency, pitching amplitude and variation in the airfoil thickness. Several main conclusions can be drawn from the work.

  • The general characteristics of the polar data can be predicted by all investigated dynamic stall models. Despite that, only the Adema model and the present IAG model are able to demonstrate the higher harmonic effects among the three investigated models.

  • The exact time signal imposed based on the measurement campaign improves the prediction accuracy of the IAG model in comparison with the idealized sinusoidal motion.

  • The dynamic forces reconstructed by the IAG model are in a sound agreement with the experimental data under various flow conditions by variation in α, k and Δα and for four different airfoils by changing only the values of the critical angle of attack.

  • The amplitudes at low-frequency domains increase with increasing k and can be attributed to the effects of inflow/external unsteadiness. The amplitudes at high-frequency domains decrease with increasing k values, which are driven by flow separation effects.

  • When the airfoil operates at a high α within the stall regime, a small Δα leads to increased vibrations for lift. The opposite is true for the pitching moment.

5 Recommendations for future work

The present paper evaluates the newly developed IAG model under various flow conditions for four different airfoils. The following aspects are encouraged for future work.

  • In the present studies, the assessment was mainly carried out for the S801 airfoil having a relative thickness of 13.5 %. This airfoil is mainly characterized by leading edge separation, which is very challenging for validating the accuracy of a dynamic stall model. However, typical modern wind turbine blades usually employ airfoils with no less than 18 % relative thickness and a much higher Reynolds number. Therefore, future investigations shall be done for thicker airfoils at various flow conditions as well.

  • The above statement is also true for the current available experimental data. Therefore, experiments on dynamic stall for thick airfoils at a much higher Reynolds number are encouraged.

  • Three-dimensional effects (Himmelskamp or tip loss effects) for a rotating blade can alter the loads significantly even under a steady inflow condition. Further consideration and examination of the model under this condition shall be carried out.

  • Further tests and recalibration of the model for deep-stall conditions at extremely large angles of attack are encouraged, which can be relevant for a turbine in standstill.

Appendix A: List of symbols
Variables
s Nondimensional time (s)
V Incoming wind speed (m s−1)
t Time (s)
c Chord (m)
CN Normal force coefficient ()
CT Tangential force coefficient ()
CL Lift force coefficient ()
CD Drag force coefficient ()
CM Pitching moment coefficient ()
CNP Total inviscid CN ()
CNP1 Time-lagged total inviscid CN ()
CNI Impulsive inviscid CN ()
CNC Circulatory inviscid CN ()
CNf Viscous CN ()
CTf Viscous tangential force coefficient ()
CMf Viscous pitching moment coefficient ()
CMC Circulatory pitching moment coefficient ()
CNV Vortex lift CN ()
CNCRIT Critical CN ()
CPf Stepping parameter moment ()
f Frequency (Hz)
f0 Pitching frequency (Hz)
fn, f1, f2 Separation factor ()
F1, F2 First- and second-order forcing term ()
k Reduced frequency (k=πf0c/V) ()
Kf Stiffness coefficient ()
ks Constant close to Strouhal number ()
M Mach number ()
X, Y, D Deficiency functions ()
c1, c2, c3 Curve-fitting constants ()
a1, a2, a3 Curve-fitting constants ()
S1, S2, S3, α1 Curve-fitting constants ()
A1, A2, b1, b2 Model constants ()
Kα, Tp, Tf, Tv, Tvl Model constants ()
Kv, η, KfC, TmU, TmD Model constants ()
Greek letters
α Angle of attack (rad (unless stated otherwise))
α0 Zero lift α (rad (unless stated otherwise))
αe Effective α (rad (unless stated otherwise))
αf Time-lagged αe (rad (unless stated otherwise))
αCRIT Critical α (rad (unless stated otherwise))
β Mach-number-dependent parameter ()
τv Nondimensional vortex time ()
τ Time constant ()
ζv Vortex lift drag limiting factor ()
Subscripts
n Present sampling time ()
f Viscous lagged value ()
v Vortex lift affected value ()
Superscripts
INV Static inviscid ()
VISC Static viscous ()
I Impulsive ()
CRIT Critical ()
D Dynamic loading ()
D1 First-order correction ()
D2 Second-order correction ()
Data availability

The raw data of the simulation results can be shared by contacting the corresponding author of the paper.

Author contributions

GB developed the new model, designed the studies and conducted the analyses. TL and MA supported the research and provided suggestions and discussion about the manuscript.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

The authors gratefully acknowledge Wobben Research and Development GmbH for providing the research funding through the collaborative joint work DSWind. The measurement data provided from The Ohio State University are highly appreciated.

Financial support

This research has been supported by the Wobben Research and Development GmbH (grant no. DSWind).

This open-access publication was funded
by the University of Stuttgart.

Review statement

This paper was edited by Alessandro Bianchini and reviewed by Khiem V. Truong and Gerard Schepers.

References

Adema, N., Kloosterman, M., and Schepers, G.: Development of a second-order dynamic stall model, Wind Energ. Sci., 5, 577–590, https://doi.org/10.5194/wes-5-577-2020, 2020. a, b, c, d, e, f, g, h, i

Bangga, G.: Three-Dimensional Flow in the Root Region of Wind Turbine Rotors, Kassel University Press GmbH, Kassel, https://doi.org/10.19211/KUP9783737605373, 2018. a

Bangga, G.: Numerical studies on dynamic stall characteristics of a wind turbine airfoil, J. Mech. Sci. Technol., 33, 1257–1262, https://doi.org/10.1007/s12206-019-0225-1, 2019. a

Beddoes, T.: Practical computation of unsteady lift, in: 8th European Rotorcraft Forum, Aix-en-Provence, France, 1982. a, b, c

Carr, L. W., McAlister, K. W., and McCroskey, W. J.: Analysis of the development of dynamic stall based on oscillating airfoil experiments, Tech. rep., NASA TN D-8382, National Aeronautics and Space Administration, Washington, D.C., USA, 1977. a

Dat, R. and Tran, C.: Investigation of the stall flutter of an airfoil with a semi-empirical model of 2 D flow, TP no. 1981-146, ONERA, p. 11, 1981. a, b

Elgammi, M. and Sant, T.: A Modified Beddoes–Leishman Model for Unsteady Aerodynamic Blade Load Computations on Wind Turbine Blades, J. Sol. Energ. Eng., 138, 051009, https://doi.org/10.1115/1.4034241, 2016. a, b, c

Faber, M.: A comparison of dynamic stall models and their effect on instabilities, MS thesis, Delft University of Technology, Delft, available at: https://repository.tudelft.nl/islandora/object/uuid:0001b1eb-c19f-48c3-973d-57eca4996a91 (last access: 5 January 2020), 2018. a, b, c

Galbraith, R.: Return from Airfoil Stall During Ramp-Down Pitching Motions, J. Aircraft, 44, 1856–1864, 2007. a

Gupta, S. and Leishman, J. G.: Dynamic stall modelling of the S809 aerofoil and comparison with experiments, Wind Energy, 9, 521–547, 2006. a, b, c, d

Hansen, M. H., Gaunaa, M., and Madsen, H. A.: A Beddoes–Leishman type dynamic stall model in state-space and indicial formulations, Tech. rep., Risø-R-1354, Risø National Laboratory, Roskilde, Denmark, 2004. a, b, c

Hoffman, M., Ramsay, R., and Gregorek, G.: Effects of grit roughness and pitch oscillations on the NACA 4415 airfoil, Tech. rep., National Renewable Energy Lab., Golden, CO, USA, 1996. a, b, c, d, e

Holierhoek, J., De Vaal, J., Van Zuijlen, A., and Bijl, H.: Comparing different dynamic stall models, Wind Energy, 16, 139–158, 2013. a, b

Janiszewska, J., Ramsay, R., Hoffman, M., and Gregorek, G.: Effects of grit roughness and pitch oscillations on the S814 airfoil, Tech. rep., National Renewable Energy Lab., Golden, CO, USA, 1996. a, b, c, d, e

Khan, M. A.: Dynamic Stall Modeling for Wind Turbines, MS thesis, Delft University of Technology, Delft, available at: https://repository.tudelft.nl/islandora/object/uuid:f1ee9368-ca44-47ca-abe2-b816f64a564f (last access: 5 January 2020), 2018. a

Larsen, J. W., Nielsen, S. R., and Krenk, S.: Dynamic stall model for wind turbine airfoils, J. Fluids Struct., 23, 959–982, 2007. a, b, c, d

Leishman, J.: Validation of approximate indicial aerodynamic functions for two-dimensional subsonic flow, J. Aircraft, 25, 914–922, 1988. a, b, c

Leishman, J. G. and Beddoes, T.: A Semi-Empirical model for dynamic stall, J. Am. Helicopt. Soc., 34, 3–17, 1989. a, b, c, d, e, f, g, h, i, j, k

Martin, J., Empey, R., McCroskey, W., and Caradonna, F.: An experimental analysis of dynamic stall on an oscillating airfoil, J. Am. Helicopt. Soc., 19, 26–32, 1974. a

McAlister, K. W., Carr, L. W., and McCroskey, W. J.: Dynamic stall experiments on the NACA 0012 airfoil, Tech. rep., NASA Technical Paper 1100, National Aeronautics and Space Administration, California, USA, 1978. a

Øye, S.: Dynamic stall simulated as time lag of separation, in: Proceedings of the 4th IEA Symposium on the Aerodynamics of Wind Turbines, Rome, Italy, 1991. a

Pereira, R., Schepers, G., and Pavel, M.: Validation of the Beddoes-Leishman Dynamic Stall Model for Horizontal Axis Wind Turbines using MEXICO data, in: 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 4–7 January 2011, Orlando, USA, AIAA 2011-151, American Institute of Aeronautics and Astronautics (AIAA), Orlando, 2011. a

Petot, D.: Differential equation modeling of dynamic stall, La Recherche Aerospatiale, English Edn., ONERA, 59–72, 1989. a, b

Ramsay, R., Hoffman, M., and Gregorek, G.: Effects of grit roughness and pitch oscillations on the S809 airfoil, Tech. rep., National Renewable Energy Lab., Golden, CO, USA, 1995. a, b, c, d, e

Ramsay, R., Hoffman, M., and Gregorek, G.: Effects of grit roughness and pitch oscillations on the S801 airfoil, Tech. rep., National Renewable Energy Lab., Golden, CO, USA, 1996. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v

Sheng, W., Galbraith, R. M., and Coton, F.: A new stall-onset criterion for low speed dynamic-stall, J. Sol. Energ. Eng., 128, 461–471, 2006. a

Sheng, W., Galbraith, R., and Coton, F.: A modified dynamic stall model for low Mach numbers, J. Sol. Energ. Eng., 130, 031013, https://doi.org/10.1115/1.2931509, 2008. a

Sheng, W., Galbraith, R. A. M., and Coton, F. N.: Applications of low-speed dynamic-stall model to the NREL airfoils, J. Sol. Energ. Eng., 132, 011006, https://doi.org/10.1115/1.4000329, 2010. a

Snel, H.: Heuristic modelling of dynamic stall characteristic, in: Proceedings of the European Wind Energy Conference, Dublin, Ireland, 1997. a, b, c, d, e, f, g, h, i

Tarzanin, F.: Prediction of control loads due to blade stall, J. Am. Helicopt. Soc., 17, 33–46, 1972.  a

Tran, C. and Petot, D.: Semi-empirical model for the dynamic stall of airfoils in view of the application to the calculation of responses of a helicopter blade in forward flight, in: 6th European Rotorcraft Forum, Bristol, UK, 1980. a, b, c

Truong, K.-V.: Modeling aerodynamics for comprehensive analysis of helicopter rotors, in: Proceedings of 42nd European Rotorcraft Forum, Lille, France, 2016. a

Truong, V.: A 2-d dynamic stall model based on a hopf bifurcation, in: Proceedings of 19th European Rotorcraft Forum, Cernobbio, Italy, 1993. a, b

Wang, Q. and Zhao, Q.: Modification of Leishman–Beddoes model incorporating with a new trailing-edge vortex model, Proc. Inst. Mech. Eng. Pt G, 229, 1606–1615, 2015. a

Download
Short summary
Robust and accurate dynamic stall modeling remains one of the most difficult tasks in wind turbine load calculations despite its long research effort in the past. The present paper describes a new second-order dynamic stall model for wind turbine airfoils. The new model is robust and improves the prediction for the aerodynamic forces and their higher-harmonic effects due to vortex shedding but also provides improved predictions for pitching moment and drag.