1. Modelling dynamic stall
The phenomenon of dynamic stall refers to the massive flow separation experienced by an airfoil due to an unsteady change in angle of attack beyond the critical static stall angle (McAlister, Carr & Mccroskey Reference McAlister, Carr and Mccroskey1978; McCroskey Reference McCroskey1981; Carr Reference Carr1988). Dynamic stall is associated with a delay in stall onset and stall recovery with respect to static flow separation and boundary-layer reattachment. The delay in stall onset is related to an increase in the maximum attainable lift but when separation occurs, it leads to a rapid aft movement of the centre of pressure, large nose-down pitching moments on the airfoil section, and an increase in torsional loads on the wings (McAlister, Lambert & Petot Reference McAlister, Lambert and Petot1984). The unsteady dynamic stall load response reduces the aerodynamic efficiency, limits operational and performance boundaries, and decreases structural stability of lifting surfaces operating at high angles of attack or under highly unsteady flow conditions. Accurate prediction of the unsteady load response is key to design efficient measures to eliminate or control massive flow separation and improve the design and operation of wind turbines, unmanned aerial vehicles, rapidly manoeuvring aircraft and helicopters (Leishman Reference Leishman2002; Reich et al. Reference Reich, Eastep, Altman and Albertani2011; Uhlig & Selig Reference Uhlig and Selig2017; Smith et al. Reference Smith, Jones, Ayancik, Mulleners and Naughton2020).
The characteristic features identified within a full dynamic stall life-cycle include attached flow, the emergence and spreading of flow reversal on the airfoil's suction side, the formation of a large-scale coherent structure referred to as the dynamic stall vortex, the separation of the first dynamic stall vortex indicating stall onset, massive flow separation and potentially flow reattachment (Carr, McAlister & McCroskey Reference Carr, McAlister and McCroskey1977; Shih et al. Reference Shih, Lourenco, Dommelen and Krothapalli1992; Mulleners & Raffel Reference Mulleners and Raffel2013). Each of these phenomena is characterised by different time scales and governed by different parameters (Deparday & Mulleners Reference Deparday and Mulleners2019; Ansell & Mulleners Reference Ansell and Mulleners2020). The succession of different flow stages and the interplay between the associated characteristic phenomena and time scales make the modelling and accurate prediction of the dynamic stall load response challenging.
One of the most prominent dynamic stall models in literature, the Beddoes–Leishman model, manages to represent the complexity of the dynamic stall through four interconnected submodels that predict the aerodynamic coefficients during different phases of dynamic stall (Leishman & Beddoes Reference Leishman and Beddoes1989). An attached flow submodel represents the unsteady but linear aerodynamic forces prior to stall onset as a sum of non-circulatory and circulatory loading. A separated flow submodel models the separation of the boundary layer, considering the formation and evolution of the leading-edge vortex. It represents the leading-edge and the trailing-edge separation as a function of losses in the circulation and non-linear forces. A dynamic stall onset model predicts the timing of stall initiation and a dynamic stall vortex model considers the vortex induced aerodynamic forces. The Beddoes-Leishman model is widely used for helicopter and wind turbine applications and it accurately predicts the dynamic stall load response while offering a physics-based representation of the stall process for incompressible flows.
Various attempts have been made in the past decades to improve or modify individual aspects of the Beddoes-Leishman model. In a series of papers, Sheng, Galbraith & Coton (Reference Sheng, Galbraith and Coton2006) introduced and improved a new stall-onset criteria for low-Mach-number flows (${{Ma}} < 0.3$) (Sheng et al. Reference Sheng, Galbraith and Coton2006; Sheng, Galbraith & Coton Reference Sheng, Galbraith and Coton2008). Hansen, Gaunaa & Madsen (Reference Hansen, Gaunaa and Madsen2004) neglected the compressibility effects and the leading-edge flow separation of the Beddoes-Leishman model to develop a model that is better suited to wind turbine applications. Larsen, Nielsen & Krenk (Reference Larsen, Nielsen and Krenk2007) proposed a computationally cheaper version of the Beddoes-Leishman model by expressing the static-lift curve as a function of the fully attached flow-lift value and the degree of attachment. The Beddoes-Leishman model and derivatives provide a detailed representation of the flow physics and predict the dynamic stall load response but they typically rely on more than 10 empirical parameters that need to be extracted and tuned based on experiments or numerical simulations.
For active or closed-loop flow control applications, the level of detail provided by Beddoes-Leishman-type stall models might not be necessary and state-space models with fewer empirical parameters and reduced complexity are more desirable. Data-driven approaches combined with system identification techniques are also better suited for modern control application as they allow for a broad range of input conditions, and can be designed to deal well with aggressive pitch-angle manoeuvres (Brunton, Dawson & Rowley Reference Brunton, Dawson and Rowley2014; Hemati, Dawson & Rowley Reference Hemati, Dawson and Rowley2017; Le Provost, Williams & Brunton Reference Le Provost, Williams and Brunton2018). However, most purely data-driven models do not provide much physical insight into the dynamic stall process. The state-space Goman–Khrabrov model (Goman & Khrabrov Reference Goman and Khrabrov1994) is a first-order model with two empirical parameters that provides a good compromise between interpretability and simplicity. The model has received renewed interest in the past decade for active control applications (Williams et al. Reference Williams, An, Iliev, King and Reißner2015, Reference Williams, Greenblatt, Müller-Vahl, Santra and Reißner2019; Sedky, Jones & Lagor Reference Sedky, Jones and Lagor2020; An et al. Reference An, Williams, Eldredge and Colonius2021). It accurately models dynamic stall load hysteresis over a wide range of reduced frequencies and flow conditions where dynamic stall is characterised by a progressive increase of trailing-edge separation. This is the case for most occurrences of dynamic stall on airfoils at moderate and high Reynolds numbers. The Goman–Khrabrov model is not applicable to flat-plate leading-edge stall cases.
2. Goman–Khrabrov dynamic stall model
The Goman–Khrabrov model is a nonlinear state-space model that represents the aerodynamic forces and moments produced by unsteady flows with trailing-edge separation and vortex breakdown on wings at high angles of attack (Goman & Khrabrov Reference Goman and Khrabrov1994). The model predicts the degree of flow attachment on the upper surface of the airfoil or the position of the vortex breakdown described by a single internal state variable $X$. The value $X = 1$ represents completely attached flow, and $X = 0$ corresponds to fully separated flow. The state equation for $X$ is given by
Here, $ {{X}_{{0}}} (\alpha )$ is the location of the separation point as a function of the angle of attack under static stall conditions, $ {{\tau }_{{1}}}$ is the relaxation time constant related to the relaxation of an unsteady force response to a steady state, $ {{\tau }_{{2}}}$ corresponds to the stall delay time constant, and $\dot {\alpha }$ indicates the instantaneous pitch rate of the airfoil. The left-hand side of (2.1) is a first-order linear differential model that describes the evolution of the degree of flow attachment during unsteady manoeuvres of the airfoil. The right-hand side is a nonlinear function that returns the value of $X$ for the different angles of attack based on the static response. Once $X$ is found, the lift response is obtained using Kirchhoff's law:
where ${{{\rm d}} {{C}_{{l}}}}/{{{\rm d}} \alpha }|_{0}$ is the static lift slope.
The Goman–Khrabrov model includes two empirical parameters $ {{\tau }_{{1}}}$ and $ {{\tau }_{{2}}}$ and requires the static evolution of the separation point location with angle of attack $ {{X}_{{0}}}(\alpha )$ as input. The separation point location $ {{X}_{{0}}}(\alpha )$ is typically determined by fitting the Kirchhoff law (2.2) to the static lift polar $ {{C}_{{l}}}(\alpha )$ of an airfoil at a given Reynolds number. Once $ {{X}_{{0}}}(\alpha )$ is obtained, $ {{\tau }_{{1}}}$ and $ {{\tau }_{{2}}}$ are empirically determined by applying (2.1) to unsteady experimental or numerical training data. The best-fit values of $ {{\tau }_{{1}}}$ and $ {{\tau }_{{2}}}$ are considered to be constant for a large range of pitching frequencies for a specific airfoil shape (Williams et al. Reference Williams, An, Iliev, King and Reißner2015; Le Provost et al. Reference Le Provost, Williams and Brunton2018; Williams & King Reference Williams and King2018; An et al. Reference An, Williams, Eldredge and Colonius2021).
We have applied this procedure to dynamic stall data of a sinusoidally pitching OA209 airfoil at a chord-based Reynolds number $Re = 9.2 \times 10^{5}$ and Mach number ${{Ma}}={0.14}$. These data have been previously presented in various publications (Mulleners & Raffel Reference Mulleners and Raffel2012, Reference Mulleners and Raffel2013; Deparday & Mulleners Reference Deparday and Mulleners2019; Ansell & Mulleners Reference Ansell and Mulleners2020). The data has been obtained in an open test section of a closed-loop wind tunnel. For this study, the data are corrected with the classical open test section corrections specified in Ewald (Reference Ewald1998). Pressure-based lift coefficients per unit span are obtained by numerically integrating the instantaneously measured chord-wise pressure distributions in the airfoil's mid-span cross-section. The lift coefficient per unit span is defined as: $ {{C}_{{l}}} = L/(0.5\rho {{{U}}_{{\infty }}}^{2} c)$, where $L$ is the lift per unit span, $\rho$ is the air density, $ {{{U}}_{{\infty }}}$ is the wind tunnel velocity and $c$ is the airfoil chord length. The static lift polar was obtained by quasi-steadily increasing the angle of attack starting from well below the critical stall angle.
An exemplary result of the Goman–Khrabrov lift coefficient prediction for a sinusoidal pitching motion around a mean angle of attack $ {{\alpha }_{{0}}}=20^{\circ }$, with an amplitude $ {{\alpha }_{{1}}}=8^{\circ }$ and a reduced frequency $k={0.05}$ is presented in figure 1. The lift coefficient is presented as a function of convective time, with the origin corresponding to the time, $ {{t}_{{ss}}}$, at which the static stall angle, $ {{\alpha }_{{ss}}}$, is exceeded. Selected angle-of-attack values are indicated in the axis on top. The Goman–Khrabrov prediction is the best fit of the phase-averaged lift evolution that was obtained treating $ {{\tau }_{{1}}}$ and $ {{\tau }_{{2}}}$ as fitting parameters in (2.1). The model result is compared with the experimentally measured lift evolution in figure 1. The black line corresponds to the experimental data for a single pitching cycle, the coloured line to the predicted lift coefficient, and the shaded area to the envelope covered by all measured pitching cycles. The Goman–Khrabrov model predicts well the attached flow regime, the magnitude and timing of the first lift peak, and the general decaying lift trend post dynamic stall, which occurs at $t= {{t}_{{ds}}}$. It is a first-order model that neglects higher-order features in the flow field and cannot reproduce the higher harmonics of the load fluctuations due to subsequent vortex shedding seen in figure 1 (Culler & Farnsworth Reference Culler and Farnsworth2019; An et al. Reference An, Williams, Eldredge and Colonius2021).
The values of $ {{\tau }_{{1}}}$ and $ {{\tau }_{{2}}}$ that yield the best-fit predictions of the lift evolution are presented in figure 1(b,c) for various pitching kinematics. The different pitching kinematics are characterised by their normalised effective unsteadiness $ {{\dot {\alpha }}_{{ss}}} c / (2 {{{U}}_{{\infty }}})$. Here, $ {{\dot {\alpha }}_{{ss}}}$ is the rate of change of the angle of attack when the static stall angle is exceeded. This instantaneous pitch rate serves as representative effective pitch rate to characterise the time scales of dynamic stall development for sinusoidal motions (Mulleners & Raffel Reference Mulleners and Raffel2012; Deparday & Mulleners Reference Deparday and Mulleners2019; Kissing et al. Reference Kissing, Kriegseis, Li, Feng, Hussong and Tropea2020; Le Fouest, Deparday & Mulleners Reference Le Fouest, Deparday and Mulleners2021). The values for $ {{\tau }_{{1}}}$ that yield the best-fit results do not show any dependence on the motion unsteadiness and reach an average value of 0.026(9) s for all pitching motions considered here. The values for $ {{\tau }_{{2}}}$ show a general decrease with increasing effective unsteadiness of the driving motion. In most prominent previous work that uses the Goman–Khrabrov model, the time coefficients have been manually tuned and constant values have been selected for $ {{\tau }_{{1}}}$ and $ {{\tau }_{{2}}}$ for a given airfoil profile for all pitching frequencies (Williams et al. Reference Williams, An, Iliev, King and Reißner2015). The empirical selection of the time constants is the main drawback of the Goman–Khrabrov model, which limits it range of applicability.
Following the construction of the Goman–Khrabrov model, the time constants $ {{\tau }_{{1}}}$ and $ {{\tau }_{{2}}}$ are physically interpreted as a relaxation time and a stall delay, respectively. Yet, so far, no solutions have been proposed to replace the empirical parameters in the Goman–Khrabrov dynamic stall model by physically derived time scales. Here, we aim to close that gap by introducing functional relationships to determine the parameters $ {{\tau }_{{1}}}$ and $ {{\tau }_{{2}}}$ based solely on the kinematic input parameters.
3. Novel physics-based time constants
We start with the second time constant $ {{\tau }_{{2}}}$, which is the stall delay constant. Instead of empirically fixing this constant for a range of reduced frequencies, we will directly connected $ {{\tau }_{{2}}}$ to the physical dynamic stall delay between the time at which the static stall angle is exceeded and the onset of dynamic stall. The stall delay constant represents the delay in the dynamic stall angle of attack due to an unsteady pitching motion to:
The difference between the static and the dynamic stall angle can also be found by integrating the pitching motion during the dynamic stall delay, such that
with $ {{{\rm \Delta} t}_{{ds}}}$ the temporal dynamic stall delay. For constant pitch rate ramp-up motions, (3.2) yields
and the Goman–Khrabrov stall delay constant equals the temporal dynamic stall delay (Ericsson & Reding Reference Ericsson and Reding1988). For sinusoidal pitching motions with effective unsteadiness ${{\dot{\alpha}}_{{ss}}}$, (3.2) yields
For sinusoidal pitching motions, the Goman–Khrabrov stall delay constant is again a function of the temporal dynamic stall delay and three kinematic input parameters: the pitching amplitude, the pitching frequency and the effective unsteadiness.
Equations (3.3) and (3.4) would only truly aid to reduce the empiricism of the Goman–Khrabrov model if the dynamic stall delay is known. Recently, Le Fouest et al. (Reference Le Fouest, Deparday and Mulleners2021) compared the stall delays for three different airfoils, a NACA0015, a NACA0018 and an OA209, at different Reynolds numbers ranging from $7.5 \times 10^{4}$ to $10^{6}$, and covering linear ramp-up and sinusoidal pitching motions. The data points in figure 2 have different colours and marker shapes referring to the different airfoil shapes, Reynolds number values, and motion kinematics according to the legend on the right of the figure. All measured dynamic stall delays decrease with increasing motion unsteadiness and converge towards a minimum value for high values of the unsteadiness. This overall decay of the stall delay is fitted by a power law:
for all deep dynamic stall cases considered, regardless of the airfoil geometries, motion kinematics and Reynolds number (figure 2). The $R^2$-value of the fit is 0.978.
The overall stall delay is the sum of the delays attributed to the two stages that are identified in the dynamic stall development process (Mulleners & Raffel Reference Mulleners and Raffel2013; Deparday & Mulleners Reference Deparday and Mulleners2019). Dynamic stall development includes a reaction stage, which is characterised by accumulation of vorticity in the shear layer, followed by a vortex formation stage, during which the shear layer rolls up into a large-scale dynamic stall vortex. These two stages are characterised by different time responses. The duration of the first stage decreases with increasing motion unsteadiness. The duration of the second stage is independent of the motion unsteadiness and indicates the minimum time required for a dynamic stall vortex to form. The lower limit of the stall response $ {{{\rm \Delta} t}_{{ds}}} ( {{\dot {\alpha }}_{{ss}}}\rightarrow \infty )$ corresponds to the duration of the vortex formation stage and is of the order of the general characteristic vortex formation numbers observed for a variety of vortex generators (Dabiri Reference Dabiri2009).
Using the universal expression for the dynamic stall delay provided by (3.5), we can find $ {{\tau }_{{2}}}$ solely based on the input parameters. This reduces the empiricism in the Goman–Khrabrov dynamic stall model and generalises it.
The time constant $ {{\tau }_{{1}}}$ is the relaxation time constant which represents the relaxation time for a perturbed flow to recover to a steady state (Goman & Khrabrov Reference Goman and Khrabrov1994). This steady state can be an attached flow state or a fully developed separated flow state with bluff body vortex shedding. The post-stall regime of a deep dynamic stall cycle is characterised by repeated large-scale vortex shedding which leads to large fluctuations in the lift. The post-stall lift fluctuations due to vortex shedding are most clearly visible in the lift evolution of a single pitching cycle represented by the black curve in figure 1. Increased flow perturbations post stall cause fluctuations or jitter in the post-stall vortex formation and shedding between successive pitching cycles. This leads to the so-called cycle-to-cycle variations in the lift evolution. The shaded area in figure 1 indicates the region between the minimum and maximum values observed across a total of 40 pitching cycles. When the flow is attached, at the beginning of the cycle up to dynamic stall onset at $ {{t}_{{ds}}}$ and after reattachment, there are little to no cycle-to-cycle variations. The cycle-to-cycle variations appear post stall and can lead to broadening or even elimination of the post-stall lift fluctuations in the phase-averaged lift cycle.
The Goman–Khrabrov is a first-order model and cannot predict the post-stall lift fluctuations. To provide insight into the characteristic time scales that govern the post-stall lift decay, we present the magnitude of the post-stall lift peaks and their timing for multiple pitching cycles and different sinusoidal pitching motions of the OA209 airfoil at $Re = 9.2\times 10^{5}$ in figure 3. The timing of the lift peaks is presented with respect to the onset of dynamic stall. The colour of the markers indicates the normalised effective unsteadiness that characterises the individual sinusoidal pitching kinematics ($ {{\alpha }_{{0}}}\in \{18^{\circ },20^{\circ },22^{\circ }\}$, $ {{\alpha }_{{1}}}\in \{6^{\circ },8^{\circ }\}$, $k\in \{0.05,0.075,0.1\}$).
The lift peaks corresponding to different sinusoidal pitching motions align surprising well despite the large variation in the angle of attack at which they occur. Independent of the motion unsteadiness and the local angle of attack, the local lift peaks decay with increasing convective time post stall. The decay has two parts. The initial decay is steep and covers approximately 12 convective times. The lift peaks during this initial phase seem to cluster in three groups with a spacing of approximately four to five convective times. This corresponds to the lower limit of the stall response $ {{{\rm \Delta} t}_{{ds}}} ( {{\dot {\alpha }}_{{ss}}}\rightarrow \infty )=4.24 c/ {{{U}}_{{\infty }}}$ in figure 2. An interval of 4.24 convective times is equivalent to a Strouhal number (St) of 0.235. The decay of the lift peaks becomes less steep but remains more or less linear for longer convective times post stall. As the post-stall vortex shedding is independent of the motion unsteadiness and governed by convective time scales, we propose to replace the relaxation time constant $ {{\tau }_{{1}}}$ by the lower limit of the stall response or
For the data presented in figure 4, this would yield $ {{\tau }_{{1}}} = 0.0254\ {\rm s}$ which is in the middle of the best-fit values.
To validate the performance of the Goman–Khrabrov model with our physics-based time constants, we have applied the model to the three experimental data sets mentioned above. The resulting predictions of the Goman–Khrabrov model are compared with a selected ensemble-averaged response for the ramp-up motion and phase-averaged lift responses for the pitching motions in figure 4(a–c). The model based on our proposed time scales predicts well the overall evolution of the lift, the timing of the first lift peak, and the general post-stall lift decay. The magnitude of the first lift peak is well predicted for the sinusoidal motions, but under-predicted for the ramp-up motion. Yet, this would be the same for the best-fit Goman–Khrabrov prediction.
To qualitatively assess the performance of the newly introduced time constants, we calculated two metrics of interest for the best-fit model and our modified version for all dynamic stall cases within our three experimental data sets. The first metric is the coefficient of determination, or $R^{2}$-value. This metric gives us a general measure of how well the experimental data are replicated by the model. The second metric is the relative time difference between the timing of the main lift peak in the experimental data and the model prediction $(\Delta t{\rm Cl},\max)$. This is a more local measure that is of particular importance for the prediction of dynamic stall. The results of both metrics are summarised as a function of the effective unsteadiness, which characterises the motion, in figure 4(d,e). Purple markers indicate results from the best-fit Goman–Khrabrov model and orange markers indicate the results of the model with the novel physics-based time constants. The shapes of the markers refer to the different data sets as indicated in figures 2 and 4(a–c).
In general, with $R^{2}$-values above 0.85, both Goman–Khrabrov predictions provide an accurate model of all the experimental data considered here. The predictions using the time constants based on (3.6) and (3.2) have $R^{2}$-values that are only marginally lower than the best-fit results but they provide a better prediction of the timing of the first lift peak.
4. Conclusion
We have generalised the Goman–Khrabrov dynamic stall model and reduced its empiricism by introducing physically derived time scales to replace the empirical parameters in the original model. The first time constant represents the shedding period of the post-stall load fluctuations. The post-stall decay rate is independent of the motion kinematics and the time constant is obtained assuming Strouhal number ${{St}}=0.25$. The second time constant represents the temporal dynamic stall delay. We presented a generally applicable functional relationship between the non-dimensional stall delay and the normalised instantaneous pitch rate at static stall. This relationship is largely independent of the type of the motion, the Reynolds number, and the airfoil geometry based on three available dynamic stall data sets. Using this general expression for the dynamic stall delay, we compute the second time constant directly from the input parameters. The Goman–Khrabrov model with our newly defined time constants gives excellent predictions of the unsteady lift evolutions for a large set of experimental dynamic stall data covering different airfoil shapes, motion types, a range of Reynolds numbers from $7.5 \times 10^{4}$ and $10^{6}$ and a large range of reduced frequencies. The use of our physics-based time constants generalises the Goman–Khrabrov dynamic stall model and opens new opportunities for closed-loop flow control applications.
Funding
The work presented is supported by the SNSF Assistant Professor energy grant number PYAPP2_173652.
Declaration of interests
The authors report no conflict of interest.