Development and Validation of the IAG Dynamic Stall Model in State-Space Representation for Wind Turbine Airfoils

1. Introduction

Dynamic stall can occur when there is a strong variation of the inflow conditions due to environmental effects and turbine operation strategy, e.g., yaw misalignment, wind turbulence, shear & gusts, tower shadow and strong aeroelastic effects of the blade. Studies clearly highlighted [1,2,3,4,5,6,7,8,9,10] that the aerodynamic loads can be significantly different than the stationary conditions. The dynamic stall phenomenon is usually initiated by an increase of lift with increasing angle of attack ($\alpha$) past the corresponding static stall angle. This is associated with the generation of a leading edge vortex (LEV). This vortical structure creates a suction effect which enhances the circulation effect on the airfoil [9]. LEV is further convected downstream toward the trailing edge. This causes a further increase of the lift and drag forces, at the same time the pitching moment becomes more negative before stall occurs. The dynamic loads associated with the dynamic stall effects can significantly alter wind turbine loads and have a strong influence on the aerodynamic damping. As a direct consequence, considering the dynamic stall effects in engineering calculations is essential for correcting the aerodynamic loads acting on wind turbines, both during power production and stand-still cases. A correct representation of the aerodynamic force hysteresis is especially critical to accurately model the aerodynamic damping that has a significant impact on the turbine aeroelastic stability. This is especially important for modern large rotors because the blade is longer, more flexible and slender compared to older turbine designs.

Engineering modeling relies on dynamic stall models to include this complex phenomenon. These models can produce reasonable results while still maintaining low computational burden. The models from Beddoes and Leishman (BL) and its derivatives [11,12,13,14,15,16,17,18,19,20] are considered to be the industry standard for wind turbine design process. The models are commonly derived by combining a time delay of the angle of attack with an approximate solution of flow separation during the dynamic conditions. The original BL model was dedicated for high speed aerodynamics applications and included the Mach number effects in most flow states. This model was further simplified by Hansen et al. [12] by removing the compressiblity effects and transforming the model into a state-space representation. However, this model ignores the effect of the leading edge vortex by arguing that the airfoil thickness for wind turbine is no less than 15%. Examples of the other available models are provided by Øye [21], Tran & Petot (ONERA model) [22], Tarzanin (Boeing-Vertol model) [23] and Snel model [24]. Furthermore, the interest to model a higher harmonic effect of the dynamic stall polar increases nowadays. Snel model [24] is one of the few models which can be done to include the effect. This model was further extended in [25]. Bangga et al. [13] combined the state-of-the-art BL model and the second order term of the Snel model to predict the unsteady characteristics of wind turbine airfoils, hereby referred to as the “IAG dynamic stall model”. Several improvements were made to enhance the accuracy of the first order and the second order terms [13]. The first order term of the IAG model was built based on the indicial formulations while the second order term was derived using a state-space representation.

The wind turbine design tool Bladed [26] traditionally employs the incompressible version of the BL model with some modifications [12,26]. Although the adopted model works at moderate angles of attack, the performance of this type of model starts to deteriorate when it is used to predict deep stall conditions, which may lead to an incorrect prediction of the aerodynamic damping. The first order term of the IAG model was recently implemented in a development version of Bladed and this paper is intended to report the implementation and verification processes. The IAG model is transformed from the indicial type of formulation into a state-space representation, allowing the model to be linearized to enable steady-state stability analysis. In this paper, a consistent step-by-step procedure will be presented and the physical justification for the modifications applied will be highlighted. The new model is compared with other dynamic stall models in Bladed, including the widely used incompressible BL model [12,26] and Øye [21], and validated against measurement data of wind turbine airfoils with different relative thicknesses and flow conditions, allowing a full assessment of its suitability for different purposes.

The paper is organized as following. Section 2 describes the mathematical formulation of the IAG model in the original indicial and the transformed state-space representations. An automated procedure to determine the polar gradient is also documented. Then, in Section 3 assessments are carried out on the performance of the IAG dynamic stall model in comparison with measurement data and other models. The new model is further tested at various flow conditions and airfoils without further calibrating the constants to examine its robustness for different situations. Finally, all results will be concluded in Section 4.

2. Mathematical Foundation

The IAG model [13] is comprised of a first order and second order model. The first order term is based the Beddoes-Leishman model [27,28] with some improvements. The second order term adopts the improved second order Snel model to accommodate the higher harmonic effects to be modeled. In this paper, the first order term will be the focus of the investigation because this determines the overall shape of the hysteresis curve. The modeling strategy of the first order term is divided into several aerodynamic states: (1) attached flow state, (2) separated flow state and (3) vortex lift state. The mathematical foundations to both implementation types are further described in Section 2.1 and Section 2.2.

2.1. IAG model in indicial formulation

The loads of the attached flow state are assumed to originate from two main sources; (a) the circulatory loading which quickly builds up to the steady state value and (b) impulsive loading which is derived from the piston theory. For a lifting surface exposed to unsteady change of the angle of attack ${\alpha }_n$, there is a time delay between the actual effective angle of attack (${\alpha }_{e_n}$) seen by the structure. This aspect can be represented as:

$$\begin{array}{c}\hfill {\alpha }_{e_n}={\alpha }_n-X_n-Y_n\end{array}$$

(1)

with n being the current sample time (current value of the actual angle o attack). The deficiency functions are described by:

$$\begin{array}{c}{X}_n=X_{n-1}\exp \left(-b_1{\beta }^2\Delta s\right)+A_1\Delta {\alpha }_n\exp \left(-b_1{\beta }^2\Delta s/2\right)\end{array}$$

(2)

$$\begin{array}{c}{Y}_n=Y_{n-1}\exp \left(-b_2{\beta }^2\Delta s\right)+A_2\Delta {\alpha }_n\exp \left(-b_2{\beta }^2\Delta s/2\right)\end{array}$$

(3)

Note that variable time is represented by a dimensionless parameter $s=2Vt/c$. This defines the relative distance traveled by the airfoil in terms of semi chord. The variables V, t and c describe the free stream wind speed, time and chord length, respectively. In this sense, $\Delta {\alpha }_n$ and $\Delta s$ are given by:

$$\begin{array}{c}\Delta {\alpha }_n={\alpha }_{n+1}-{\alpha }_n\end{array}$$

(4)

$$\begin{array}{c}\Delta s=s_n-s_{n-1}\end{array}$$

(5)

$$\begin{array}{c}\Delta t=t_n-t_{n-1}.\end{array}$$

(6)

Furthermore, $b_1$, $b_2$, $A_1$ and $A_2$ are defined as the attached flow constants. The compressibility effects are included in the formulation by adopting $\beta =\sqrt{1-M^2}$, with M being the Mach number.

The attached flow normal force coefficient is obtained by:

$$C_{N_n}^P=C_{N_n}^C+C_{N_n}^I,$$

(7)

which combines the circulatory and attached flow components. The circulatory normal force can be obtained using

$$C_{N_n}^C={\displaystyle \frac{d{C}_N}{d\alpha }}({\alpha }_{e_n}-{\alpha }_0^{INV})$$

(8)

with ${\alpha }_{0_{INV}}$ being the angle of attack for zero inviscid normal force and $d{C}_N/d\alpha$ being the normal force gradient of the polar data within the linear regime. This is in contrast with the original implementation of the Beddoes-Leishman model [27,28] which disregarded the use of ${\alpha }_0^{INV}$. However, this term is important when the airfoil has a finite camber. This has been pointed out as well by Hansen et al. [12]. The impulsive normal force is calculated based on

$$C_{N_n}^I={\displaystyle \frac{4K_{\alpha }{T}_I}{M}}\left({\displaystyle \frac{\Delta {\alpha }_n}{\Delta t}}-D_n\right).$$

(9)

with $T_I=Mc/V$. $D_n$ represents the deficiency function and is given by

$$D_n=D_{n-1}\exp \left({\displaystyle \frac{-\Delta t}{K_{\alpha }{T}_I}}\right)+\left({\displaystyle \frac{\Delta {\alpha }_n-\Delta {\alpha }_{n-1}}{\Delta t}}\right)\exp \left({\displaystyle \frac{-\Delta t}{2K_{\alpha }{T}_I}}\right),$$

(10)

with $K_{\alpha }$ being a constant for determining the strength of the impulsive normal force effect.

To obtain the value of the unsteady viscous normal force, the Kirchhoff equation is used to reconstruct the aerodynamic properties as:

$$C_{N_n}^f={\displaystyle \frac{d{C}_N}{d\alpha }}{\left({\displaystyle \frac{1+\sqrt{f_{2_n}}}{2}}\right)}^2({\alpha }_{e_n}-{\alpha }_0^{VISC})+C_{N_n}^I.$$

(11)

Note that the impulsive normal force effect is included in the formulation. Also notice that the zero normal force angle of attack for viscous polar data (${\alpha }_{0_{VISC}}$) is used here. In the original formulation of the IAG model [13], the attached flow data is obtained from the inviscid panel-method XFOIL calculations, thus the zero normal force angle of attack can be different with the value obtained from the viscous polar data.

The unsteady trailing edge separation point ($f_{2_n}$) in Equation (11) is obtained from:

$$f_{2_n}=f_n-D_{f_n},$$

(12)

which depends on two parameters $f_n$ and $D_{f_n}$. The first parameter is defined by an inverse Kirchhoff equation as:

$$f_n={\left(2\sqrt{{\displaystyle \frac{C_{N_n}^{VISC}}{{\displaystyle \frac{d{C}_N}{d\alpha }}({\alpha }_{f_n}-{\alpha }_0^{VISC})}}}-1.0\right)}^2,$$

(13)

and can be computed based on the pressure lag response of the angle of attack:

$${\alpha }_{f_n}={\alpha }_0^{INV}+\left({\displaystyle \frac{C_{N_n}^{P1}}{d{C}_N/d\alpha }}\right).$$

(14)

with

$$C_{N_n}^{P1}=C_{N_n}^P-D_{p_n},$$

(15)

and $D_{p_n}$ being represented as:

$$D_{p_n}=D_{p_{n-1}}\exp \left(-{\displaystyle \frac{\Delta s}{T_p}}\right)+\left(C_{N_n}^P-C_{N_{n-1}}^P\right)\exp \left(-{\displaystyle \frac{\Delta s}{2T_p}}\right).$$

(16)

The last term of Equation (12) is given by the following time dependent equation (which also depends on $f_n$):

$$D_{f_n}=D_{f_{n-1}}\exp \left(-{\displaystyle \frac{\Delta s}{T_f}}\right)+\left(f_n-f_{n-1}\right)\exp \left(-{\displaystyle \frac{\Delta s}{2T_f}}\right).$$

(17)

The parameters $T_p$ and $T_f$ are the constants for determining the pressure lag and flow separation effects, respectively.

The last part of the first order term represents the vortex lift effect which can be computed by:

$$C_{N_n}^V=\left\{\begin{array}{c}{C}_{N_{n-1}}^V\exp \left(-{\displaystyle \frac{\Delta s}{T_v}}\right)+\left(C_{V_n}-C_{V_{n-1}}\right)\exp \left(-{\displaystyle \frac{\Delta s}{2T_v}}\right);\hfill \\ \hfill & \mathrm{if}0<{\tau }_{v_n}<T_{vl}\hfill \\ C_{N_{n-1}}^V\exp \left(-{\displaystyle \frac{\Delta s}{T_v}}\right);\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(18)

Notice that Equation (18) depends on the value of $C_{V_n}$ and this can be estimated from the circulatory normal force effect as:

$$C_{V_n}=C_{N_n}^C\left(1-{\displaystyle \frac{1}{4}}{\left(1+\sqrt{f_{2_n}}\right)}^2\right)$$

(19)

Variables $T_v$ and $T_{vl}$ represent the vortex decay and vortex travel time constants, respectively. The non-dimensional vortex time (${\tau }_{v_n}$) itself is important for determining the vortex lift influence and can be calculated using

$${\tau }_{v_n}=\left\{\begin{array}{cc}{\tau }_{v_{n-1}}+0.45{\displaystyle \frac{\Delta t}{c}}V;\hfill & \mathrm{if}{C}_{N_n}^{P1}>C_N^{CRIT}\hfill \\ 0;\hfill & \mathrm{if}{C}_{N_n}^{P1}<C_N^{CRIT}\mathrm{and}\Delta {\alpha }_n\ge 0\hfill \\ {\tau }_{v_{n-1}};\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(20)

$C_N^{CRIT}$ is an important parameter which defines the inviscid critical static normal force. This is usually indicated by the break of the pitching moment polar at the critical angle of attack ${\alpha }_n^{CRIT}$. The magnitude of $C_N^{CRIT}$ can be computed as:

$$C_N^{CRIT}={\displaystyle \frac{d{C}_N}{d\alpha }}({\alpha }_n^{CRIT}-{\alpha }_0^{INV}).$$

(21)

The total normal force contribution can be obtained by computing

$$C_{N_n}^D=C_{N_n}^f+C_{N_n}^V.$$

(22)

In the IAG model, the tangential (chordwise) force is simply obtained from the static polar data at the time-lagged angle of attack ${\alpha }_{f_n}$ by:

$$C_{T_n}^D=C_T^{VISC}\left({\alpha }_{f_n}\right).$$

(23)

Finally, the lift and drag forces can be computed as:

$$\begin{array}{c}{C}_{L_n}^D=C_{N_n}^D\cos {\alpha }_n-C_{T_n}^D\sin {\alpha }_n\end{array}$$

(24)

$$\begin{array}{c}{C}_{D_n}^D=C_{N_n}^D\sin {\alpha }_n+C_{T_n}^D\cos {\alpha }_n\end{array}$$

(25)

I n the IAG model, a limiter is applied to the drag force to correct the hysteresis effects. This is done through a parameter ${\zeta }_n$ which is defined as:

$${\zeta }_n=\frac{1}{\pi }{\displaystyle \frac{d{C}_N}{d\alpha }}{\left({\displaystyle \frac{1+\sqrt{f_n}}{2}}\right)}^2.$$

(26)

Drag hysteresis is observed to occur when ${\zeta }_n\lesssim {\zeta }_v$, with ${\zeta }_v$ being a constant with a value of 0.76. The adopted correction reads:

$$C_{D_n}^D=\left\{\begin{array}{cc}1.2C_{D_n}^{VISC};\hfill & \mathrm{if}{C}_{D_n}^D>1.2C_{D_n}^{VISC}\mathrm{and}\left(C_{N_n}^P-C_{N_{n-1}}^P\right)\ge 0.0\hfill \\ \\ C_{D_n}^{VISC};\hfill & \mathrm{if}\left(C_{N_n}^P-C_{N_{n-1}}^P\right)<0.0\hfill \\ \\ C_{D_n}^D;\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(27)

As for the pitching moment coefficient, the total dynamic response obtained from all contributions is obtained by computing

$$C_{M_n}^D=C_{M_n}^f+C_{M_n}^V+C_{M_n}^C$$

(28)

where the separated, vortex and circulatory contributions of the pitching moments are defined by $C_{M_n}^f$, $C_{M_n}^V$ and $C_{M_n}^C$, respectively.

In the original Beddoes-Leishman formulation, the pitching moment was obtained by a curve fitting procedure. This is not straightforward as the user needs to perform curve fitting of the polar data. In the IAG model, the separated moment coefficient is easily obtained from the static viscous polar data at ${\alpha }_{f_n}$, that reads

$$C_{M_n}^f=C_M^{VISC}\left({\alpha }_{f_n}\right).$$

(29)

In this sense, the moment coefficient can be reconstructed easily without the need to adjust the parameters in advance, minimizing the user error. The vortex contribution of the dynamic moment coefficient can be formulated as

$$C_{M_n}^V=-C_{P{v}_n}{C}_{N_n}^V.$$

(30)

This depends on the idealized variation of the center of pressure and can be estimated as

$$C_{P{v}_n}=K_v\left(1-\cos \left({\displaystyle \frac{\pi {\tau }_v}{T_{vl}}}\right)\right).$$

(31)

The last term represents the circulatory moment. To model this, a relatively simple approach is introduced by applying a time delay on the circulatory moment response as:

$$C_{M_n}^C=\left\{\begin{array}{l}{C}_{M_{n-1}}^C\exp \left(-{\displaystyle \frac{\Delta s}{T_M^U}}\right)-C_{P{f}_n}\left(C_{V_n}-C_{V_{n-1}}\right)\exp \left(-{\displaystyle \frac{\Delta s}{2T_M^U}}\right);\\ \mathrm{if}{\tau }_{v_n}<T_{vl}\mathrm{and}\Delta {\alpha }_n\ge 0\hfill \\ \\ C_{M_{n-1}}^C\exp \left(-{\displaystyle \frac{\Delta s}{T_M^D}}\right)-C_{P{f}_n}\left(C_{V_n}-C_{V_{n-1}}\right)\exp \left(-{\displaystyle \frac{\Delta s}{2T_M^D}}\right);\\ \mathrm{if}\Delta {\alpha }_n<0\\ \\ C_{M_{n-1}}^C;\mathrm{otherwise}\end{array}\right.$$

(32)

2.2. IAG model in state-space representation

In the present paper, the IAG model is transformed into a state space formulation. The main modeling strategy is done similar to the incompressible Beddoes-Leishman model [12]. However, the formulations are defined in forms of the normal force coefficient, in contrast to Ref. [12] which used the lift coefficient directly. The same principles are adopted as in the indicial formulations, i.e., dividing the solutions into attached flow, separated flow and vortex lift states.

The attached flow solutions are represented into two ordinary differential equations (ODEs) which are comparable to the deficiency functions already presented in Section 2.1. These are described by:

$$\begin{array}{c}\left(\frac{c}{2V}\right){\dot{x}}_1\left(t\right)=b_1{A}_1{\alpha }_n^{3/4}-\left[b1+\left(\frac{c\dot{V}}{V^2}\right)\right]x_1\left(t\right)\end{array}$$

(33)

$$\begin{array}{c}\left(\frac{c}{2V}\right){\dot{x}}_2\left(t\right)=b_2{A}_2{\alpha }_n^{3/4}-\left[b2+\left(\frac{c\dot{V}}{V^2}\right)\right]x_2\left(t\right)\end{array}$$

(34)

Here, $x_{1,2}\left(t\right)$ and ${\dot{x}}_{1,2}\left(t\right)$ represent the state and the state derivative of the ODEs, respectively. Note that the definition of the angle of attack being adopted in the ODEs uses the angle of attack measured at the collocation point (3/4 of the chord measured from the leading edge). This has no influence for two dimensional pitching airfoil but can be different with the quarter chord position for wind turbine blade sections affected by structural flexibility. As suggested in [12], the added mass term $\left(\frac{c\dot{V}}{V^2}\right)$ is added to include the effects for varying wind speed in wind turbine cases. The effective angle of attack is calculated using:

$$\begin{array}{c}\hfill {\alpha }_{e_n}={\alpha }_n^{3/4}(1-A_1-A_2)+x_1\left(t\right)+x_2\left(t\right)\end{array}$$

(35)

The circulatory normal force can be obtained using

$$C_{N_n}^C={\displaystyle \frac{d{C}_N}{d\alpha }}sin({\alpha }_{e_n}-{\alpha }_0^{VISC}).$$

(36)

Note that Equation (36) looks different compared to Equation (8). One of it is the usage of ${\alpha }_0^{VISC}$ instead of ${\alpha }_0^{INV}$. Although the magnitude of these two parameters can be different depending on the airfoils, the usage of the inviscid value is not convenient because users need to provide two polar data inputs, one the static viscous data and the other is the inviscid data (e.g., from XFOIL). The second key difference is the usage of a sinusoidal form to determine the force. This is beneficial for high angle of attack calculations since data demonstrates that a linear model does not capture the nonlinear inviscid effects. An illustration of the effects for lift is given in Figure 1.

To add the contribution for the impulsive effect, Equation (9) can be used. Two options can be adopted in this sense, either by transforming the deficiency function in Equation (10) as a state-space representation or completely disregarding the Mach number. The latter is chosen because this implies that $D_n$ will be zero. As a consequence, it reduces the number of states to be solved by the integrator, i.e., faster computational speed. This is formulated as:

$$C_{N_n}^I=\frac{4K_{\alpha }c}{V}{\dot{\alpha }}_n.$$

(37)

Finally, the attached flow normal force coefficient is obtained by:

$$C_{N_n}^P=C_{N_n}^C+C_{N_n}^I.$$

(38)

The next set of ODE is used to determine the time-lagged pressure response of the normal force, as similarly done in Equation (15). This is formulated as:

$$\begin{array}{c}\left(\frac{c}{2V}\right)T_p{\dot{x}}_3\left(t\right)=-x_3\left(t\right)+C_{N_n}^P.\end{array}$$

(39)

Variable $x_3\left(t\right)$ is the solution for the pressure-lagged normal force or $C_{N_n}^{P1}$ in Equation (15). This parameter is used to calculate the delayed angle of attack

$${\alpha }_{f_n}={\alpha }_0^{VISC}+\left({\displaystyle \frac{x_3\left(t\right)}{d{C}_N/d\alpha }}\right).$$

(40)

By using the inverse of the Kirchoff equation, the position of separation point at ${\alpha }_{f_n}$ can be determined by:

$$f_n={\left(2\sqrt{{\displaystyle \frac{C_{N_n}^{VISC}}{{\displaystyle \frac{d{C}_N}{d\alpha }}sin({\alpha }_{f_n}-{\alpha }_0^{VISC})}}}-1.0\right)}^2.$$

(41)

Again, a sinusoidal approximation is also adopted in Equation (41). This leads to the fourth state-space equation to further delay the obtained separation position as:

$$\begin{array}{c}\left(\frac{c}{2V}\right)T_f{\dot{x}}_4\left(t\right)=-x_4\left(t\right)+f_n,\end{array}$$

(42)

which is comparable to Equation (12).

The unsteady viscous normal force can be reconstructed using the Kirchhoff equation as already done in Equation (11). However, here a sinusoidal approximation is consistently adopted instead of using a linear model and is denoted as:

$$C_{N_n}^f={\displaystyle \frac{d{C}_N}{d\alpha }}{\left({\displaystyle \frac{1+\sqrt{x_4\left(t\right)}}{2}}\right)}^{2}sin({\alpha }_{e_n}-{\alpha }_0^{VISC})+C_{N_n}^I.$$

(43)

Note that $x_4\left(t\right)$ in Equation (11) is comparable to $f_{2_n}$ in Equation (43).

The last term to consider is the vortex lift effect due to the presence of the leading edge vortex. Similar to the description in Section 2.1, this is calculated by the following ODE:

$$\begin{array}{c}\left(\frac{c}{2V}\right)T_v{\dot{x}}_5\left(t\right)=-x_5\left(t\right)+\left(\frac{c}{2V}\right)T_v{\dot{C}}_{V_n},\end{array}$$

(44)

with $x_5\left(t\right)$ specifying the normal force due to the vortex lift effect, which equals to $C_{N_n}^V$ in Equation (18). To solve this equation, information about the rate of change of $C_{V_n}$ is needed. To do so, by applying a chain rule, a time derivative of Equation (19) is calculated as:

$${\dot{C}}_{V_n}={\displaystyle \frac{d{C}_N}{d\alpha }}{\dot{\alpha }}_{e_n}\left(1-{\displaystyle \frac{1}{4}}{\left(1+\sqrt{x_4\left(t\right)}\right)}^2\right)-{\displaystyle \frac{1}{4}}{C}_{N_n}^C\left(1+\sqrt{x_4\left(t\right)}\right)\frac{{\dot{x}}_4\left(t\right)}{\sqrt{x_4\left(t\right)}}.$$

(45)

Equation (45) is only applied when there is a dynamic vortex state being active, otherwise ${\dot{C}}_{V_n}=0$. Vortex state itself is defined when $x_3\left(t\right)$ is greater/smaller than the critical normal force:

$$C_N^{CRIT}=\left\{\begin{array}{cc}{C}_N^{MAX};\hfill & \mathrm{if}\mathrm{Positive}\mathrm{Stall}\hfill \\ C_N^{MIN};\hfill & \mathrm{if}\mathrm{Negative}\mathrm{Stall}\hfill \end{array}\right.$$

(46)

which refers to the maximum normal force (or minimum depending on the case) obtained from the static data. This is different with Equation (21) and the effect will be discussed in Section 3.2. Another important consideration to activate ${\dot{C}}_{V_n}$ is that it has to be during the upstroke motion. This implies that:

$$\mathrm{Vortex}\mathrm{state}=\left\{\begin{array}{cc}\mathrm{True};\hfill & \mathrm{if}\mathrm{Upstroke}\&x_3\left(t\right)>C_N^{CRIT}\&0<{\tau }_{v_n}<T_{vl}\hfill \\ \mathrm{False};\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(47)

The upstroke state can be determined by $\Delta {\alpha }_n$, but this will be different when the case being considered is positive stall ($\Delta {\alpha }_n>0$) or negative stall ($\Delta {\alpha }_n<0$). A positive stall case can be identified when $x_3\left(t\right)\ge 0$, while a negative stall case may be detected when $x_3\left(t\right)<0$. For the backwinded case (when $\parallel \alpha \parallel$ is greater than 90°), the definition of upstroke is flipped. Now, the non-dimensional vortex time may be calculated using

$${\tau }_{v_n}=\left\{\begin{array}{cc}{\tau }_{v_{n-1}}+0.45{\displaystyle \frac{\Delta t}{c}}V;\hfill & \mathrm{if}\mathrm{Vortex}\mathrm{state}\hfill \\ {\tau }_{v_{n-1}}exp\left(-{\displaystyle \frac{\Delta t}{c}}V\right);\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(48)

The second term is only adopted to avoid a sharp jump of the value numerically, and may be set to zero in the implementation if desired. For very large angle of attack approaching 90°, it is assumed that the vortex lift effect vanishes. Therefore, the value is slowly relaxed toward zero from 45° to 75°. This consideration is based on the fact that lift generally increases again after stall up to a full flow separation at around 60°, see e.g. [29].

Similar to the indicial formulation in Section 2.1, the total normal force contribution can be obtained by computing

$$C_{N_n}^D=C_{N_n}^f+C_{N_n}^V.$$

(49)

On top of that, the tangential (chordwise) force is obtained directly from the static polar data at the time-lagged angle of attack ${\alpha }_{f_n}$ by:

$$C_{T_n}^D=C_T^{VISC}\left({\alpha }_{f_n}\right).$$

(50)

In the present model, only lift coefficient is calculated from $C_{N_n}^D$ and $C_{T_n}^D$ transformation as:

$$\begin{array}{c}{C}_{L_n}^D=C_{N_n}^D\cos {\alpha }_n-C_{T_n}^D\sin {\alpha }_n.\end{array}$$

(51)

The calculation for drag coefficient is done differently by applying the following relation:

$$\begin{array}{c}{C}_{D_n}^D=C_{D_n}^{VISC}+\left({\alpha }_n-{\alpha }_{e_n}\right)C_{N_n}^C\hfill \\ \hfill +\left(C_{D_n}^{VISC}-C_{D_0}^{VISC}\right)\left[{\left(\frac{1-\sqrt{x_4\left(t\right)}}{2}\right)}^2-{\left(\frac{1-\sqrt{f^{VISC}}}{2}\right)}^2\right]+x_5\left(t\right)sin{\alpha }_n\end{array}$$

(52)

In Equation (52), variables $C_{D_n}^{VISC}$, $C_{D_0}^{VISC}$ and $f^{VISC}$ represent the static value of drag coefficient, drag level at zero normal force angle of attack and static separation position, respectively. Furthermore, it is assumed that drag force shows strong hysteresis only near stall regime. Therefore, the magnitude of drag is returned to the static value once the angle of attack is greater than 30°. This is done through a linear blending up to 45°. Note that a sudden drag increase occurs during stall and this takes place at $\parallel \alpha \parallel <{30}^{\circ }$ for most airfoils.

Lastly, the pitching moment coefficient is calculated in a similar manner as the indicial representation in Equation (53) that reads:

$$C_{M_n}^D=C_{M_n}^f+C_{M_n}^V+C_{M_n}^C.$$

(53)

However, the last term is not computed based on two deficiency function as in Equation (32) to avoid unnecessary computational effort. This is modeled as an added mass instead, which can be defined as:

$$C_{M_n}^C=-\frac{\pi c{\dot{\alpha }}_n}{4V}.$$

(54)

For backwinded airfoil orientation (when $\parallel \alpha \parallel$ is greater than 90°), the sign of Equation (54) becomes positive. The calculations for $C_{M_n}^f$ and $C_{M_n}^V$ are exactly the same as done in Equation (29) and Equation (30), respectively. Similar as for the drag coefficient, the magnitude of the pitching moment is returned to the static value once the angle of attack is greater than 30° through a linear blending with the static data.

2.3. Adopted constants

The following constants are applied in the tested IAG dynamic stall model. These values are kept constant throughout the paper and are listed in Table 1. Note that the state-space IAG model has fewer number of constants compared to the indicial one (no $K_f^C$, $T_m^U$, $T_m^D$ and ${\zeta }_s$). For the indicial IAG dynamic stall model, the critical angle of attack plays a major role and the values are given explicitly in Table 2. This corresponds to the angle when the viscous pitching moment breaks. In contrast, this information is not required in the newly implemented state-space IAG model because the values are determined automatically based on the maximum viscous normal force coefficient.

2.4. Automatic determination of the normal force gradient

Traditionally, Bladed determines the normal force gradient by using a two-point approximation. This approach takes the gradient computed between two distinct points of the polar data: (1) when $C_N=0$ and (2) when $\alpha =0$. This approach works for most airfoils but the validity breaks when there is a gradient change between these two distinct points as illustrated in Figure 2. To cover this aspect, a new method called the "linear fit approach" is adopted by incorporating more data points for determining the normal force gradient. First of all, local gradients for all data point between ${\alpha }_0^{VISC}$ and $\alpha =7^{\circ }$ are computed. Within this area, it is assumed that the flow is mostly attached. From the obtained gradients, all data points are binned and only the points with gradients between $1.8\pi$ and $2.5\pi$ are retained. A linear fit is then searched using the obtained points to find the most appropriate gradient of the polar data.

3. Results and Discussion

3.1. Test cases and treatment of the input data

The present studies were carried out based on the experimental data obtained from the measurement campaign conducted at the Ohio State University (OSU) [4,5,6,7]. The calculations were conducted mainly on Selig airfoil family developed specially for wind turbine applications. Most discussions in this paper are focused on the S809 airfoil having a 21% of relative thickness ($\delta /c$), but are also extended for the S801 ($\delta /c=13.5\%$) and S814 ($\delta /c=24\%$) airfoils to evaluate the influence of the airfoil thickness. The evaluations of the results are performed on the airfoil with the Reynolds numbers ($Re$) of 750K, 1000K and 1250K according to the available measurement data. Note that unless stated otherwise, the figures presented in subsequent sections correspond to the Reynolds number of 750K case. Most test cases are for the airfoils employed with a leading edge grit (turbulator) to enable the "soiled" effects on a wind turbine blade unless stated otherwise.

Generally, it is a common practice to replicate the experimental campaign of dynamic stall by applying a sinusoidal pitching motion in the simulation environment. However, this contains some uncertainty because it is not always possible to have a perfect sinusoidal motion in the experiment. This has been demonstrated by Bangga et al. [13] and they devised that this can lead to an additional discrepancy contribution. On the other hand, the oscillator orientation angle is not always recorded in fine resolution, making it not possible to supply the data directly into the simulation environment. A cubic fitting on the inflow data was found best representing the inflow time series [13], and this is also highlighted in Figure 3. Therefore, the same procedure is adopted in the present paper to obtain the best consistent setup with the experimental campaign.

3.2. Consistency with the original indicial formulation

To evaluate the consistency of the newly implemented IAG dynamic stall model in Bladed, a comparison against the original indicial formulation in the B-GO code [13,30,31] is carried out. The results of the dynamic stall calculations in Bladed using the state-space version of the IAG model are presented in Figure 4. It can be seen that the implemented state-space model using the original constant shows a reasonable consistency with the original model. Some deviations are observed especially near the maximum angle of attack. This is mainly caused by the way the critical angle of attack is defined in both models. In the newly implemented model, this is automatically determined based on the maximum normal force location, while for the original model it was observed manually at a location where the pitching moment breaks (pitching moment stall). Furthermore, the critical normal force in the indicial formulation is defined as the "inviscid" normal force based on the critical angle of attack itself. This is different with the present implementation which uses the maximum "viscous" normal force directly obtained from the polar data. When the critical normal force in B-GO is set to use the maximum viscous normal force, the agreement between both codes is significantly improved as shown in Figure 5.

Another minor source of discrepancy is the linear assumption of the attached flow polar reconstruction of the original model against sinusoidal assumption used in the present implementation. Drag deviation is a direct consequence for the normal force modeling. As for the pitching moment, no circulatory pitching moment model as in Equation (32) is used in the state-space formulation. This term is replaced with the added mass effect formulated using Equation (54), which yields a similar impact without having to provide two additional constants as in the original IAG model.

Nevertheless, the newly developed model is generally in a reasonable agreement with the original model as well as the experimental data. Note that the higher harmonic effects in the experimental data are not captured in both models because the second order model is not implemented. Only the first order term is the focus in the present studies since this term governs the main shape of the hysteresis curve.

3.3. The effects of mean angle of attack

In this section, Bladed calculations employing the IAG dynamic stall model are compared against the experimental data at three different values of the mean angle of attack, as depicted in Figure 6. Bladed calculations employing the incompressible Beddoes-Leishman model and the Øye model are also presented for comparison. The constants of the other models remain default and are documented in [26].

It can be seen that the IAG model outperforms the other available dynamic stall model in Bladed for all three different mean angles of attack, but this is especially pronounced in the deep stall conditions. The IAG model is able to predict the hysteresis effects well, not only for the lift coefficient but also for the drag coefficient and the pitching moment coefficient. Two peaks of the lift force at $\overline{\alpha }={20}^{\circ }$ are not captured by the model because this is inherently the second order effects due to vortex shedding. This is not modeled in the first order term and requires the inclusion of the second order term.

3.4. The effects of reduced frequency

The verification of the IAG dynamic stall model at various reduced frequencies is given in this section. This parameter holds a very important position in dynamic stall effects because it determines the degree of unsteadiness for the case. The results are presented in Figure 7 for three different reduced frequencies: $k=0.033$, $k=0.064$ and $k=0.098$. It can be clearly seen that the IAG dynamic stall model scales well when the reduced frequency is varied, where the stronger reduced frequency creates a stronger hysteresis loop. The same behavior is observed for the experimental data which align well with the IAG dynamic stall results. In contrast, the existing dynamic stall models in Bladed (BL and Øye) fail to predict such characteristics accurately.

3.5. The effects of excitation amplitude

Figure 8 presents the dynamic polar comparison between Bladed calculations and experimental data for two different excitation amplitudes ($\Delta \alpha =5.5^{\circ }$ and $\Delta \alpha =5.5^{\circ }$). The highest available mean angle of attack ($\overline{\alpha }={20}^{\circ }$) is selected for the investigation because the strongest dynamic stall effect is seen at the post stall region. The already implemented BL and Øye models in Bladed are observed to underestimate the dynamic hysteresis effects significantly. For $\Delta \alpha =5.5^{\circ }$, it is even observed that these models do not predict any force unsteadiness. This characteristic does not agree well with the experimental data, which shows force hysteresis for both investigated amplitudes. The IAG model, again, shows a good agreement with the experimental data for all force components. As already observed in preceding sections, the higher harmonic effects are not captured without including the second order term for vortex shedding modeling.

3.6. Sensitivity to airfoil thickness

To demonstrate the generality of the developed IAG dynamic stall model, the studies are extended to cover airfoils with different relative thicknesses. The results for all three force components are presented in Figure 9. It can be clearly seen that the airfoil thickness has a strong influence on the dynamic stall characteristics. The stall for the thinnest airfoil (S801) is the strongest among the considered airfoils in Figure 9. This is in agreement with the studies carried out in [10]. The IAG model is able to predict the stall angle location in response to the airfoil thickness to some degree. The generic shape of the polar hysteresis is predicted by the IAG model for the investigated airfoils, which clearly outperforms the BL and Øye models.

3.7. Model performance for negative stall

In the real operation, wind turbine airfoil sections are exposed not only to positive angle of attack regimes but also negative regimes. This is especially true when the turbine is operating in stand-still conditions at a large yaw misalignment on top of strong inflow turbulence. In this section, the performance of the IAG model in negative stall prediction will be examined. To be able to do so, negative polar data needs to be generated because the OSU experimental campaign only covers the negative stall regimes. The available static and dynamic polar data sets are transformed into the negative data by applying Equation (55) to Equation (58). The synthesized static negative polar for the S809 airfoil is displayed in Figure 10.

$$\begin{array}{c}{\alpha }^{Neg}=-\alpha +2{\alpha }^{CL=0}\end{array}$$

(55)

$$\begin{array}{c}{C}_L^{Neg}=-C_L\end{array}$$

(56)

$$\begin{array}{c}{C}_D^{Neg}=C_D\end{array}$$

(57)

$$\begin{array}{c}{C}_M^{Neg}=-C_M+2C_M^{CL=0}\end{array}$$

(58)

Using the reconstructed negative data, Bladed calculations were carried out by adopting the IAG, BL and Øye models. Two mean angles of attack were considered in the simulations, having $k=0.064$, $\Delta \alpha ={10}^{\circ }$ and $Re=750K$. It is demonstrated that the IAG dynamic stall model implemented in Bladed consistently shows a reasonable agreement with the experimental data, similar to its positive stall prediction accuracy. This highlights that the model is applicable both for positive and negative polar regimes. The mathematical foundation presented in Section 2.2 covers the negative stall implementation strategy and all assumptions made, allowing a full data reproduction for follow up future studies.

3.8. Sensitivity to airfoil surface roughness

Investigations carried out in preceding sections were done for the airfoil equipped with a surface roughness element (leading edge grit type) which promotes earlier flow transition. This stimulates soiled/rough conditions in real wind turbine operation due to dirt or insects. For new turbines or when the blades are cleaned, the turbine performance is usually higher because the aerodynamic characteristic of the airfoil is better. This is represented by "clean" airfoil measured data in Figure 12. It can be seen that the maximum lift value is higher and the stall angle of attack is slightly delayed.

To assess the performance of Bladed dynamic stall modeling under clean conditions, studies were carried out by simulating the S809 airfoil data without the surface roughness effects. The results are presented in Figure 13. Three mean angles of attack were considered in the simulations, having $k=0.064$, $\Delta \alpha ={10}^{\circ }$ and $Re=750K$. Both the IAG model and the BL model show reasonable agreement with the experimental data for $\overline{\alpha }=8^{\circ }$ and $\overline{\alpha }={14}^{\circ }$. The Øye model underestimates the hysteresis effects for both mean angles of attack. When the airfoil is operating at a deep stall condition at $\overline{\alpha }={20}^{\circ }$, both the BL model and the Øye model fail to predict the hysteresis. At this condition, only the IAG model is able to reconstruct the dynamic polar reasonably well.

3.9. Sensitivity to Reynolds number

Finally, the discussions are extended to cover the performance of Bladed dynamic stall modeling for three different Reynolds numbers ($Re=750K$, $Re=1000K$ and $Re=1250K$) in this section. The results are depicted in Figure 14. Interestingly, the effects of dynamic stall seem to scale down with increasing Reynolds number for the S809 airfoil, although the higher harmonic effects are still present for all considered Reynolds number. This effect is probably caused by the improved flow momentum by the increased Reynolds number. As a consequence, separation is delayed and the dynamic stall vortex is weakened. It can be seen from Figure 14 that the newly implemented IAG model is demonstrated to work and scaled well when the Reynolds number is increased.

4. Conclusions and Outlook

Comprehensive studies have been performed in this paper to reformulate the IAG dynamic stall model into a state-space representation and to validate the results against experimental data. The implementation was carried out in the wind turbine design tool Bladed. The validations were done for various scenarios, including the effects of unsteady parameters, airfoil thickness, airfoil surface roughness and Reynolds number. The following aspects can be derived from the present paper:

• The state-space representation of the IAG model has been successfully formulated.
• All the governing equations and assumptions are presented in a consistent manner to allow replication for future studies.
• The newly implemented IAG model clearly shows superior performance compared to the standard incompressible Beddoes-Leishman model and the Øye model in Bladed.
• The IAG model demonstrates a good accuracy against experimental data under various unsteady parameters, including the effects of mean angle of attack, reduced frequency and excitation amplitude.
• The IAG model results are able to be generalized to airfoils having different relative thicknesses.
• The performance of the IAG model is well validated for both clean airfoil and airfoil under the effects of surface roughness.
• The IAG model agrees well with the measurement data for three different values of the Reynolds number.
• The incompressible BL model prediction is deemed sufficient for small to moderate angle of attack cases.

The implementation of the IAG model in Bladed will allow a more accurate calculations of future wind turbine load assessments. Investigations on the effects of the IAG dynamic stall model on the design load cases are highly suggested for future studies. Furthermore, dynamic stall models shall be validated for extremely high angle of attack, far beyond the static stall angle. This becomes important for wind turbine aerodynamics because the angle of attack range can reach [-180°,+180°]. It is recommended that follow up works to be aimed at this aspect.