Prerpints.org logo
Preprint
Article

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

Altmetrics

Downloads

269

Views

330

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

27 March 2023

Posted:

28 March 2023

You are already at the latest version

Alerts
Abstract
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, and impacts significantly the turbine aeroelastic stability. The employed dynamic stall model needs to be accurate and robust for a wide range of airfoils and range of angle of attack. The present studies are intended to demonstrate the performance of a recently implemented "IAG dynamic stall" model in a wind turbine design tool Bladed. The model is transformed from the indicial type of formulation into a state-space representation. The new model is validated against measurement data and other dynamic stall models in Bladed for various flow conditions and airfoils. It is demonstrated that the new model is able to reproduce the measured dynamic polar accurately without airfoil specific parameter calibration and has a superior performance compared to other models in Bladed.
Keywords: 
Subject: Engineering  -   Other

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 ( α ) 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 α n , there is a time delay between the actual effective angle of attack ( α e n ) seen by the structure. This aspect can be represented as:
α e n = α n − X n − Y n
with n being the current sample time (current value of the actual angle o attack). The deficiency functions are described by:
X n = X n − 1 exp − b 1 β 2 Δ s + A 1 Δ α n exp − b 1 β 2 Δ s / 2
Y n = Y n − 1 exp − b 2 β 2 Δ s + A 2 Δ α n exp − b 2 β 2 Δ s / 2
Note that variable time is represented by a dimensionless parameter s = 2 V t / 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, Δ α n and Δ s are given by:
Δ α n = α n + 1 − α n
Δ s = s n − s n − 1
Δ t = t n − t n − 1 .
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 β = 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 ,
which combines the circulatory and attached flow components. The circulatory normal force can be obtained using
C N n C = d C N d α ( α e n − α 0 I N V )
with α 0 I N V being the angle of attack for zero inviscid normal force and d C N / d α 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 α 0 I N V . 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 = 4 K α T I M Δ α n Δ t − D n .
with T I = M c / V . D n represents the deficiency function and is given by
D n = D n − 1 exp − Δ t K α T I + Δ α n − Δ α n − 1 Δ t exp − Δ t 2 K α T I ,
with K α 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 = d C N d α 1 + f 2 n 2 2 ( α e n − α 0 V I S C ) + C N n I .
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 ( α 0 V I S C ) 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 ,
which depends on two parameters f n and D f n . The first parameter is defined by an inverse Kirchhoff equation as:
f n = 2 C N n V I S C d C N d α ( α f n − α 0 V I S C ) − 1.0 2 ,
and can be computed based on the pressure lag response of the angle of attack:
α f n = α 0 I N V + C N n P 1 d C N / d α .
with
C N n P 1 = C N n P − D p n ,
and D p n being represented as:
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 .
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 − Δ s T f + f n − f n − 1 exp − Δ s 2 T f .
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 = 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 v l C N n − 1 V exp − Δ s T v ; otherwise
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 1 − 1 4 1 + f 2 n 2
Variables T v and T v l represent the vortex decay and vortex travel time constants, respectively. The non-dimensional vortex time ( τ v n ) itself is important for determining the vortex lift influence and can be calculated using
τ v n = τ v n − 1 + 0.45 Δ t c V ; if C N n P 1 > C N C R I T 0 ; if C N n P 1 < C N C R I T and Δ α n ≥ 0 τ v n − 1 ; otherwise
C N C R I T 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 α n C R I T . The magnitude of C N C R I T can be computed as:
C N C R I T = d C N d α ( α n C R I T − α 0 I N V ) .
The total normal force contribution can be obtained by computing
C N n D = C N n f + C N n V .
In the IAG model, the tangential (chordwise) force is simply obtained from the static polar data at the time-lagged angle of attack α f n by:
C T n D = C T V I S C ( α f n ) .
Finally, the lift and drag forces can be computed as:
C L n D = C N n D cos α n − C T n D sin α n
C D n D = C N n D sin α n + C T n D cos α n
I n the IAG model, a limiter is applied to the drag force to correct the hysteresis effects. This is done through a parameter ζ n which is defined as:
ζ n = 1 π d C N d α 1 + f n 2 2 .
Drag hysteresis is observed to occur when ζ n ≲ ζ v , with ζ v being a constant with a value of 0.76. The adopted correction reads:
C D n D = 1.2 C D n V I S C ; if C D n D > 1.2 C D n V I S C and C N n P − C N n − 1 P ≥ 0.0 C D n V I S C ; if C N n P − C N n − 1 P < 0.0 C D n D ; otherwise
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
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 α f n , that reads
C M n f = C M V I S C ( α 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. The vortex contribution of the dynamic moment coefficient can be formulated as
C M n V = − C P v n C N n V .
This depends on the idealized variation of the center of pressure and can be estimated as
C P v n = K v 1 − cos π τ v T v l .
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 = 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 v l 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

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:
c 2 V x ˙ 1 ( t ) = b 1 A 1 α n 3 / 4 − b 1 + c V ˙ V 2 x 1 ( t )
c 2 V x ˙ 2 ( t ) = b 2 A 2 α n 3 / 4 − b 2 + c V ˙ V 2 x 2 ( t )
Here, x 1 , 2 ( t ) and x ˙ 1 , 2 ( t ) 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 c V ˙ V 2 is added to include the effects for varying wind speed in wind turbine cases. The effective angle of attack is calculated using:
α e n = α n 3 / 4 ( 1 − A 1 − A 2 ) + x 1 ( t ) + x 2 ( t )
The circulatory normal force can be obtained using
C N n C = d C N d α sin ( α e n − α 0 V I S C ) .
Note that Equation (36) looks different compared to Equation (8). One of it is the usage of α 0 V I S C instead of α 0 I N V . 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 = 4 K α c V α ˙ n .
Finally, the attached flow normal force coefficient is obtained by:
C N n P = C N n C + C N n I .
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:
c 2 V T p x ˙ 3 ( t ) = − x 3 ( t ) + C N n P .
Variable x 3 ( t ) is the solution for the pressure-lagged normal force or C N n P 1 in Equation (15). This parameter is used to calculate the delayed angle of attack
α f n = α 0 V I S C + x 3 ( t ) d C N / d α .
By using the inverse of the Kirchoff equation, the position of separation point at α f n can be determined by:
f n = 2 C N n V I S C d C N d α sin ( α f n − α 0 V I S C ) − 1.0 2 .
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:
c 2 V T f x ˙ 4 ( t ) = − x 4 ( t ) + f n ,
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 = d C N d α 1 + x 4 ( t ) 2 2 sin ( α e n − α 0 V I S C ) + C N n I .
Note that x 4 ( t ) 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:
c 2 V T v x ˙ 5 ( t ) = − x 5 ( t ) + c 2 V T v C ˙ V n ,
with x 5 ( t ) 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:
C ˙ V n = d C N d α α ˙ e n 1 − 1 4 1 + x 4 ( t ) 2 − 1 4 C N n C 1 + x 4 ( t ) x ˙ 4 ( t ) x 4 ( t ) .
Equation (45) is only applied when there is a dynamic vortex state being active, otherwise C ˙ V n = 0 . Vortex state itself is defined when x 3 ( t ) is greater/smaller than the critical normal force:
C N C R I T = C N M A X ; if Positive Stall C N M I N ; if Negative Stall
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 C ˙ V n is that it has to be during the upstroke motion. This implies that:
Vortex state = True ; if Upstroke & x 3 ( t ) > C N C R I T & 0 < τ v n < T v l False ; otherwise
The upstroke state can be determined by Δ α n , but this will be different when the case being considered is positive stall ( Δ α n > 0 ) or negative stall ( Δ α n < 0 ). A positive stall case can be identified when x 3 ( t ) ≥ 0 , while a negative stall case may be detected when x 3 ( t ) < 0 . For the backwinded case (when ∥ α ∥ is greater than 90°), the definition of upstroke is flipped. Now, the non-dimensional vortex time may be calculated using
τ v n = τ v n − 1 + 0.45 Δ t c V ; if Vortex state τ v n − 1 exp − Δ t c V ; otherwise
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 .
On top of that, the tangential (chordwise) force is obtained directly from the static polar data at the time-lagged angle of attack α f n by:
C T n D = C T V I S C ( α f n ) .
In the present model, only lift coefficient is calculated from C N n D and C T n D transformation as:
C L n D = C N n D cos α n − C T n D sin α n .
The calculation for drag coefficient is done differently by applying the following relation:
C D n D = C D n V I S C + α n − α e n C N n C + C D n V I S C − C D 0 V I S C 1 − x 4 ( t ) 2 2 − 1 − f V I S C 2 2 + x 5 ( t ) sin α n
In Equation (52), variables C D n V I S C , C D 0 V I S C and f V I S C 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 ∥ α ∥ < 30 ∘ 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 .
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 = − π c α ˙ n 4 V .
For backwinded airfoil orientation (when ∥ α ∥ 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 ζ 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 α = 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 α 0 V I S C and α = 7 ∘ 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 π and 2.5 π 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 ( δ / c ), but are also extended for the S801 ( δ / c = 13.5 % ) and S814 ( δ / 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 ( R e ) 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 α ¯ = 20 ∘ 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 ( Δ α = 5 . 5 ∘ and Δ α = 5 . 5 ∘ ). The highest available mean angle of attack ( α ¯ = 20 ∘ ) 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 Δ α = 5 . 5 ∘ , 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.
α N e g = − α + 2 α C L = 0
C L N e g = − C L
C D N e g = C D
C M N e g = − C M + 2 C M C L = 0
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 , Δ α = 10 ∘ and R e = 750 K . 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 , Δ α = 10 ∘ and R e = 750 K . Both the IAG model and the BL model show reasonable agreement with the experimental data for α ¯ = 8 ∘ and α ¯ = 14 ∘ . The Øye model underestimates the hysteresis effects for both mean angles of attack. When the airfoil is operating at a deep stall condition at α ¯ = 20 ∘ , 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 ( R e = 750 K , R e = 1000 K and R e = 1250 K ) 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.

Author Contributions

Galih Bangga: Conceptualization, data curation, formal analysis, investigation, methodology, resources, validation, writing - original draft, writing - review and editing. Steven Parkinson: Writing - review and editing. William Collier: Writing - review and editing.

Funding

This research was funded by DNV research project grant "IEA Wind Task 47" and "IAG Dynamic Stall Model".

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data can be made available by contacting the corresponding author.

Acknowledgments

The authors would like to thank all Bladed team members who provided assistance in the testing process of the IAG dynamic stall model, especially to Dr. Ali Bakhshandeh Rostami, Mohit Mukhija, Sigurd Jensen, Hassan Moharram and Clint Kallada. The measurement data sets provided by the Ohio State University are highly appreciated.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Martin, J.; Empey, R.; McCroskey, W.; Caradonna, F. An experimental analysis of dynamic stall on an oscillating airfoil. Journal of the American Helicopter Society 1974, 19, 26–32. [Google Scholar] [CrossRef]
  2. Carr, L.W.; McAlister, K.W.; McCroskey, W.J. Analysis of the development of dynamic stall based on oscillating airfoil experiments. Technical report, NASA TN D-8382, National Aeronautics and Space Administration, United States, 1977.
  3. McAlister, K.W.; Carr, L.W.; McCroskey, W.J. Dynamic stall experiments on the NACA 0012 airfoil. Technical report, NASA Technical Paper 1100, National Aeronautics and Space Administration, United States, 1978.
  4. Ramsay, R.; Hoffman, M.; Gregorek, G. Effects of grit roughness and pitch oscillations on the S801 airfoil. Technical report, National Renewable Energy Lab., Golden, CO (United States), 1996.
  5. Hoffman, M.; Ramsay, R.; Gregorek, G. Effects of grit roughness and pitch oscillations on the NACA 4415 airfoil. Technical report, National Renewable Energy Lab., Golden, CO (United States), 1996.
  6. Ramsay, R.; Hoffman, M.; Gregorek, G. Effects of grit roughness and pitch oscillations on the S809 airfoil. Technical report, National Renewable Energy Lab., Golden, CO (United States), 1995.
  7. Janiszewska, J.; Ramsay, R.; Hoffman, M.; Gregorek, G. Effects of grit roughness and pitch oscillations on the S814 airfoil. Technical report, National Renewable Energy Lab., Golden, CO (United States), 1996.
  8. Bangga, G.; Hutomo, G.; Wiranegara, R.; Sasongko, H. Numerical study on a single bladed vertical axis wind turbine under dynamic stall. Journal of Mechanical Science and Technology 2017, 31, 261–267. [Google Scholar] [CrossRef]
  9. Bangga, G. Numerical studies on dynamic stall characteristics of a wind turbine airfoil. Journal of Mechanical Science and Technology 2019, 33, 1257–1262. [Google Scholar] [CrossRef]
  10. Bangga, G.; Hutani, S.; Heramarwan, H. The Effects of Airfoil Thickness on Dynamic Stall Characteristics of High-Solidity Vertical Axis Wind Turbines. Advanced Theory and Simulations 2021, 4, 2000204. [Google Scholar] [CrossRef]
  11. Leishman, J.G.; Beddoes, T. A Semi-Empirical model for dynamic stall. Journal of the American Helicopter society 1989, 34, 3–17. [Google Scholar]
  12. Hansen, M.H.; Gaunaa, M.; Madsen, H.A. A Beddoes-Leishman type dynamic stall model in state-space and indicial formulations. Technical report, Risø-R-1354, Risø National Laboratory, Denmark, 2004.
  13. Bangga, G.; Lutz, T.; Arnold, M. An improved second-order dynamic stall model for wind turbine airfoils. Wind Energy Science 2020, 5, 1037–1058. [Google Scholar] [CrossRef]
  14. Larsen, J.W.; Nielsen, S.R.; Krenk, S. Dynamic stall model for wind turbine airfoils. Journal of Fluids and Structures 2007, 23, 959–982. [Google Scholar] [CrossRef]
  15. Gupta, S.; Leishman, J.G. Dynamic stall modelling of the S809 aerofoil and comparison with experiments. Wind Energy 2006, 9, 521–547. [Google Scholar] [CrossRef]
  16. Elgammi, M.; Sant, T. A Modified Beddoes–Leishman Model for Unsteady Aerodynamic Blade Load Computations on Wind Turbine Blades. Journal of Solar Energy Engineering 2016, 138, 051009. [Google Scholar] [CrossRef]
  17. Wang, Q.; Zhao, Q. Modification of Leishman–Beddoes model incorporating with a new trailing-edge vortex model. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering 2015, 229, 1606–1615. [Google Scholar] [CrossRef]
  18. Sheng, W.; Galbraith, R.M.; Coton, F. A new stall-onset criterion for low speed dynamic-stall. Journal of solar energy engineering 2006, 128, 461–471. [Google Scholar] [CrossRef]
  19. Galbraith, R. Return from Airfoil Stall During Ramp-Down Pitching Motions. Journal of Aircraft 2007, 44, 1856–1864. [Google Scholar]
  20. Sheng, W.; Galbraith, R.; Coton, F. A modified dynamic stall model for low Mach numbers. Journal of Solar Energy Engineering 2008, 130. [Google Scholar] [CrossRef]
  21. Øye, S. Dynamic stall simulated as time lag of separation. In Proceedings of the Proceedings of the 4th IEA Symposium on the Aerodynamics of Wind Turbines. Rome, Italy, 1991.
  22. Tran, C.; 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 Proceedings of the 6th European Rotorcraft Forum, 1980.
  23. Tarzanin, F. Prediction of control loads due to blade stall. Journal of the American Helicopter Society 1972, 17, 33–46. [Google Scholar] [CrossRef]
  24. Snel, H. Heuristic modelling of dynamic stall characteristic. In Proceedings of the Proceedings of the European Wind Energy Conference, 1997.
  25. Adema, N.; Kloosterman, M.; Schepers, G. Development of a Second Order Dynamic Stall Model. Wind Energy Science Discussions 2019, 2019, 1–18. [Google Scholar] [CrossRef]
  26. UK, D.S. Bladed Theory Manual 4.14; DNV Services UK, 2022.
  27. Beddoes, T. Practical computation of unsteady lift. In Proceedings of the 8th European Rotorcraft Forum, 1982.
  28. Leishman, J. Validation of approximate indicial aerodynamic functions for two-dimensional subsonic flow. Journal of Aircraft 1988, 25, 914–922. [Google Scholar] [CrossRef]
  29. Garbaruk, A.; Shur, M.; Strelets, M.; Travin, A. NACA0021 at 60 deg. incidence. Notes on Numerical Fluid Mechanics and Multidisiplinary Design 2009, 103, 127–139. [Google Scholar]
  30. Bangga, G. Comparison of blade element method and CFD simulations of a 10 MW wind turbine. Fluids 2018, 3, 73. [Google Scholar] [CrossRef]
  31. Bangga, G.; Lutz, T. Aerodynamic modeling of wind turbine loads exposed to turbulent inflow and validation with experimental data. Energy 2021, 223, 120076. [Google Scholar] [CrossRef]
Preprints 70226 g001
Figure 1. Comparison between inviscid normal force reconstruction using a linear function, a sinusoidal function and the inviscid XFOIL data for NACA 0021.
Figure 1. Comparison between inviscid normal force reconstruction using a linear function, a sinusoidal function and the inviscid XFOIL data for NACA 0021.
Preprints 70226 g001
Preprints 70226 g002
Figure 2. Determination of the normal force gradient required for dynamic stall calculations.
Figure 2. Determination of the normal force gradient required for dynamic stall calculations.
Preprints 70226 g002
Preprints 70226 g003
Figure 3. Comparison of the measured angle of attack data points with a cubic fitted time series for the S809 airfoil at α ¯ = 20 ∘ , k = 0.064 , Δ α = 10 ∘ and R e = 750 K . Variable k represents the reduced frequency ω c / ( 2 V ) .
Figure 3. Comparison of the measured angle of attack data points with a cubic fitted time series for the S809 airfoil at α ¯ = 20 ∘ , k = 0.064 , Δ α = 10 ∘ and R e = 750 K . Variable k represents the reduced frequency ω c / ( 2 V ) .
Preprints 70226 g003
Preprints 70226 g004
Figure 4. Consistency between the IAG model in state-space and indicial formulations. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , Δ α = 10 ∘ , k = 0.064 and R e = 750 K . The critical normal force in B-GO was provided as the inviscid normal force.
Figure 4. Consistency between the IAG model in state-space and indicial formulations. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , Δ α = 10 ∘ , k = 0.064 and R e = 750 K . The critical normal force in B-GO was provided as the inviscid normal force.
Preprints 70226 g004
Preprints 70226 g005
Figure 5. Consistency between the IAG model in state-space and indicial formulations. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , Δ α = 10 ∘ , k = 0.064 and R e = 750 K . The critical normal force in B-GO was provided as the viscous normal force obtained directly from the polar data.
Figure 5. Consistency between the IAG model in state-space and indicial formulations. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , Δ α = 10 ∘ , k = 0.064 and R e = 750 K . The critical normal force in B-GO was provided as the viscous normal force obtained directly from the polar data.
Preprints 70226 g005
Preprints 70226 g006
Figure 6. Dynamic polar reconstruction using Bladed for three different mean angles of attack and the comparison against measurement data. The dynamic polar was simulated under the following conditions: Δ α = 10 ∘ , k = 0.064 and R e = 750 K .
Figure 6. Dynamic polar reconstruction using Bladed for three different mean angles of attack and the comparison against measurement data. The dynamic polar was simulated under the following conditions: Δ α = 10 ∘ , k = 0.064 and R e = 750 K .
Preprints 70226 g006
Preprints 70226 g007
Figure 7. Dynamic polar reconstruction using Bladed for three different reduced frequencies and the comparison against measurement data. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , Δ α = 10 ∘ and R e = 750 K .
Figure 7. Dynamic polar reconstruction using Bladed for three different reduced frequencies and the comparison against measurement data. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , Δ α = 10 ∘ and R e = 750 K .
Preprints 70226 g007
Preprints 70226 g008
Figure 8. Dynamic polar reconstruction using Bladed for two different excitation amplitudes and the comparison against measurement data. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , k = 0.064 and R e = 750 K .
Figure 8. Dynamic polar reconstruction using Bladed for two different excitation amplitudes and the comparison against measurement data. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , k = 0.064 and R e = 750 K .
Preprints 70226 g008
Preprints 70226 g009
Figure 9. Dynamic polar reconstruction using Bladed for three airfoils having different relative thicknesses and the comparison against measurement data. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , k = 0.064 , Δ α = 10 ∘ and R e = 750 K . The considered airfoils are: S801 ( δ / c = 13.5 % ), S809 ( δ / c = 21 % ) and S814 ( δ / c = 24 % ). Note that the exact value of the reduced frequency is different for each airfoil, but they are close to each other.
Figure 9. Dynamic polar reconstruction using Bladed for three airfoils having different relative thicknesses and the comparison against measurement data. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , k = 0.064 , Δ α = 10 ∘ and R e = 750 K . The considered airfoils are: S801 ( δ / c = 13.5 % ), S809 ( δ / c = 21 % ) and S814 ( δ / c = 24 % ). Note that the exact value of the reduced frequency is different for each airfoil, but they are close to each other.
Preprints 70226 g009
Preprints 70226 g010
Figure 10. Negative polar reconstruction based on the available positive data of the S809 airfoil.
Figure 10. Negative polar reconstruction based on the available positive data of the S809 airfoil.
Preprints 70226 g010
Preprints 70226 g011
Figure 11. Dynamic polar reconstruction using Bladed for negative stall cases and the comparison against measurement data. The dynamic polar was simulated under the following conditions: k = 0.064 , Δ α = 10 ∘ and R e = 750 K .
Figure 11. Dynamic polar reconstruction using Bladed for negative stall cases and the comparison against measurement data. The dynamic polar was simulated under the following conditions: k = 0.064 , Δ α = 10 ∘ and R e = 750 K .
Preprints 70226 g011
Preprints 70226 g012
Figure 12. Comparison between the clean and soiled static polar data of the S809 airfoil.
Figure 12. Comparison between the clean and soiled static polar data of the S809 airfoil.
Preprints 70226 g012
Preprints 70226 g013
Figure 13. Dynamic polar reconstruction using Bladed for clean airfoil cases and the comparison against measurement data. The dynamic polar was simulated under the following conditions: k = 0.064 , Δ α = 10 ∘ and R e = 750 K .
Figure 13. Dynamic polar reconstruction using Bladed for clean airfoil cases and the comparison against measurement data. The dynamic polar was simulated under the following conditions: k = 0.064 , Δ α = 10 ∘ and R e = 750 K .
Preprints 70226 g013
Preprints 70226 g014
Figure 14. Dynamic polar reconstruction using Bladed for three different Reynolds numbers and the comparison against measurement data. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , k = 0.064 , Δ α = 10 ∘ . The considered Reynolds numbers are: 750K, 1000K and 1250K.
Figure 14. Dynamic polar reconstruction using Bladed for three different Reynolds numbers and the comparison against measurement data. The dynamic polar was simulated under the following conditions: α ¯ = 20 ∘ , k = 0.064 , Δ α = 10 ∘ . The considered Reynolds numbers are: 750K, 1000K and 1250K.
Preprints 70226 g014
Table 1. Constants applied for the IAG models.
Table 1. Constants applied for the IAG models.
Model A 1 A 2 b 1 b 2 K α T p T f T v T v l K v K f C T m U T m D ζ s
Indicial IAG 0.3 0.7 0.7 0.53 0.75 1.7 3.0 6.0 6.0 0.2 0.1 1.5 1.5 0.76
State-Space IAG 0.3 0.7 0.7 0.53 0.75 1.7 3.0 6.0 6.0 0.2 N/A N/A N/A N/A
Table 2. Critical angle of attack ( α n C R I T ) applied for the IAG models.
Table 2. Critical angle of attack ( α n C R I T ) applied for the IAG models.
Model S801 S809 S814
Indicial IAG 15.1° 14.1° 10°
State-Space IAG Automatic Automatic Automatic Automatic
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated