the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
An improved secondorder dynamic stall model for wind turbine airfoils
Thorsten Lutz
Matthias Arnold
Robust and accurate dynamic stall modeling remains one of the most difficult tasks in wind turbine load calculations despite its long research effort in the past. In the present paper, a new secondorder dynamic stall model is developed with the main aim to model the higher harmonics of the vortex shedding while retaining its robustness for various flow conditions and airfoils. Comprehensive investigations and tests are performed at various flow conditions. The occurring physical characteristics for each case are discussed and evaluated in the present studies. The improved model is also tested on four different airfoils with different relative thicknesses. The validation against measurement data demonstrates that the improved model is able to reproduce the dynamic polar accurately without airfoilspecific parameter calibration for each investigated flow condition and airfoil. This can deliver further benefits to industrial applications where experimental/reference data for calibrating the model are not always available.
 Article
(3640 KB)  Fulltext XML
 BibTeX
 EndNote
An accurate prediction of wind turbine blade loads is influenced by many parameters including 3D and unsteady effects. The first mainly occurs in the root and tip areas of the blade due to radial flow and induced velocity influences, respectively (Bangga, 2018). The latter can occur due to variation in the inflow conditions caused by yaw misalignment, wind turbulence, shear and gusts, tower shadow, and aeroelastic effects of the blade. The abovementioned phenomena may result in dynamic stall (DS). Experimental studies (Martin et al., 1974; Carr et al., 1977; McAlister et al., 1978) showed that the aerodynamic forces can differ significantly in comparison to the static condition. DS is often initiated by the generation of a leading edge vortex (LEV), which increases the positive circulation effect on the airfoil suction side, causing delayed stall (Bangga, 2019). This intense leading edge vortex is convected downstream along the airfoil towards the trailing edge. At the same time, the lift force increases significantly, and the pitching moment becomes more negative compared to the static values. A significant drag increase is observed at large angles of attack. An example is shown in Fig. 1. Afterwards, a trailing edge vortex (TEV) with opposite rotational direction than LEV is formed, which pushes the leading edge vortex towards the wake area. This onset may result in a significant drop of the lift coefficient (C_{L}) and can be dangerous for the blade structure itself, although dynamic stall also generally enhances aerodynamic damping.
To model the behavior of the airfoil under these situations, semiempirical models can be used. The models are known to produce reasonable results without any notable increase in computational effort. Despite that, these models usually cannot reproduce higher harmonics of the load fluctuations. Furthermore, the applied constants shall be adjusted according to the flow conditions and airfoils. Leishman and Beddoes (LB) (Leishman and Beddoes, 1989) have developed a model for dynamic stall combining the flow delay effects of attached flow with an approximate representation of the development and effect of separation (Larsen et al., 2007). This model was developed for helicopter applications and therefore includes a fairly elaborate representation of the nonstationary attached flow depending on the Mach number and a rather complex structure of the equations representing the time delays (Larsen et al., 2007). Hansen et al. (2004) simplified the model for wind turbine applications by removing the consideration of compressibility effects and the leading edge separation. The latter was argued because the relative thickness of wind turbine airfoil is typically no less than 15 %. This model was called the Risø model in Larsen et al. (2007). Examples of the other models are given by Øye (Øye, 1991), Tran and Petot (Tran and Petot, 1980) (ONERA model), and Tarzanin (Boeing Vertol model) (Tarzanin, 1972). To better model the vortex shedding characteristics at large angles of attack, secondorder dynamic stall models were introduced. An example of this model was given by Snel (1997) which makes use of the difference between the inviscid and viscous static polar data as a main forcing term for the dynamic polar reconstruction, in contrast to the LB model that uses the changes of the angle of attack over the time. An improved version of the Snel model was proposed recently by Adema et al. (2020) to cover for the increased shedding effects in the downstroke phase. All abovementioned models employ the static polar data and dynamic flow parameters as the input needed for the dynamic polar reconstruction. Then, the models compute the dynamic force difference required for the reconstruction process.
Although many studies have been dedicated to dynamic stall modeling (Gupta and Leishman, 2006; Larsen et al., 2007; Adema et al., 2020; Elgammi and Sant, 2016; Wang and Zhao, 2015; Sheng et al., 2006; Galbraith, 2007; Sheng et al., 2008), engineering calculations in the industry are still relying on the very basic classical dynamic stall models such as the Leishman–Beddoes and Snel models. The reason is the simplicity to tune in the models for different airfoils and for different flow conditions. Therefore, one major key for a model to be used in industrial applications is robustness of the model itself; i.e., the model is easy to apply with a small number of welldefined user parameters. One of the purposes of this paper is to document widely used dynamic stall models in research and industries. These include the firstorder LB model and the secondorder Snel model. A very recently improved Snel model according to Adema et al. (2020) will also be evaluated. The mathematical formulations of these models will be presented in this report. Weaknesses of existing dynamic stall modeling shall be identified, and possible corrections to those limitations will be described. Finally, a new secondorder dynamic stall modeling will be proposed that is able to model not only the secondorder lift and drag forces, but also the pitching moment along with calculation examples in comparison to experimental data for different airfoils and flow conditions.
The paper is organized as follows. Section 2 describes the mathematical formulation of existing dynamic stall models and the new model developed in this work. Then, in Sect. 3 assessments are carried out on the performance of the secondorder dynamic stall models in comparison with measurement data. The new model is further tested at various flow conditions and to examine its robustness on four different airfoils without further calibrating the constants. Finally, all results will be concluded in Sect. 4.
In this section the mathematical formulations of each model are described in detail. The reasons are mainly to provide information on how each model was employed and to gain deeper insights for further developing the new model. Note that each existing model was developed by different authors; thus different symbols and formulation methods were adopted in those publications (Beddoes, 1982; Leishman, 1988; Leishman and Beddoes, 1989; Snel, 1997; Adema et al., 2020). In this paper, all models are described in a consistent way for clarity and for an easier interpretation and implementation process.
2.1 Leishman–Beddoes model
The original Leishman–Beddoes model is composed of three main contributions representing various flow regimes: (1) unsteady attached flow, (2) unsteady separated flow and (3) dynamic stall. The present section will elaborate the mathematical description and its physical interpretation of each module. Figure 2 illustrates several main parameters needed for modeling the dynamic stall characteristics.
2.1.1 Unsteady attached flow
In this module, the unsteady aerodynamic response of the loads is represented by the time delay effects. The indical formulas were constructed based on the work of Beddoes (1982) and have been refined by Leishman (1988). The loads are assumed to originate from two main sources: one for the initial noncirculatory loading from the piston theory and another for the circulatory loading which builds up quickly to the steadystate value (Leishman and Beddoes, 1989). In the formulation, the relative distance traveled by the airfoil in terms of semichords is represented by $s=\mathrm{2}Vt/c$, which can also be used to describe the nondimensional time. Note that V, t and c are freestream wind speed, time and chord length, respectively. For a continuously changing angle of attack α_{n}, the effective angle of attack (${\mathit{\alpha}}_{{\mathrm{e}}_{n}}$) can be represented as
where n is the current sample time. The last two terms describe the deficiency functions that are given by
where
In these equations, b_{1}, b_{2}, A_{1} and A_{2} are constants. The variable β represents the compressibility effects and is formulated as $\mathit{\beta}=\sqrt{\mathrm{1}{M}^{\mathrm{2}}}$. Because information about the previous cycle is needed in the formulations, initializations are required. The solution needs to develop for a certain time until convergence of the resulting unsteady loads is obtained.
The circulatory normal force due to an accumulating series of step inputs in angle of attack can be obtained using
The variable ${\mathit{\alpha}}_{{\mathrm{0}}_{\mathrm{INV}}}$ is the angle of attack for zero inviscid normal force. The original formulation of the model disregarded the use of ${\mathit{\alpha}}_{\mathrm{0}}^{\mathrm{INV}}$. However, this term is important when the airfoil has a finite camber. This has been pointed out as well by Hansen et al. (2004).
The noncirculatory (impulsive) normal force is obtained by
where T_{I} is given by ${T}_{\mathrm{I}}=Mc/V$. The deficiency function D_{n} is given by
and $\mathrm{\Delta}t={t}_{n}{t}_{n\mathrm{1}}$.
The total normal force coefficient under attached flow conditions is given by the sum of circulatory and noncirculatory components as
2.1.2 Unsteady separated flow
Leishman and Beddoes (1989) stated that the onset of leading edge separation is the most important aspect in dynamic stall modeling. The condition under which leading edge stall occurs is controlled by a critical leading edge pressure coefficient that is linked into the formulation by defining a lagged normal force coefficient ${C}_{{\mathrm{N}}_{n}}^{\mathrm{P}\mathrm{1}}$ as
where ${D}_{{\mathrm{p}}_{n}}$ is given by
It has been investigated by Leishman and Beddoes (1989) that the calibration time constant T_{p} is largely independent of the airfoil shape. The substitute value of the effective angle of attack incorporating the leading edge pressure lag response may be obtained using
In most airfoil shapes, the progressive trailing edge separation causes loss of circulation and introduces nonlinear effects on the lift, drag and pitching moment, especially on cambered airfoils. This is even more important for wind turbine airfoils because the relative thickness is large. To derive a correlation between the normal force coefficient and the separation location (f_{n}), the relation based on the flat plate from Kirchhoff and Helmholtz can be used, which reads
The location of the separation point is usually obtained by a curvefitting procedure in literature. For example, Leishman and Beddoes (1989) proposed the following correlation:
The coefficients S_{1} and S_{2} define the static stall characteristic while α_{1} defines the static stall angle. The derivation was based on the NACA 0012, HH02 and SC1095 airfoils that have a single break point of the static lift force coefficient. Gupta and Leishman (2006) proposed the formulation for the S809 airfoil as
that has two break points (α_{1} and α_{2}) of the static lift force coefficient, where c_{1}, c_{2}, c_{3}, a_{1}, a_{2} and a_{3} are constants.
The additional effects of the unsteady boundary layer response may be represented by application of a firstorder lag to the value of f_{n} to produce the final value for the unsteady trailing edge separation point ${f}_{{\mathrm{2}}_{n}}$ (Leishman and Beddoes, 1989). This can be represented as
where ${D}_{{\mathrm{f}}_{n}}$ is given by
and T_{f} is a constant. Then, the unsteady viscous normal force coefficient for each sample time can be obtained using
The tangential component of the force can be obtained by (Leishman and Beddoes, 1989)
Note that positive ${C}_{{\mathrm{T}}_{n}}^{\mathrm{f}}$ is defined in the direction of the trailing edge while η is a constant.
According to Leishman and Beddoes (1989) and Gupta and Leishman (2006), a general expression for the pitching moment behavior cannot be obtained from Kirchhoff theory, and an alternative empirical relation must be formulated. Gupta and Leishman (2006) proposed the following formulation for the S809 airfoil
where ${C}_{{\mathrm{M}}_{\mathrm{0}}}$ defines the moment coefficient at zero normal force and K_{0} is the mean offset of the aerodynamic center from the quarter chord position. K_{1}, K_{2}, K_{3}, K_{4} and m are constants.
2.1.3 Dynamic stall
The third part of the model describes the poststall characteristics where the vortical disturbances near the leading edge become stronger. The effect of vortex shedding is given by defining the vortex lift as the difference between the linearized value of the unsteady circulatory normal force and the unsteady nonlinear normal force obtained from the Kirchhoff approximation, which reads
where K_{n} is given by
The normal force is allowed to decay, but it is updated with a new increment in the normal force based on prior forcing conditions, which can be defined as
where T_{v} and T_{vl} are the vortex decay and center of pressure travel time constants, respectively. The nondimensional vortex time is given by (Pereira et al., 2011; Elgammi and Sant, 2016)
with ${C}_{\mathrm{N}}^{\mathrm{CRIT}}$ being the inviscid critical static normal force, usually indicated by the break of the (viscous) moment polar at the critical angle of attack ${\mathit{\alpha}}_{n}^{\mathrm{CRIT}}$. This can be formulated as
The idealized variation in the center of pressure with the convection of the leading edge vortex can be modeled by
The dynamic moment coefficient can be formulated as
Therefore, the total dynamic loading on the airfoil from all modules can be written as
and by converting these forces into lift and drag, one obtains
2.1.4 Note to present implementation
In Eqs. (14) and (15), a curvefitting procedure is usually adopted in literature. In this sense, the parameters or even the formulation need to be adjusted when the airfoil is different. Therefore, in the present implementation, the separation point is derived directly from the static polar data using the inversion of Eq. (13) as
The same approach was used for example by Hansen et al. (2004). This way, the user can avoid dealing with curve fitting adjustment (which requires changes on the constants for different airfoils and flow conditions) as long as the static polar data are available.
In the original formulation, the pitching moment is also obtained by a curvefitting procedure in Eq. (20). Again, this kind of approach is not straightforward as the user needs to perform curve fitting of the polar data. In the present implementation, the moment coefficient is easily obtained from the static viscous polar data by interpolating the value at the effective angle of attack, incorporating the leading edge pressure time lag ${\mathit{\alpha}}_{{\mathrm{f}}_{n}}$, which reads
In this sense, the moment coefficient can be reconstructed easily without the need to adjust the parameters in advance, minimizing the user error.
Furthermore, to avoid discontinuity in the downstroke phase for Eq. (24), an additional condition is applied in the present implementation as
2.2 Snel secondorder model
The history of the Snel secondorder model (Snel, 1997) dates back to 1993 based on Truong's observation on dynamic lift coefficient characteristics (Truong, 1993). Truong proposed that the difference between the static and dynamic lift can be divided into two terms: the forcing frequency response and the higherfrequency dynamics of a selfexcited nature. The total dynamic response of the airfoil is formulated as
with ^{D1} and ^{D2} being the first and secondorder corrections, respectively. The first correction is modeled using an ordinary differential equation (ODE) by applying a springdamping like function as
The frequency of the firstorder corrected lift follows the frequency of the forcing term F_{1}. This term is based on the time derivative of the difference between the steady inviscid coefficient ${C}_{{\mathrm{L}}_{n}}^{\mathrm{INV}}$ and viscous lift coefficient ${C}_{{\mathrm{L}}_{n}}^{\mathrm{VISC}}$ of an airfoil ($\mathrm{\Delta}{C}_{{\mathrm{L}}_{n}}^{\mathrm{INV}}$) as
with n and dC_{L}∕dα as the current sample time and inviscid lift gradient, respectively. The time constant τ in the above equation represents the time required for the flow to travel half a chord distance as
The “stiffness” coefficient of the firstorder term $K{f}_{{\mathrm{10}}_{n}}$ can be expressed as
As shown in Faber (2018), the above equation becomes numerically unstable if $\dot{{\mathit{\alpha}}_{n}}$ is large (increasing reduced frequency above 0.1) for $\dot{{\mathit{\alpha}}_{n}}{C}_{{\mathrm{L}}_{n}}^{\mathrm{INV}}\le \mathrm{0}$. The reason is that the denominator goes to zero and then negative, causing numerical integration instability. Thus, based on pure intuition the denominator value was set to a minimum of 2.0 in Faber (2018). In the present implementation, a similar approach is adopted but the limit differs. Instead, the minimum denominator value is limited to $\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$, because it yields more physical results for several cases tested by the authors.
To incorporate the higherorder frequency dynamics, a secondorder ODE is used to describe the secondorder correction term. The general form may be written as
similar to the firstorder correction. The frequency of the higherorder dynamics is determined by the forcing term ${F}_{{\mathrm{2}}_{n}}$, defined as
It is noted that the value 0.1 as a constant was chosen according to Adema et al. (2020). This is not a fixed value and can be adjusted based on the evaluated cases as seen in literature (Adema et al., 2020; Snel, 1997; Holierhoek et al., 2013; Faber, 2018; Khan, 2018). Variable k_{s} is a constant with a typical value of 0.2. The spring coefficient is given by
and the damping coefficient as
2.3 Adema–Snel secondorder model
The recently developed model of Adema et al. (2020) improves the original Snel model (Snel, 1997) in several aspects. Instead of using the lift coefficient (C_{L}), the normal force coefficient (C_{N}) is used, similar to the LB model (Leishman and Beddoes, 1989). The total dynamic response of the airfoil is formulated as follows.
The model introduces some modifications of the original model in terms of (1) projected k_{s}, (2) the firstorder coefficient and (3) the secondorder coefficient. The mathematical formulation of the firstorder term of the model is listed as
and for the secondorder correction term as
One may notice that Eq. (56) contains τ in the second term of the righthand side (RHS). This is intended to remove the dependency of the model on the velocity as the input parameter. The other main difference with the original model is also observed in Eq. (57), where k_{s} is projected by sin α_{n}. At last, the downstroke motion of the secondorder term of Eq. (58) is modified to enable vortex shedding effects.
To sum up the characteristics of the abovediscussed stateoftheart dynamic stall models, Table 1 lists the properties of each model and in which aspects the model can be improved further.
2.4 New secondorder IAG model
The proposed IAG model is developed based on knowledge gained from four different models: Leishman–Beddoes, Snel, Adema–Snel and ONERA (Tran and Petot, 1980; Dat and Tran, 1981; Petot, 1989) models with modifications. Similar to the modern models like those from Snel (and ONERA) and its derivatives, the present model is constructed by two main terms: the firstorder and secondorder corrections. The total dynamic response of the airfoil is formulated as
with ^{D1} and ^{D2} being the first and secondorder corrections, respectively. Below the description of the modifications made for the new model will be discussed in detail.
2.4.1 Firstorder correction
Based on the Hopf bifurcation model of Truong (1993) that used the LB model as the starting point of the firstorder correction, the present model operates similarly. Despite that, the LB model is not transferred into the statespace formulation, but it is retained as the indical formulation. The model applies the superposition of the solution using a finitedifference approximation to Duhamel's integral to construct the cumulative effect on an arbitrary time history of angle of attack. The LB model described in Sect. 2.1.1 to 2.1.3 will be used with the following modifications.
In the above LB model, predictions for drag are not accurate as will be shown in Sect. 3.1. This inaccuracy lies in the determination of η in Eq. (19) for the tangential force component because drag is more sensitive to tangential force than the lift force is. Also, to maintain simplicity, parameter η is removed, and the tangential force is obtained from the static data at the timelagged angle of attack ${\mathit{\alpha}}_{{\mathrm{f}}_{n}}$ by
If one uses this formulation directly, at some point drag still becomes lower than the static drag value by a significant amount. By evaluating the experimental data for several airfoils and various flow conditions, this is not physical at small angles of attack, especially in the downstroke regime, where it usually just returns to the static value. In fact, those experimental data infer that strong drag hysteresis occurs only at high angles of attack beyond stall. Similarly, in the upstroke regime the drag value increases only slightly (approximately only 20 %). In Fig. 3, one can see that drag hysteresis occurs when ζ_{n}≲ζ_{v}, with
and ζ_{v} is a constant with a value of 0.76. Based on these observations, a simple drag limiting factor is adopted when ζ_{n}≥0.76 as
Note that for the purpose of numerical implementation, it is always recommended in practice to adopt relaxation to avoid any discontinuity which may present in the above formulation. Furthermore, the value of ζ_{v} may be chosen differently for different airfoils depending on the vortex lift influence on drag. Further studies on this aspect are encouraged. The effects of these modifications are displayed in Fig. 4.
It will also be shown in Sect. 3.1 that predicting moment coefficient directly from the static polar data by means of the timelagged angle of attack has its drawback in the correct damping effect calculation. One may obtain better results by using the “fitting function” as in Eq. (20) instead of using Eq. (34). However, this limits the usability for different airfoils, since the fitting has to be done again for each individual airfoil. For wind turbine simulations, this is fairly impractical because a wind turbine blade is usually constructed by several different airfoils, not to mention the interpolated shapes in between each airfoil position. To solve for this issue, a relatively simple approach is introduced by applying a time delay to the circulatory moment response as
where
with ${K}_{\mathrm{f}}^{\mathrm{C}}$, ${T}_{\mathrm{M}}^{\mathrm{U}}$ and ${T}_{\mathrm{M}}^{\mathrm{D}}$ being constants relatively insensitive to airfoils. Furthermore, the second condition of Eq. (35) is modified to avoid discontinuity which occurs at a large reduced frequency (e.g., k=0.2) for any LBbased models without recalibration of the time constant as
The effects of these modifications are displayed in Fig. 5.
The total firstorder dynamic response of the airfoil is formulated as
where the definition and description of each variable were given in Sect. 2.1 for the LB model. Thus the firstorder lift and drag responses can be obtained by
2.4.2 Secondorder correction
The secondorder correction takes the form of the nonlinear ordinary differential equation according to the secondorder correction of the Snel (Snel, 1997) and Adema–Snel (Adema et al., 2020) models. Generally, the basis of implementation of the present model mostly uses the Adema–Snel (Adema et al., 2020) model where the normal force is used instead of the lift force as for the original Snel model (Snel, 1997) as
This way, shedding effects apply not only to the lift force but also to the drag force. Note that τ is not present in Eq. (73) in contrast to the original formulation in Eqs. (44) and (55). The equation is changed because the time derivatives in the above equation are no longer with respect to time but to $s=\mathrm{2}Vt/c$, similar to the ONERA model (Tran and Petot, 1980; Dat and Tran, 1981; Petot, 1989). This is done to nondimensionalize the equations.
In Eq. (57), the k_{s} was projected as a function of the angle of attack by sin α_{n}. This modification causes problems when the hysteresis effect takes place in both positive and negative angles because Eq. (57) will be zero and then negative, causing instability of the ODE. Thus, the original form of the Snel model (Snel, 1997) is retained in the present model, but the constant is modified as
The idea for the downstroke damping as in Eq. (58) is adopted in the present model; the following form and constants are used:
Note again that τ is not present in the above equation. The original formulation in Eq. (58) replaces the damped oscillator when $\dot{{\mathit{\alpha}}_{n}}\le \mathrm{0}$ for a selfexcited oscillator of Van der Pol type with more damping. This is in contrast with the implementation done in Truong (2016), where the selfexcited oscillator is only replaced by the damped oscillator, when the flow is reattached on the return cycle. Under such circumstances, the oscillatory behavior still subsists in the return cycle, albeit with smaller amplitude. To accommodate this aspect, the last term of Eq. (75) is applied when the angle is smaller than ${\mathit{\alpha}}_{n}^{\mathrm{CRIT}}$. As for the forcing term, the original form of the Snel model (Snel, 1997) is adopted as
To facilitate the inclusion of the higher harmonic effects for the pitching moment, the idealized center of pressure obtained in the firstorder correction given in Eq. (26) is passed into the secondorder model. Thus, the dynamic moment reaction takes the form
Regarding the tangential force, a similar assumption is made as in Eq. (49) where the influence of $\mathrm{\Delta}{C}_{{\mathrm{T}}_{n}}^{\mathrm{D}\mathrm{2}}$ is neglected in the formulation, but $\mathrm{\Delta}{C}_{{\mathrm{T}}_{n}}^{\mathrm{D}\mathrm{1}}$ is considered in the firstorder term correction described above. Finally, the secondorder term of the lift ($\mathrm{\Delta}{C}_{{\mathrm{L}}_{n}}^{\mathrm{D}\mathrm{2}}$) and drag ($\mathrm{\Delta}{C}_{{\mathrm{D}}_{n}}^{\mathrm{D}\mathrm{2}}$) can be calculated easily. The effects of the secondorder term are displayed in Fig. 6.
2.5 Constants applied for the investigated dynamic stall models
The following constants are applied in the implemented dynamic stall models. These values are kept constant throughout the paper. The constants for the Leishman–Beddoes model and for the proposed IAG model are given in Tables 2 and 3, respectively. For any model developed based on the Leishman–Beddoes type, the critical angle of attack plays a major role. This can be obtained as the angle where the viscous pitching moment breaks or when the drag increases significantly (or the stall angle). The applied critical angles are given in Table 4. The validation is done by comparing the calculations with experimental data obtained at The Ohio State University for the S801 airfoil (13.5 % relative thickness) (Ramsay et al., 1996), NACA4415 airfoil (15 % relative thickness) (Hoffman et al., 1996), S809 airfoil (21 % relative thickness) (Ramsay et al., 1995) and S814 airfoil (24 % relative thickness) (Janiszewska et al., 1996). All selected test cases for the airfoils are employed with a leading edge grit (turbulator) to enable the “soiled” effects on a wind turbine blade at a Reynolds number of around 750 K. Note that these polar data are different from those used for example by Sheng et al. (2010) where the natural transition cases were taken. Therefore, the critical angles of attack are also different. The results of the studies are presented in the following sections.
The three stateoftheart dynamic stall models reviewed above (Leishman–Beddoes, Snel, Adema–Snel) have been used as a basis for examining the dynamic loads of four different pitching airfoils at various flow conditions. Experience gained from those models is used to formulate a new secondorder dynamic stall model, namely the IAG model, by evaluating the weakness and strength of each model. The presented secondorder models need to solve a set of differential equations. For this purpose, the Euler–Heun forward integration method is used.
3.1 Comparison against experimental data
This section compares the predicted dynamic forces and the measurement data. For a fair comparison, all models are assessed with the same time step size of $\mathrm{\Delta}t=T/\mathrm{1440}$, with T being the pitching period. The evaluations are performed on the S801 airfoil at k=0.073. The comparison of each model is shown in Figs. 7 to 9 for the Snel, Adema and IAG models, respectively. Note that the constants of the other existing dynamic stall models are taken directly from literature without further calibration for the S801 airfoil. Therefore, it is already expected that their performance will not be optimal. The main purpose of the comparison is not to study their accuracy but to analyze the robustness of each model for a different airfoil without tuning the constants. On the other hand, the constants for the IAG model in Table 3 were obtained using the S801 airfoil. To enable a fair assessment of the model robustness, the IAG model will also be used to reconstruct the dynamic polar data of four airfoils with different relative thicknesses without changing the constants in Sect. 3.6.
The original Snel models cannot predict the drag and moment coefficients in the original formulations. Thus, only the static polar data are shown. The Snel model actually shows an acceptable accuracy even though the constants are taken as found in literature. The higher harmonic effects are unfortunately not captured by this model. This is further refined by the Adema model which was developed as an improvement for the Snel model. The model performs fairly well for the lift and drag predictions, though the drag value at small angles of attack is a bit off. The pitching moment prediction is also not included in its formulations. These disadvantages are better treated in the proposed IAG model. Not only the prediction of the lift coefficient but also the accuracy of drag prediction are improved significantly. The modifications described in Sect. 2.4 result in the improvement at low and high angleofattack regimes. The model is also able to reconstruct the pitching moment polar accurately, which is important for aeroelastic calculations of wind turbine blades.
For the following sections, the proposed IAG model will be tested under various flow conditions and for several airfoils at various relative thicknesses in comparison with measurement data. Note that these calculations are performed without changing the constants to assess the robustness of the model at different flow conditions.
3.2 Effects of time signal deviation
The actual pitching motion within The Ohio State University (OSU) measurement differs slightly from the intended motion. The actual time series of the angle of attack is included in the experimental data (Ramsay et al., 1996; Hoffman et al., 1996; Ramsay et al., 1995; Janiszewska et al., 1996). To assess the effects of this time signal deviation on the aerodynamic response, the calculations using these time signal data were performed applying the IAG model. The results are compared with the experimental data and the results of the IAG model presented in Sect. 3.1. Note that these time signal data are fairly coarse and can cause problems for secondorder dynamic stall models because the gradient of α change can be extremely large. To cover for this issue, the time signal is interpolated in between each available point using a thirdorder cubicspline interpolation. The time signals are discretized by $\mathrm{\Delta}t=T/\mathrm{1440}$ over a single pitching period. The first period of oscillation is used for initialization of the time integration; thus a constant angle of attack is applied as shown in Fig. 10.
Figure 11 presents the influence of the time signal variation on the aerodynamic performance in terms of C_{L}, C_{D} and C_{M}. TS labels the exact time signals in the experimental campaign. Although the time signal difference has almost no influence on the global prediction characteristics, some deviation from the idealized sinusoidal motion can be noticed clearly. For example, one can see the increased lift buildup in the upstroke regime before stall and the location of the lift stall. Some deviations in the drag and pitching moment coefficients are observed as well. For the rest of the paper, the prediction procedure using the actual time signal from the experimental data is used for the best consistency with the experimental campaign.
3.3 Performance of the model for different mean angles of incidence
In this section, the effects of the mean angle of attack are evaluated. Three different angles of attack at the same inflow conditions are selected for this purpose. These are $\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}=\mathrm{8}$, 14 and 20^{∘}. Note that these mean angles of attack are only approximations since the actual time signal data from the experimental campaign are used. The selected mean angles represent the regime of attached flow, partly separated flow and fully separated flow conditions. These are helpful to assess the model performance under various flow situations.
Figure 12 presents the results for the lift coefficient under these three investigated mean angles of attack. The model performs very well for these different cases. The maximum lift is a bit overestimated in the model for the lowest $\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}$, but in general all unsteady lift characteristics in the measurement data are reproduced in a sound agreement with the experimental data. A similar behavior is shown for the drag prediction depicted in Fig. 13. The proposed model captures the increased drag effect and its shedding characteristics well. The simple modifications applied in Sect. 2.4 result in a good prediction of the drag coefficient compared with the experimental data. In Fig. 14, the prediction for pitching moment is shown. Here the predicted moment coefficient is in a good agreement with the measured values.
3.4 Performance of the model for different reduced frequencies
The effects of pitching frequency on the aerodynamic response will be discussed in this section. Three different reduced frequencies are examined, namely k=0.036, 0.073 and 0.111. The stall regime is shown here, where the prediction is the most challenging. The actual time signals as of the measurement campaign are used, following the procedure described in Sect. 3.2.
Figure 15 displays the results for the dynamic lift coefficient response. The lowest reduced frequency of 0.036 is dominated by the viscous effects. It represents the case where the “delayed” angleofattack response is the weakest. It can be seen that the maximum attained lift coefficient increases with increasing k. The gradient of the lift polar in the upstroke and downstroke phases also increases as well. These characteristics are present in both experimental data and predictions delivered by the IAG model. A similar behavior is also displayed in drag and pitching moments in Figs. 16 and 17, respectively. It is obvious that stall occurs much earlier for a smaller k value. One can see that the maximum amplitude of all three force components increases with increasing k. This can be dangerous for the structural stability, since the amplitude determines the fatigue loads.
To better investigate the effects of k, the IAG model is used to reconstruct the dynamic polar data at various k values by applying an idealized sinusoidal motion as presented in Fig. 18. Only the last DS cycle is shown for clarity of the observation. While the maximum amplitude of all three force components at lowfrequency domains increases with increasing k (blue and green markers), the amplitudes for all three forces at highfrequency domains show different characteristics as shown in the Fourier transformation in Fig. 19, albeit with much smaller values. The higher harmonic amplitudes are attributed to flow separation effects, while for lowfrequency domains are driven by the pitching motion (i.e., external unsteadiness or inflow).
3.5 Performance of the model for different pitching amplitudes
In this section, the effects of pitching amplitude on the aerodynamic response of a pitching airfoil are investigated. The mean angle of attack is fixed at $\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}=\mathrm{20}{}^{\circ}$. Note again that $\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}$ is only an approximation because the actual time signal data from the measurement campaign are applied. This large mean angle of attack is purposely selected because the poststall characteristic is of interest and is well known for its violent vibration, even for the static condition. The small amplitude in this case means that the whole pitch oscillation occurs within the stall regime.
Figures 20 to 22 display the dynamic force responses due to pitching motion of the airfoil predicted by the IAG model in comparison with the experimental data. The model accurately reconstructs the dynamic forces despite the predicted case being challenging within the poststall regime. Interesting to note is that the smallpitchingamplitude case induces stronger shedding effects for lift than the largeramplitude case. This can be explained as follows. As described by Leishman in his papers (Beddoes, 1982; Leishman, 1988; Leishman and Beddoes, 1989), the airfoil sees a lagged force response compared to the imposed disturbance. Therefore, in his model, a timelagged angle of attack is introduced as the “effective” angle actually seen by the airfoil section. When the pitching motion takes place partly within the fully separated region (in the static case) and partly in the attached/partly separated flow region, the airfoil still sees the lower angle (where the flow is still attached) even though the pitching motion already reaches the poststall regime. This effect stops/decreases when the effective angle is larger than the critical angle defined in Table 4. As the critical angle for the S801 airfoil is defined at 15.1^{∘}, the loweramplitude case is fully operating within the stall regime, where the attached flow effect is not present.
3.6 Performance of the model for different airfoils
In this section, the performance and robustness of the proposed IAG model are assessed for airfoils with different relative thickness. All model constants in Table 3 remain the same for all calculations. The difference from one airfoil calculation to the other lies only in the critical angle of attack value as shown in Table 4. The value was obtained simply by looking at the static polar data where the viscous pitching moment breaks or when the drag increases significantly.
Despite the increased airfoil thickness from 13.5 % to 24 %, Figs. 23 to 25 demonstrate that the reconstructed dynamic forces are in a good agreement with the experimental data, not only for the general trend but also the higher harmonic effects. As is also the case for the Leishman–Beddoes model, it is important to select the appropriate value for the critical angle of attack. The simple approach used in the present paper has shown its usefulness and potentially reduces the complexity of parameter tuning for industrial applications. Elgammi and Sant (2016) for example defined two different critical angles of attack, one for C_{N} and the other for C_{T} that were shown to improve the prediction accuracy. Although their attempt might be beneficial, this is not followed in this work because one main aim of the studies is to reduce parameter tuning required for one case to the other cases.
3.7 Predictions of the center of pressure
To further complement the analyses conducted in Sect. 3.6, the location of the actual pressure center is calculated in this section as
which indicates the distance of the pressure point to the quarter chord position where C_{M} is defined.
A correct location of the pressure point is important for determining the stability on aeroelastic simulations of wind turbine blades. The results of the calculations both for the experimental data and for the proposed IAG model are presented in Fig. 26 for all four investigated airfoils both as time series and as the polar plot. It can be seen clearly that the agreement between the experimental data and the present predictions is excellent for all investigated airfoils.
3.8 L_{2} norm of error analyses
Holierhoek et al. (2013) introduced a way of quantifying the absolute error between the experimental data and modeled lift coefficient. The general formulation reads
with ϕ being the variable of interest, i the current sample and N the total number of samples. In their paper, however, only lift was considered. Here all three force components will be shown for four different airfoils. Figure 27 displays the quantified error for two different flow categories, attached and deep stall. The time series of the angle of attack was obtained from the measured data by applying a thirdorder cubicspline interpolation in between each available point. One can see that generally the attached flow case is predicted very well, while the error increases as the flow condition becomes more complicated. Interestingly, especially for lift, it seems that the error decreases with increasing airfoil thickness. The reason for the larger error obtained for the thinner airfoil is attributed to the complex characteristics of the leading edge stall, causing severe load variations, especially with increasing angle of attack. Thus, it makes the prediction more challenging. Furthermore, the quantification of the error was also performed on two other dynamic stall models, Snel and Adema–Snel. The same approach for the angle of attack signal was applied. One can see that the IAG model shows its improved prediction in particular for the deepstall case for all three force components.
Comprehensive studies on the accuracy of several stateoftheart dynamic models to predict the aerodynamic loads of a pitching airfoil have been conducted. From the studies, the strength and weaknesses of each model were highlighted. This information was then transferred to develop a new secondorder dynamic stall model proposed in this paper. The new model improves the prediction for the aerodynamic forces and their higherharmonic effects due to vortex shedding, developed for robustness to improve its usability in practical wind turbine calculations. Details on the model characteristics, modifications and treatment for numerical implementation were summarized in the present paper. The studies were conducted by examining the influence of the time step size, time signal deviation, mean angle of attack, reduced frequency, pitching amplitude and variation in the airfoil thickness. Several main conclusions can be drawn from the work.

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

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

The dynamic forces reconstructed by the IAG model are in a sound agreement with the experimental data under various flow conditions by variation in $\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}$, k and Δα and for four different airfoils by changing only the values of the critical angle of attack.

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

When the airfoil operates at a high $\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}$ within the stall regime, a small Δα leads to increased vibrations for lift. The opposite is true for the pitching moment.
The present paper evaluates the newly developed IAG model under various flow conditions for four different airfoils. The following aspects are encouraged for future work.

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

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

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

Further tests and recalibration of the model for deepstall conditions at extremely large angles of attack are encouraged, which can be relevant for a turbine in standstill.
Variables  
s  Nondimensional time (s) 
V  Incoming wind speed (m s^{−1}) 
t  Time (s) 
c  Chord (m) 
C_{N}  Normal force coefficient (–) 
C_{T}  Tangential force coefficient (–) 
C_{L}  Lift force coefficient (–) 
C_{D}  Drag force coefficient (–) 
C_{M}  Pitching moment coefficient (–) 
${C}_{\mathrm{N}}^{\mathrm{P}}$  Total inviscid C_{N} (–) 
${C}_{\mathrm{N}}^{\mathrm{P}\mathrm{1}}$  Timelagged total inviscid C_{N} (–) 
${C}_{\mathrm{N}}^{\mathrm{I}}$  Impulsive inviscid C_{N} (–) 
${C}_{\mathrm{N}}^{\mathrm{C}}$  Circulatory inviscid C_{N} (–) 
${C}_{\mathrm{N}}^{\mathrm{f}}$  Viscous C_{N} (–) 
${C}_{\mathrm{T}}^{\mathrm{f}}$  Viscous tangential force coefficient (–) 
${C}_{\mathrm{M}}^{\mathrm{f}}$  Viscous pitching moment coefficient (–) 
${C}_{\mathrm{M}}^{\mathrm{C}}$  Circulatory pitching moment coefficient (–) 
${C}_{\mathrm{N}}^{\mathrm{V}}$  Vortex lift C_{N} (–) 
${C}_{\mathrm{N}}^{\mathrm{CRIT}}$  Critical C_{N} (–) 
C_{Pf}  Stepping parameter moment (–) 
f  Frequency (Hz) 
f_{0}  Pitching frequency (Hz) 
f_{n}, f_{1}, f_{2}  Separation factor (–) 
F_{1}, F_{2}  First and secondorder forcing term (–) 
k  Reduced frequency ($k=\mathit{\pi}{f}_{\mathrm{0}}c/V$) (–) 
K_{f}  Stiffness coefficient (–) 
k_{s}  Constant close to Strouhal number (–) 
M  Mach number (–) 
X, Y, D  Deficiency functions (–) 
c_{1}, c_{2}, c_{3}  Curvefitting constants (–) 
a_{1}, a_{2}, a_{3}  Curvefitting constants (–) 
S_{1}, S_{2}, S_{3}, α_{1}  Curvefitting constants (–) 
A_{1}, A_{2}, b_{1}, b_{2}  Model constants (–) 
K_{α}, T_{p}, T_{f}, T_{v}, T_{vl}  Model constants (–) 
K_{v}, η, ${K}_{\mathrm{f}}^{\mathrm{C}}$, ${T}_{\mathrm{m}}^{\mathrm{U}}$, ${T}_{\mathrm{m}}^{\mathrm{D}}$  Model constants (–) 
Greek letters  
α  Angle of attack (rad (unless stated otherwise)) 
α_{0}  Zero lift α (rad (unless stated otherwise)) 
α_{e}  Effective α (rad (unless stated otherwise)) 
α_{f}  Timelagged α_{e} (rad (unless stated otherwise)) 
α^{CRIT}  Critical α (rad (unless stated otherwise)) 
β  Machnumberdependent parameter (–) 
τ_{v}  Nondimensional vortex time (–) 
τ  Time constant (–) 
ζ_{v}  Vortex lift drag limiting factor (–) 
Subscripts  
_{n}  Present sampling time (–) 
_{f}  Viscous lagged value (–) 
_{v}  Vortex lift affected value (–) 
Superscripts  
^{INV}  Static inviscid (–) 
^{VISC}  Static viscous (–) 
^{I}  Impulsive (–) 
^{CRIT}  Critical (–) 
^{D}  Dynamic loading (–) 
^{D1}  Firstorder correction (–) 
^{D2}  Secondorder correction (–) 
The raw data of the simulation results can be shared by contacting the corresponding author of the paper.
GB developed the new model, designed the studies and conducted the analyses. TL and MA supported the research and provided suggestions and discussion about the manuscript.
The authors declare that they have no conflict of interest.
The authors gratefully acknowledge Wobben Research and Development GmbH for providing the research funding through the collaborative joint work DSWind. The measurement data provided from The Ohio State University are highly appreciated.
This research has been supported by the Wobben Research and
Development GmbH (grant no. DSWind).
This
openaccess publication was funded
by the University of
Stuttgart.
This paper was edited by Alessandro Bianchini and reviewed by Khiem V. Truong and Gerard Schepers.
Adema, N., Kloosterman, M., and Schepers, G.: Development of a secondorder dynamic stall model, Wind Energ. Sci., 5, 577–590, https://doi.org/10.5194/wes55772020, 2020. a, b, c, d, e, f, g, h, i
Bangga, G.: ThreeDimensional Flow in the Root Region of Wind Turbine Rotors, Kassel University Press GmbH, Kassel, https://doi.org/10.19211/KUP9783737605373, 2018. a
Bangga, G.: Numerical studies on dynamic stall characteristics of a wind turbine airfoil, J. Mech. Sci. Technol., 33, 1257–1262, https://doi.org/10.1007/s1220601902251, 2019. a
Beddoes, T.: Practical computation of unsteady lift, in: 8th European Rotorcraft Forum, AixenProvence, France, 1982. a, b, c
Carr, L. W., McAlister, K. W., and McCroskey, W. J.: Analysis of the development of dynamic stall based on oscillating airfoil experiments, Tech. rep., NASA TN D8382, National Aeronautics and Space Administration, Washington, D.C., USA, 1977. a
Dat, R. and Tran, C.: Investigation of the stall flutter of an airfoil with a semiempirical model of 2 D flow, TP no. 1981146, ONERA, p. 11, 1981. a, b
Elgammi, M. and Sant, T.: A Modified Beddoes–Leishman Model for Unsteady Aerodynamic Blade Load Computations on Wind Turbine Blades, J. Sol. Energ. Eng., 138, 051009, https://doi.org/10.1115/1.4034241, 2016. a, b, c
Faber, M.: A comparison of dynamic stall models and their effect on instabilities, MS thesis, Delft University of Technology, Delft, available at: https://repository.tudelft.nl/islandora/object/uuid:0001b1ebc19f48c3973d57eca4996a91 (last access: 5 January 2020), 2018. a, b, c
Galbraith, R.: Return from Airfoil Stall During RampDown Pitching Motions, J. Aircraft, 44, 1856–1864, 2007. a
Gupta, S. and Leishman, J. G.: Dynamic stall modelling of the S809 aerofoil and comparison with experiments, Wind Energy, 9, 521–547, 2006. a, b, c, d
Hansen, M. H., Gaunaa, M., and Madsen, H. A.: A Beddoes–Leishman type dynamic stall model in statespace and indicial formulations, Tech. rep., RisøR1354, Risø National Laboratory, Roskilde, Denmark, 2004. a, b, c
Hoffman, M., Ramsay, R., and Gregorek, G.: Effects of grit roughness and pitch oscillations on the NACA 4415 airfoil, Tech. rep., National Renewable Energy Lab., Golden, CO, USA, 1996. a, b, c, d, e
Holierhoek, J., De Vaal, J., Van Zuijlen, A., and Bijl, H.: Comparing different dynamic stall models, Wind Energy, 16, 139–158, 2013. a, b
Janiszewska, J., Ramsay, R., Hoffman, M., and Gregorek, G.: Effects of grit roughness and pitch oscillations on the S814 airfoil, Tech. rep., National Renewable Energy Lab., Golden, CO, USA, 1996. a, b, c, d, e
Khan, M. A.: Dynamic Stall Modeling for Wind Turbines, MS thesis, Delft University of Technology, Delft, available at: https://repository.tudelft.nl/islandora/object/uuid:f1ee9368ca4447caabe2b816f64a564f (last access: 5 January 2020), 2018. a
Larsen, J. W., Nielsen, S. R., and Krenk, S.: Dynamic stall model for wind turbine airfoils, J. Fluids Struct., 23, 959–982, 2007. a, b, c, d
Leishman, J.: Validation of approximate indicial aerodynamic functions for twodimensional subsonic flow, J. Aircraft, 25, 914–922, 1988. a, b, c
Leishman, J. G. and Beddoes, T.: A SemiEmpirical model for dynamic stall, J. Am. Helicopt. Soc., 34, 3–17, 1989. a, b, c, d, e, f, g, h, i, j, k
Martin, J., Empey, R., McCroskey, W., and Caradonna, F.: An experimental analysis of dynamic stall on an oscillating airfoil, J. Am. Helicopt. Soc., 19, 26–32, 1974. a
McAlister, K. W., Carr, L. W., and McCroskey, W. J.: Dynamic stall experiments on the NACA 0012 airfoil, Tech. rep., NASA Technical Paper 1100, National Aeronautics and Space Administration, California, USA, 1978. a
Øye, S.: Dynamic stall simulated as time lag of separation, in: Proceedings of the 4th IEA Symposium on the Aerodynamics of Wind Turbines, Rome, Italy, 1991. a
Pereira, R., Schepers, G., and Pavel, M.: Validation of the BeddoesLeishman Dynamic Stall Model for Horizontal Axis Wind Turbines using MEXICO data, in: 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 4–7 January 2011, Orlando, USA, AIAA 2011151, American Institute of Aeronautics and Astronautics (AIAA), Orlando, 2011. a
Petot, D.: Differential equation modeling of dynamic stall, La Recherche Aerospatiale, English Edn., ONERA, 59–72, 1989. a, b
Ramsay, R., Hoffman, M., and Gregorek, G.: Effects of grit roughness and pitch oscillations on the S809 airfoil, Tech. rep., National Renewable Energy Lab., Golden, CO, USA, 1995. a, b, c, d, e
Ramsay, R., Hoffman, M., and Gregorek, G.: Effects of grit roughness and pitch oscillations on the S801 airfoil, Tech. rep., National Renewable Energy Lab., Golden, CO, USA, 1996. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v
Sheng, W., Galbraith, R. M., and Coton, F.: A new stallonset criterion for low speed dynamicstall, J. Sol. Energ. Eng., 128, 461–471, 2006. a
Sheng, W., Galbraith, R., and Coton, F.: A modified dynamic stall model for low Mach numbers, J. Sol. Energ. Eng., 130, 031013, https://doi.org/10.1115/1.2931509, 2008. a
Sheng, W., Galbraith, R. A. M., and Coton, F. N.: Applications of lowspeed dynamicstall model to the NREL airfoils, J. Sol. Energ. Eng., 132, 011006, https://doi.org/10.1115/1.4000329, 2010. a
Snel, H.: Heuristic modelling of dynamic stall characteristic, in: Proceedings of the European Wind Energy Conference, Dublin, Ireland, 1997. a, b, c, d, e, f, g, h, i
Tarzanin, F.: Prediction of control loads due to blade stall, J. Am. Helicopt. Soc., 17, 33–46, 1972. a
Tran, C. and Petot, D.: Semiempirical model for the dynamic stall of airfoils in view of the application to the calculation of responses of a helicopter blade in forward flight, in: 6th European Rotorcraft Forum, Bristol, UK, 1980. a, b, c
Truong, K.V.: Modeling aerodynamics for comprehensive analysis of helicopter rotors, in: Proceedings of 42nd European Rotorcraft Forum, Lille, France, 2016. a
Truong, V.: A 2d dynamic stall model based on a hopf bifurcation, in: Proceedings of 19th European Rotorcraft Forum, Cernobbio, Italy, 1993. a, b
Wang, Q. and Zhao, Q.: Modification of Leishman–Beddoes model incorporating with a new trailingedge vortex model, Proc. Inst. Mech. Eng. Pt G, 229, 1606–1615, 2015. a