A review of propeller stall flutter

Abstract

Research on propeller performance has been reinvigorated by the development of new classes of vehicles, ranging from electrically powered fixed-wing aircraft, to multi-rotor electrical Vertical Take-off/Landing (eVTOL) vehicles and tilt-rotor aircraft. These types of aircraft utilise a range of modern propellers, often with more advanced planforms and features such as anhedral, and operate in flight envelopes that are outwith the traditional bands of performance. The use of advanced materials (mostly composites), high geometrical sweeps and variable angular velocities are the source of unsteady aerodynamics, that is often coupled with the blade’s structural response. Data from experimental investigations is mostly historic, with the majority of studies conducted before 1960, when aviation shifted rapidly towards jet propulsion. These studies lack in flutter boundary assessment. Modern propellers are likely to be pushed toward their flutter boundaries, but the experimental database published to-date is insufficient to provide flutter boundary assessment. This review examines the value of the available experimental research and the status of the state-of-the-art numerical methods, in order to establish the requirements for modern research on propeller stall flutter.

Nomenclature

Latin Symbols

Ar

non-dimensional rate of angle-of-attack

a

axis of rotation for a pitch and plunge aerofoil

$a_{on}$

ONERA empirical parameter

b

aerofoil semi-chord

[C]

damping matrix

C(k)

theodorsen function

$C_{Am}$

Gangwani moment stall empirical parameter

$C_{At}$

Gangwani trailing egde empirical parameter

$C_{L_0}$

zero-angle lift coefficient

$C_{L1}$

ONERA linear aerodynamics lift coefficient

$C_{L2}$

ONERA non-linear aerodynamics lift coefficient

$C_{L_c}$

circulatory lift coefficient

$C_{L\ell}$

static lift coefficient

$C_{L_{nc}}$

non-circulatory lift coefficient

$C_M$

pitching moment coefficient

$C_{M_0}$

zero-angle pitching moment coefficient

$C_{M1}$

ONERA linear aerodynamics pitching moment coefficient

$C_{M2}$

ONERA non-linear aerodynamics pitching moment coefficient

$C_{M_c}$

circulatory pitching moment coefficient

$C_{M\ell}$

static pitching moment coefficient

$C_{M_{nc}}$

non-circulatory pitching moment coefficient

$C_N$

normal force coefficient

$C_{wR}$

Gangwani reattachment empirical parameter

$C_{Wm}$

Gangwani moment stall empirical parameter

c

chord

$c_{\ell}$

lift coefficient

$c_{l_{\alpha}}$

lift curve slope

$c_{m,0.25c}$

pitching moment coefficient at quarter chord point

E

Young’s modulus

EF

stress excitation factor

$E_{on}$

ONERA empirical parameter

e

distance between mass and elastic axes

GJ

torsional stiffness

g

structural damping coefficient

$H_n^{(2)}$

Hankel function of the second kind

h

aerofoil plunging motion

$I_1$

chordwise moment of inertia

$I_y$

mass moment of inerta

[K]

stiffness matrix

$K_1$

Gormont pitch rate coefficient

k

reduced frequency

$k_0$

non-dimensional aerodynamic centre offset from the quarter chord

$k_{1/2}$

Leishman-Beddoes semi-empirical coefficients

$k_y$

blade area radii of gyration

$k_{xx},k_{zz}$

blade mass radii of gyration

L

aerodynamic lift

LF

Kirchoff-Helmholtz unsteady lift factor

M

aerodynamic pitching moment

${\rm M}_a$

Mach number

m

section mass

$m_i$

generalised mass for i-th mode

n

number of modes

PF

Kirchoff-Helmholtz unsteady pitching moment factor

Q

number of blades

$Q_i$

generalised force for i-th mode

$q_i$

modal coordinate for i-th mode

R

blade radius

$R_{hub}$

hub radius

Re

Reynolds number

r

non-dimensional blade radius

$r_{on}$

ONERA empirical parameter

s

aerodynamic time variable, Wagner function

$s_{1/2}$

Leishman-Beddoes semi-empirical coefficients

$s_{{1/2}_n}$

negative angle Leishman-Beddoes semi-empirical coefficients

$s_{mt}$

time at which the vortex detaches from the trailing edge

$s_{on}$

ONERA empirical parameter

Ta

blade axial tension

t

time

$t/c$

thickness/Chord ratio

U

freestream velocity

$V_{eq}$

equivalent airspeed

$x_a$

distance from axis of rotation to aerodynamic centre

$x_s$

non-dimensional static upper surface separation point

$x_{s_E}$

Leishman-Beddoes effective separation point

$y_{CG}$

centre of gravity offset

$y_{EA}$

elastic axis offset

z

blade flapping motion

Greek symbols

$\alpha$

angle-of-attack

$\alpha_0$

initial angle-of-attack

$\alpha_{Dm}$

pitching moment stall angle

$\alpha_E$

effective angle-of-attack

$\alpha_{QS}$

quasi-steady angle-of-attack

$\alpha_s$

static-stall angle-of-attack

$\alpha_{s_n}$

negative static-stall angle-of-attack

$\alpha_W$

unsteady decay parameter

$\alpha_{Wm}$

unsteady decay parameter at moment stall

$\beta$

Prandtl-Glauert compressibility correction

$\gamma$

blade taper ratio

$\gamma_g$

Gormont non-dimensional empirical parameter

$\Delta \alpha_{DS}$

incremental dynamic stall angle

$\Delta C_L$

difference in lift coefficient

$\Delta C_M$

difference in pitching moment coefficient

$\epsilon$

Gangwani empirical parameter

$\zeta_i$

critical damping for the i-th mode

$\theta$

aerofoil pitching motion

$\Lambda$

blade sweep

$\lambda_{1/2}$

Dowell’s Theodorsen function coefficients

$\lambda_{on}$

ONERA empirical parameter

$\rho$

air density

$\sigma$

rotor solidity

$\sigma_{on}$

ONERA empirical parameter

$\tau$

non-dimensional time

$\tau_1$

semi-empirical relaxation time constant

$\tau_2$

Leishman-Beddoes time-delay parameter

$\phi_c(s,M_a)$

compressibility corrected Wagner function

$\psi$

phasing between blades

$\psi_i$

freestream inflow angle

$\Omega$

rotational velocity

$\omega$

natural frequency

$\omega_i$

natural frequency for the i-th mode

Acronyms

BEMT

Blade Element Momentum Theory

CFD

Computational Fluid Dynamics

CG

Centre of gravity

eVTOL

electrical Vertical Take-Off/Landing

HMB3

Helicopter Multi-Block 3

NACA

National Advisory Committee for Aeronautics

NASA

National Aeronautics and Space Administration

SAS

Scale-Adaptive Simulation

SST

Shear Stress Transport

URANS

Unsteady Reynolds-Averaged Navier-Stokes

UTRC

United Technologies Research Centre

1.0 Introduction

There are at least three known types of propeller flutter: bending-torsion flutter, stall flutter and whirl flutter. Bending-torsion flutter, or classical flutter, is driven by the coupling and excitement of selected modes of vibration. For this type of flutter, the aerodynamics can often operate in the linear regime with non-linearity introduced via the coupling of modes due to the combined effects of structural and aerodynamic forces. To understand the coupling procedure, a step-by-step approach can be taken. Any change in effective angle-of-attack, whether this be driven by internal or external effects such as inflow fluctuations, results in changes to the blade lift. This results in a flapping of the blade which triggers the torsional mode. Depending on the current freestream velocity, this excitation can grow, dampen or self-sustain. As a result of the linear aerodynamics, lower fidelity, fast models can be utilised which predict the critical flutter boundaries accurately. Although a good level of accuracy is found when assessing for the critical flutter points, the models ability to predict the blade response outside of this critical flutter point loses accuracy. Moreover, such tools can be used to understand and, potentially, adjust such points to ensure that the propeller operation avoids such excitement.

Although several sources of non-linearity can be found within a propeller blade when looking at classical flutter, such as the the four sources of deflection arising from non-linear inertial loading, a type of flutter which is inherently more non-linear than classical flutter is that of whirl flutter. This is due to the combination of blade-wing-nacelle structure with the propeller gyroscopic and aerodynamic effects adding to the sources of non-linearity already found within a standard propeller. Such a response can be found in helicopters, turbo-prop aircraft and tiltrotors, with turbo-prop aircraft being the principle test case. Whirl flutter for tiltrotors is more complex than for standard turbo-props due to the ability of the tiltrotor to operate in different flight regimes, including helicopter, transition and propeller modes. As a result of the aerodynamics associated with propeller systems, non-linear aerodynamics becomes more critical and, hence, the modelling becomes more complex. Therefore, the combination of more complex structural models with non-linear aerodynamics significantly increases the sophistication of the modelling requirement. Simpler models that assume rigid propellers, for example, are available but these often produce conservative flutter boundaries with respect to the experiment. Although conservative boundaries, from an operational point-of-view and ensuring efficiency of the flight path, is not terminal as whirl flutter is typically a catastrophic oscillation and therefore, avoidance of the boundary is required at all costs.

Stall flutter takes a combination of the blade structural non-linearity with the unsteady aerodynamics associated with detached flow. Stall flutter is typically a one degree of freedom oscillation, which is triggered by detached flow. It is often found in the torsional mode for propellers due to the introduction of a pitching moment as a result of a change in the local centre of pressure. It is the detached aerodynamics that introduces significant non-linearity to this phenomenon, and the trigger of it can be a result of the blade design or flight conditions. For example, typical conditions for a propeller to undergo stall flutter is during take-off. This is due to the requirement of high thrust/pitch during the acceleration phase of the aircraft. Due to the highly dynamic and, sometimes, sudden nature of stall, significant reductions to the overall propeller flutter boundary can be found due to this phenomenon. Such results are discussed and presented in this paper.

Although stall flutter is typically a transient and short-lived response, the development of modern propeller blades with high sweep/taper and thin aerofoil sections is likely to result in changes to the flutter boundary. Furthermore, a new generation of propellers is needed to power electrically driven aircraft, both fixed- and rotary-winged. In the latter category, there is a great variety of multi-rotor configurations over a range of scales and Reynolds numbers (including eVTOL), operating with variable angular speeds. As a result, a re-investigation of the stall flutter is required to ensure safety of the new designs and to also push the boundaries of the design.

A high number of stall flutter studies focus on aerofoil dynamic stall. In this contribution, a description is presented on dynamic stall; links to aerodynamic damping and propeller stall flutter boundaries are demonstrated. Details are then shown on the propeller stall flutter. This includes historic experiments and modern numerical studies. Stall flutter associated with turbomachinary and helicopter rotors is not included in this review, except within the context of the general stall flutter description.

2.0 Stall flutter description

One of the characteristics associated with helicopter rotors in skewed flight is the stall effects on the retreating side due to the combination of the angular and freestream velocities. Similar stall effects can also be found in a propeller during a high thrust/pitch setting such as take-off. If the blade exhibits elastic behaviour, there is a potential for further blade oscillations. This type of flutter is classed as stall flutter and is due to the combination of structural dynamics and non-linear aerodynamics. It remains a challenging problem to this day.

Stall on a rotor blade due to the control pitch oscillations is termed as dynamic stall and can be studied as a two-dimensional aerofoil problem. The inclusion of dynamic stall on an aerofoil introduces a unique set of characteristics not found in static aerodynamic stall. The use of quasi-static look-up tables for lift, drag and moment allows the calculation of rotor loads for hover, vertical flight or low-speed level flight. This method has been verified with rotor load studies and is justified due to the relatively low angles of attack. However, at high speeds such quasi-static tables cannot be used; investigations such as those of Kufield and Bousman [1] pointed out the difficulty when the airloads of the UH-60A helicopter rotor were measured in flight. The airloads indicate a significant deviation from the two-dimensional aerofoil results. To understand this discrepancy, an examination of both the static-stall and dynamic stall effects on an aerofoil is presented.

2.1 Aerofoil static-stall

The range of angle-of-attack is significant in rotorcraft performance. The lift, drag and moment coefficients are examined from attached flow to deep stall; the typical trend for each case is presented in Fig. 1. This data is adapted from the experiments of Critzos et al. who studied the NACA 0012 static coefficients between $0^{\circ}$ and $180^{\circ}$ at Reynolds numbers $Re = 1.8\times10^6$ and $0.5\times10^6$ [2].

Figure 1. Typical trends in lift, drag and moment coefficient during static stall [2].

A linear increase in lift is found during the initial stages and this remains until the static stall angle ( $\alpha_s$ ). This angle is defined based upon the location of the maximum lift coefficient ( $c_{\ell_{max}}$ ) and can often be accompanied by detached flow. Following this static stall angle, the lift coefficient decreases in value for a short range before gradually increasing towards a secondary peak. This secondary peak is typically under the presence of deep stall where three-dimensional detached flow effects introduce pressure reductions causing the increase in lift.

Up to the point of static stall, the drag and moment coefficients remain fairly constant. Following this point, both the drag and moment increase sharply; the pitching moment experiences a sharp nose-down effect resulting from a rearward shift in the centre of pressure. Both coefficients continue to increase with the angle-of-attack, with the drag reaching a maximum at $90^{\circ}$ , when the aerofoil is perpendicular to the flow.

2.2 Aerofoil dynamic-stall

A common feature among aerofoils undergoing light and deep dynamic stall is the presence of hysteresis, which is caused by a time lag between forcing and aerodynamic response. This hysteresis depends on several parameters, including frequency, angle-of-attack amplitude, Mach number, turbulence, etc. Typical trends are shown in Fig. 2 for the lift and moment coefficients. One of the effects of a dynamically oscillating aerofoil is the increase in maximum lift coefficient and stall angle, Fig. 2(a). This gives an indication as to the difference in airloads between the low-angle-of-attack look-up tables and what is observed in the true blade. An increase in maximum lift is potentially beneficial in terms of the blade operation, as momentary increases can be obtained at lower angles of attack. However, an increased load can potentially lead to an increase in blade torsion, which pushes the blade towards a potential stall flutter boundary due to the introduction of a larger pitching moment.

Figure 2. Typical hysteresis trends in lift and moment coefficient during dynamic stall [3].

Likewise, a pitching moment hysteresis loop can be obtained, Fig. 2(b). A similar increase in moment above the static stall is found; this is related to the increase in lift. The number of loops included within the profile depends on the angle at which the aerofoil is oscillating around and the amount of detached flow associated with the condition. These loops can categorise the stability of an aerofoil section via the determination of aerodynamic damping in pitch. Aerodynamic damping is calculated as the summation of areas of each individual loop, with a clockwise loop indicating negative damping and anti-clockwise loop producing positive damping.

2.3 Aerodynamic damping

Aerodynamic damping is the result of aerodynamic forces and moments exerted on to a structure. Aerodynamic damping often opposes structural damping, and can potentially result in aeroelastic instabilities. Aerodynamic damping is critical to the flutter characteristics of a structure. Above the flutter velocity, the work of the given fluid on a structure is said to be negatively damped. Thus, the structure’s oscillatory motions tend to increase with time unless an external factor such as a cyclic input adjusts the given pitch angle.

Stall flutter originates from separated flow and is found to be present in helicopter rotors, propellers and other rotating wings. A stall flutter instability can only be corrected via positive structural damping, change of stiffness or a change in the aerodynamic conditions. As a result, an investigation into stall flutter can begin from the aerodynamic damping of a system. Damping estimation is often performed for aerofoils and full three-dimensional calculations are rare.

Early investigations of stall flutter were two-dimensional and experimental [4–6]. These investigations focused on the determination of the aerodynamic coefficients during dynamic stall to improve helicopter rotor performance in forward flight. Such 2D investigations of oscillating aerofoils highlighted the trends seen during differing stall regimes. These regimes are highlighted in Fig. 3, where the pitching moment coefficient is shown for a pitching NACA-0012 aerofoil.

Figure 3. Pitching moment coefficient trends for each stall regime, as observed by McCroskey for a NACA 0012 aerofoil pitching at $\alpha = \alpha_o + 10^{\circ} \sin (2 k \tau)$ , where $k=0.10$ [6].

No stall : The aerofoil motion remains below the static stall angle and the use of quasi-steady aerodynamics is sufficient enough to predict the aerofoil loading. Both the lift and pitching moment coefficients are found to circle anti-clockwise with no crossing of the downstroke and upstroke profiles. Therefore, aerodynamic pitch damping is found to be positive.

Stall onset : The aerofoil motion reaches the static stall angle. There is often found a slight reduction within the area of the anti-clockwise loop, however, no crossing of the profiles are seen and therefore quasi-steady aerodynamics can be used to estimate the loads [7, 8].

Light stall : It is within this regime that dynamic stall vortices are present. These vortices are a combination of shedding and passage events that occur across the upper surface of the aerofoil. A typical, and clear, pattern often exists where there is a build-up of detached flow at the leading edge following an attached period. This then passes over the upper surface before shedding as it reaches the trailing edge. Often, secondary events can occur under the right conditions. The cycle then returns to the attached state. During light stall, the aerofoil motion reaches values higher than the static stall angle, with the aerofoil loads characterised by a hysteresis effect and it is the dynamic stall vortices which introduce the large hysteresis. The development of separated flow regions are found to be sensitive to the aerofoil geometry, freestream Reynolds number and Mach number, reduced frequency of the aerofoil oscillation, and the mean and harmonic angles of attack. It is also within the regime that there is the highest tendency towards negative aerodynamic damping.

Deep stall : The aerofoil motion is often found to pitch entirely beyond the static stall angle and it is where the strongest effects of the dynamic stall vortex are seen. The aerofoil loading is characterised by a strong hysteresis effect, with significantly larger peak load coefficients and multiple damping loops.

Following the initial propeller stall flutter studies, a transition was made from propeller driven aircraft to turbofan engines. This resulted in a reduction in propeller research with attention paid to the more immediate issue of helicopter dynamic stall. It is from this research that a significant amount of understanding of aerodynamic damping was gleaned. The above description of dynamic stall is focused on helicopter applications, and as such, involves large amplitude variations in pitch due to the cyclic inputs of a helicopter in forward flight. In propeller dynamic stall, significantly smaller oscillations in pitch are found around the static stall angle due to the lack of cyclic inputs, and therefore the mechanism which triggers the stall flutter phenomena must be taken from the context of aerodynamic damping.

Following some of the initial experimental propeller stall flutter investigations [9], attempts were made to better understand the stall flutter boundaries with respect to aerodynamic damping. This included the experiments of Fanti et al. [10] & Lemnios [11]. However, as described by Carta & Lorber [12], such estimations were seen to be inadequate for design and deeper understanding of the mechanism due to insensitive instruments. Carta & Lorber conducted an experimental investigation into the variation in aerodynamic damping due to small amplitude oscillations in pitch of a propeller aerofoil section [12]. This study followed the successful demonstration of the damping extraction for a cascade environment [13].

Experiments were performed in the United Technologies Research Centre (UTRC) wind tunnel. Sinusoidal oscillations of the Sikorsky SC1095 aerofoil section were performed over a range of amplitudes, mean angles and reduced frequencies. Pitch amplitudes ranged from $0.5^{\circ}$ to $4^{\circ}$ , with the mean angle between $8^{\circ}$ and $10.5^{\circ}$ and reduced frequencies between 0.05 and 0.16. From initial static results, it was shown that the static stall angle was $\sim9.5^{\circ}$ for this section. Pressure transducers were installed in the section in order to determine surface pressure and load profiles. The results indicate a significant influence of these parameters. The static derivative results were shown to be sensitive to the mean value when oscillating at low angles of attack for a range of reduced frequencies. The derivatives were found to vary significantly once the mean angle approached the static stall. This variation was seen to be random, with this random nature also present within the load time histories at low oscillating amplitudes [12]. However, at higher oscillation values such as those seen in helicopter rotors, reduced frequency becomes more influential with almost no trend captured due to changes in mean angle. Therefore, in the interest of propeller stall flutter which takes into account smaller oscillation values, the mean value sensitivity is critical to the test case in question.

Pitching moment hysteresis loops similar to Fig. 3 were presented, which determine the aerodynamic damping estimations. For this experiment, low-amplitude oscillations were found to induce a dominant clockwise loop when the mean angle was selected above the reference static stall angle. The trend of aerodynamic damping from these given profiles is presented in Fig. 4. For the highest oscillation value, the mean angle-of-attack is seen to have very little influence on the aerodynamic damping value with an almost constant negative value produced. However, the greatest variation in aerodynamic damping comes from the lower values with the destabilising range being dependent on the mean angle. An overall negative damping result can be mitigated via the blade design. However, the variation in damping seen from the $0.5^{\circ}$ and $2^{\circ}$ plots, is more difficult to design for; therefore, such lower amplitude profiles become more critical.

Figure 4. Aerodynamic damping trends across the mean and oscillating magnitude of the angle-of- attack [12].

Using the data gained from experimental investigations, such as Carta & Lorber [12] and McCroskey [6], validation of aerodynamic models could be undertaken with respect to dynamic stall and aerodynamic damping. This involved theoretical and semi-empirical models, with a review conducted by Ericsson & Reding [14]. However, with the need for experimental parameters to ensure the accuracy of semi-empirical formulations, a review was conducted by Ekaterinaris & Platzer [15], which focused on the use of the Navier-Stokes equations for aerofoil aerodynamic damping extraction. This was brought on by the development of the numerical methods and increased computational power to allow for Navier-Stokes methods to be used in dynamic stall modelling.

The use of Navier-Stokes-based modelling can not only assist in the extraction of damping estimations but also determine steady airloads for the semi-empirical formulations, should no experimental data exist. Such a procedure was conducted by Beedy et al. [16] who utilised the ONERA aerodynamic model developed by Tran and Petot [17] and Dat and Tran [18] in the 1980s, with CFD coefficients for the study of aerofoil stall flutter. Such a procedure was found to qualitively correlate well to experimental data [19] in terms of flutter velocity and frequency, with a comparison shown in Fig. 5. As observed, a rapid increase in frequency is found around 30m/s, correlating with the experiment. At this point, further work was needed to assess the ability of the turbulence model to predict the required ONERA coefficients; however the results were deemed to be satisfactory and thus highlighted the ability to use Navier-Stokes CFD with semi-empirical models should typical experimental values be unavailable.

Figure 5. Prediction of flutter velocity and frequency using the ONERA aerodynamic model with Navier-Stokes coefficients [16].

Figure 6. Comparison of the NACA 0012 pitching moment coefficient for the 2D and quasi-3D simulations to experimental data [20].

Ekaterinaris & Platzer [15] discussed the effects of numerical schemes, turbulence modelling and the effect of transition, with the turbulence and transition modelling seen to have a significant effect on the experimental correlation. They concluded that compressibility, transition, and flow reattachment were key flow effects limiting the correlation between the experiments and CFD results. Compressibility effects are driven by the partially supersonic profile observed across the aerofoil section. This is typically seen on the advancing side of the rotor due to the combination of freestream and rotational velocities combining. The use of transition models was seen to have a significant effect on the leading edge separation and, with the transition models of the time utilising empirical coefficients, an improvement was required to enhance the correlation.

The final factor was related to flow reattachment, and this was seen during the aerofoil downstroke where poor agreement was shown. This is often found in the light and deep dynamic stall regimes where the influence of the three-dimensional detached flow effects become important. Two-dimensional CFD cannot capture this effect, and therefore, an increase to the computational expense is required. This comes in the form of three-dimensional calculations and/or resolved flow modelling. Resolved flow modelling allows for the shedding process associated with dynamic stall process to be captured, potentially improving the downstroke performance.

In order to improve the correlation and mitigate the effects of flow reattachment, an investigation was conducted by Higgins et al. who utilised resolved flow modelling across aerofoil sections undergoing deep dynamic stall [20]. Quasi-3D simulations of the NACA 0012, $70\%R$ and $90\%R$ Commander blade sections were conducted. The simulations were termed quasi-3D due to the use of periodic boundary conditions where typical two-dimensional planes are applied. This allowed for the use of scale-resolving techniques to be used, and for this investigation the Scale-Adaptive Simulation (SAS) method was employed. The SAS model is a hybrid URANS-Large Eddy Simulation method used to resolve the stalled flow [21]. For the SAS formulation, an additional source term is added to the $k-\omega$ Shear Stress Transport (SST) equations, which allows for the local adjustment to the von-Karman length scale and balances the contributions from resolved and statistical components. Using this method, an improvement on the NACA 0012 experimental load correlation and aerodynamic damping estimation was found against the two-dimensional URANS results. This was driven by the ability of the method of capture the shedding of vorticity during the heavily detached downstroke/reattachment phase. This emphasises the potential need for high-fidelity aerodynamic modelling during detached flow conditions.

Furthermore, many studies have been conducted in recent times which focused on two-dimensional aerofoils. Such investigations made use of pitch and plunge rigid models [22, 23, 24], cyber-physical models (where the structural response of the aerofoil is given based upon a known set of derived equations of motion) [25, 26], or if numerical, used dynamic stall based aerodynamic models [27, 28]. The use of such studies can provide fundamental analysis into aerofoil stall response, with conservative boundaries when applied to three-dimensional test cases.

3.0 Propeller stall flutter

The topic of propeller flutter is split into three components, and this includes bending-torsion, stall and whirl flutter. These types of flutter, particularly whirl and stall, require detailed modelling of the aerodynamics and structural response of a propeller. For stall flutter, the non-linear aerodynamics is a result of the detached flow-field, which triggers the aeroelastic excitation. Successful capture of such aerodynamic interactions allows for increased accuracy in surface loads, and increased accuracy in predicting the resultant flutter boundary.

3.1 Experimental studies

Propeller stall flutter research began during World War II and was primarily motivated by the use of propeller-driven aircraft such as the Supermarine Spitfire and North American P-51 Mustang. This research began with two studies: one from the USA and another by the UK. The UK study was conducted by Sterne and focused on the flutter analysis of the 4-bladed Spitfire propeller [29]. The USA investigation was conducted by Theodorsen & Regier and used a range of wind-tunnel designed propellers [30].

Focusing on the UK study [29], static experiments were conducted to determine the flutter boundary of the propeller, with the boundary presented in the form of the blade pitch versus blade rotational velocity, Fig. 7. The blade operated at 19 different pitch angles ranging from $8^{\circ}$ to $32^{\circ}$ (the pitch was measured at the $70\%R$ station). Note that a reference pitch of $8^{\circ}$ at $70\%R$ results in a $0^{\circ}$ pitch angle at the blade tip. As observed from Figure 7(a), the experiment found that the blade would flutter at a much lower velocity within a specific pitch region. This region ranged from $20^{\circ}$ to $28^{\circ}$ with a 24% reduction in flutter velocity found.

Figure 7. Spitfire and Firefly propeller flutter boundaries [29].

Sterne concluded that the distinct decrease in velocity between $20^{\circ}$ and $28^{\circ}$ was the result of the presence of stall, hence stall flutter oscillations. The remaining pitch angles were judged to be a result of classical bending-torsion flutter. This observation was also found by the work of Stüder who found a similar drop in velocity across a model wing due to stall and theorised it to be a result of negative aerodynamic damping [31].

In addition to this spin test on the Spitfire propeller, Sterne conducted an experiment to investigate the influence of the number of blades on the flutter characteristics. Those tested included a 2-, 3- and 4-bladed Firefly propeller. The Fairey Firefly blade was constructed in the same manner as the Spitfire, using compressed wood, with a slightly larger diameter of 13 feet compared to 10 feet 9 in. A comparison of the thickness/chord ratio and the solidity for the 4-bladed cases is presented in Fig. 8.

Figure 8. Spitfire and Firefly blade definitions [29].

From the results in Fig. 7 it was implied that the number of blades have an influence on the flutter boundary. In the classical flutter boundary range, the 2- and 4-bladed propellers produced similar results with the 3-bladed propeller operating at an angular velocity $\sim90 \mbox{rpm} $ higher. One difference in the experimental setup between the 2- or 4-bladed propeller and the 3-bladed was a slightly different hub. It was expected that such a difference would have an influence on the flutter boundary and, hence, the increase within the classical range could be explained via this difference. In terms of the stall flutter range, a significant reduction in stall flutter boundary was observed in the 3-bladed propeller compared to the 2- and 4-bladed. Sterne concluded that the reduction in stall velocity could be a result of the greater influence of the preceding blade stall upon the oncoming blade. With the experiment conducted in static conditions, and no thrust data known, the blade wake could very well interact with a subsequent blade. Additionally, the influence of the propeller induced velocity is unknown due to the lack of thrust data and therefore, this propeller could be operating on the edge of a wake interaction that is felt for the 3-bladed propeller but not the 4-bladed due to the expected increase in induced velocity (assuming all other parameters remain the same).

The experimental investigation of Theodorsen & Regier (the USA wartime report on propeller flutter) [30] included a range of experiments on model propellers for use, and examination, in the wind-tunnel. This investigation was conducted in connection with the construction of the Langley, Ames and Cleveland blades for use within the respective tunnels. In this study, two propeller blades were presented. As in the case of the Spitfire/Firefly blades, both were made of laminated wood (spruce) with flat-bottom Clark Y aerofoil sections. The first propeller (Propeller A) was designed as a 6-bladed propeller with a 45-inch diameter. The second propeller (Propeller B) was designed as single bladed with a reduction in the chord length and thickness/chord ratio of one third. This reduction in size was to reduce the flutter velocity in order to examine the results.

The experiments were conducted in an open wind tunnel in static conditions, with only the induced velocity present at the tunnel exit. For Propeller B, this induced velocity was enhanced by a booster fan at the rear of the motor due to being of single blade design. A stroboscope was used at the blade tip in order to measure the bending and torsion, with strain gauges also used. With the blade pitch fixed, the rotational velocity was increased until flutter was determined for a range of propeller lift coefficients. These lift coefficients were adjusted through the change in the wind tunnel exit area, i.e. essentially compressing or expanding the wake. The change in lift coefficient was confirmed via pitot tubes.

From the experiment it was observed that at the lowest flutter velocities, the oscillations were found to be pure torsional fluctuations. This was determined from the strain gauges of Propeller B. It was concluded from this experiment, as in the Sterne experiments, that the stall flutter velocities are much lower than the classical boundary. In order to mitigate this effect, Theodorsen & Regier concluded that the design blade angle should be such that the lift coefficient produces zero blade twist. Equations are presented to determine such angles, with the effect of addition twist seen as inducing further detached flow. Approaching the classical flutter boundary, less and less detached flow was required to excite the blade.

Although both studies provide a good insight into propeller stall flutter, particularly the influence of stall on the propeller flutter region, the construction of the blades make it difficult to study numerically. In comparison to modern manufacturing processes, the compressed wood construction has the potential to introduce manufacturing inconsistencies and, therefore, make the understanding of the blade structural properties potentially difficult. From this basis, the Spitfire, Firefly and NASA blades could not be used for a numerical study.

Following this observation of classical and stalled flutter within the working range of the wooden Spitfire and Firefly propellers, Ewing et al. conducted a experiment on the duralumin blades of the Hawker Tempest 5-bladed aircraft [32]. This experiment was conducted in the spinning tower of the Royal Aircraft Establishment, with strain gauges applied to the blades to measure the vibratory stress. Each propeller was pitched below and above the stalling region with only stall flutter found during the runs.

The experiments found a heavy dependence of the stress results on the phase differences between subsequent blades. This phase difference ( $\psi$ ) was given by $\psi = 2 \pi p / Q$ , where Q is the number of blades and $p = 0, 1, 2 \ldots$ . As p changed value, there was a change in the phase relationship and wave-form. The experiments found the value of p to increase with blade rotational velocity at a given pitch angle. Due to this influence of phase, an exact comparison of the stress values for the number of blades could not be made, hence general stresses are compared. For the 3-, 4- and 5-bladed propellers, the greatest stress values occurred at $24^{\circ}$ . This reduced to $20^{\circ}$ for the 2-bladed case. Either side of $24^{\circ}$ , the 4- and 5-bladed propeller stress values drop rapidly with a progressively lower gradient for the 3- and 2-bladed propellers. From this experiment, it was observed that the examination of the number of blades of the same design could only provide a rough indication of stress level due to the influence of phase. For example, the velocity at which the peak stress occurred may coincide with a phase change and thus be at a lower stress level than the other.

In 1955, [9] Baker conducted a static experimental investigation into the effects of various parameters on the flutter boundary of a single model propeller blade to determine the minimum stall flutter condition. These parameters included the effects of the blade structure, the blade geometry, and the freestream flow conditions, with the key list and range presented in Table 1. Due to the change in certain parameters, a range of blades were manufactured. For the construction, a 16-series NACA aerofoil was used with the material selected as aluminum alloy, steel or maple, depending on the model and test case. All propeller blades were tested across the full pitch angle range.

Table 1. Range of parameters tested in the experimental investigation of Baker [9].

The experiments were conducted in the Langley vacuum sphere to allow for the use of Freon-12 to achieve more realistic Reynolds numbers during the density and Mach number studies. Strain gauges were attached to the blades to measure the bending and torsion oscillations with the rotational velocity gradually increased until flutter was observed.

The experimental results showed very little difference in the boundaries for the torsional stiffness, blade taper, blade twist, length/chord ratio and density studies. Substantial changes were found for the centre-of-gravity, sweep, thickness and Mach number runs with the non-dimensional flutter velocity results presented in Fig. 9 for each. The models used for the results shown are described in Table 2.

Figure 9. The changes in the flutter boundaries based upon the experiments of Baker [9].

Table 2. Range of model propeller design parameters for the investigation of Baker [9].

^* At tip, polynomial sweep

The most significant change in the stall flutter boundary was found in the centre-of-gravity results (Fig. 9(a)), where a 14.5% shift in CG results in a 44% reduction in the minimum flutter boundary. Little difference in the thickness ratio results (Fig. 9(b)) was found between the 3% and 6% plots with the 9% model resulting in a $\sim25\%$ increase in the boundary. Although a benefit was found using a thicker aerofoil section, propeller performance dictates thin aerofoils.

The introduction of sweep (Fig. 9(c)) to the blade design has a positive influence on the stall flutter boundary, with all sweep configurations seeing an increase in the flutter velocity of greater than 25%. Examining the results from a pitch angle of $10^{\circ}$ to $35^{\circ}$ , the non-linear sweep has the greatest continuous increase in the boundary across this region. The final observation from the investigation was the increase in flutter boundary for Mach ${\rm M}_a > 0.6$ (Fig. 9(d)). The Mach number influence was tested across two models at different pitch angles. The second test (Model 6), was seen to have a more gradual increase in the boundary, with the first (Model 5) having a steeper gradient above 0.6.

Baker did note that the resultant flutter boundaries were not within the normal operating cruise conditions for typical propellers, and it was only during take-off like conditions that the blades were found to suffer from stall flutter. Following the experiment of Baker, greater insight into several aerodynamic and structural influences on the stall flutter boundary were found due to this investigation. Although some details are available for the blade used (Table 2), greater detail would be required to map the blade inertias along the radius for a numerical study.

In 1956, Hubbard et al. [33] followed on from the work conducted by Baker by focusing on three specific parameters. The primary aim was to determine the effect of Mach number, structural damping and built-in twist on the flutter boundary. For these parameters, four single-bladed model propellers were constructed (one each for the Mach number and twist study with the final two for the structural damping). Each blade was defined in the same way as the Baker propellers [9] with Table 3 presenting the blade descriptions. Model 3B was constructed using laminated steel in order to change the internal structural damping. In the same manner as the Baker experiments, the propeller models were tested in the Langley vacuum sphere in static conditions. During the experiment, strain gauges were used to measure the vibratory stress levels with runs conducted at a fixed pitch angle, only varying the rotational velocity.

Table 3. Description of each experimental model propeller for the investigation of Hubbard [33].

Figure 10 shows the experimental results for the twisted blade, Mach number and structural damping studies. Focusing on the twisted blade result (Fig. 10(a)) and using the strain gauges, the type of blade oscillations across the examined pitch range was determined. Through the low pitch range ( $-8^{\circ}$ to $8^{\circ}$ ), the blade oscillations were dominated by bending-torsion or pure bending oscillations. The highest flutter velocities were found during this region with an average value 228% above the minimum non-dimensional velocity of 1.26. The reduction of the pitch angle to higher-magnitude negative angles resulted in the activation of the torsional mode, reducing the flutter velocity. At high positive angles, the traditional shape expected from the Baker [9] and Sterne [29] experiments was found. Due to the detached flow at these angles, the flutter velocity was reduced significantly compared to the classical flutter bounds due to the active torsional mode.

Figure 10. The changes found in the flutter boundaries based upon the experiments of Hubbard [33].

The effect of Mach number (Fig. 10(b)), where the effect of Mach number is depicted by values of sound-speed coefficient which is defined as $b \omega_{\alpha} / a$ with b being the semi-chord, $\omega_{\alpha}$ the natural torsional circular frequency and a the speed of sound, on the propeller flutter velocity mirrored that found by Baker [9]. As the speed of sound was decreased, thus increasing the Mach number, the extent of the flutter region was reduced. This is likely down to a compressibility effect which limits maximum lift and delays the stall onset.

In terms of structural damping (Fig. 10(c)), an increase in the internal damping was seen to increase the flutter velocity. This was to be expected and therefore any oscillations due to a negative aerodynamic damping were mitigated by the increase in structural damping.

The combined works of Baker [9] and Hubbard et al. [33] allowed for a greater understanding of propeller stall flutter due to the volume of parameters investigated. However, the blades used in these experiments do not represent a realistic modern blade and this has a significant effect on the observed flutter boundary.

Similarly, Allis and Swihart conducted a static experimental investigation to determine the effect of blade-section camber [34]. This work followed on from the experiments of Wood & Swihart [35] who found an increase in stall flutter boundary due to an increase in camber at one specific pitch angle. This pitch angle was extended to test over a range from $16^{\circ}$ to $38^{\circ}$ . Three 2-bladed NACA propellers, of varying camber, were tested in the Langley dynamometer with the blades fitted with strain gauges. Figure 11 presents the installed blades on the dynamometer. At a fixed pitch angle for the given blade, the propeller rotational velocity was slowly increased until flutter was encountered.

Figure 11. Installed NACA propeller on the Langley dynamometer for the experimental study of Allis and Swihart [34].

Presented in Fig. 12 is the non-dimensional flutter velocity and thrust coefficient results for the range of cambered propellers. Focusing on the flutter boundary result (Fig. 12(a)), it was shown that at $16^{\circ}$ an increase in flutter velocity of 45% was achieved with the highest cambered propeller compared to the symmetrical blades. This difference slowly reduces with increased pitch angle until $30^{\circ}$ where the flutter velocity values are approximately the same. Above $30^{\circ}$ , an increase in the flutter velocity was seen by all, with the symmetrical propeller having the highest value at the largest pitch angle. Allis and Swihart observed that at $16^{\circ}$ , only the inboard radial stages were stalled and this increased to the entire blade at $30^{\circ}$ .

Figure 12. Flutter boundary and thrust coefficient changes for each cambered propeller blade [34].

One of the key observations from the experiment was the significant difference in the flutter boundary at $16^{\circ}$ . The flutter boundary was not only driven by the combination of detached flow and structural damping, but also the total airloads applied to the blade. At higher thrust, higher fluctuations in blade deformation was expected. This is what was observed at the $16^{\circ}$ pitch angle, with this confirmed in Fig. 12(b) where the highest camber blade was shown to have a 36% increase in thrust in comparison to the symmetrical blade. With greater thrust comes greater induced velocity and therefore a reduced effective angle-of-attack at a given rotational velocity. In order to assess the true effect of blade-section camber, a trimming process would be required to ensure the same thrust was generated from the blade.

Rogallo & Yaggy conducted an experimental campaign on a 3-bladed 10-foot propeller to determine the effect of positive/negative thrust and thrust axis inclination on the stall flutter characteristics [36]. The experiments were conducted in the 40-by-80 foot NASA Ames wind tunnel facility, with the installed blades presented in Fig. 13(a). The blades were instrumented with strain gauges to determine the torsional and bending stresses. The positioning of the strain gauges, and definition of the blade profiles, shown in Fig. 13(b). These gauges were accompanied by flow-field rakes upstream and downstream of the propeller, with the rakes shown in Fig. 13(a). The use of the NASA Ames wind tunnel facility allowed for the assessment of freestream velocity, the first published dataset to take into account its effect.

Figure 13. Propeller blade installation and instrumentation for the stall flutter study of Rogallo & Yaggy [36].

In order to capture positive and negative thrust conditions, the blade pitch angle was varied throughout the experiment from $-17.5 ^{\circ}$ to $+42^{\circ}$ . The exact range was dependent on the desired condition and also the freestream velocity. The freestream velocity was varied from zero to 183.5 ft/s. This maximum velocity represents 50% of the maximum achieved tip velocity. For the majority of the tests, the rotational velocity was fixed with the blade pitch angle varied until stall flutter was encountered. When examining the effect of the thrust axis inclination (yaw), the positive thrust pitch range was examined along with the inclusion of a freestream velocity. The axis angle was varied by $10^{\circ}$ and $15^{\circ}$ , with only the $15^{\circ}$ yaw at the highest velocity presented.

Presented in Fig. 14 is a selection of the stall flutter results for the Rogallo & Yaggy experimental investigation [36]. A selection using different formats is presented to highlight the key observations of Rogallo & Yaggy. For the positive thrust tests, similar trends were found in the torsional stress results with varying freestream velocity. To highlight the sharp increase in stress observed at the flutter boundary, the results in static conditions are shown in Fig. 14(a). With increasing rotational velocity, lower pitch angles are required before entering stall flutter. For rotational velocities between 1,800 and 1,200rpm, a sharp rise in stress was captured at a given pitch with very few additional angles tested before reaching the stain gauge limits. This sharp rise trend transitioned at 1,000rpm, in static conditions, to a more gradual gradient with the lower velocity runs able to reach the full pitch angle range.

Figure 14. Stall flutter results for the experimental investigation of Rogallo & Yaggy [36].

This transition velocity, from very high to lower stress gradients, reduced with increasing freestream. The overall effect of increasing freestream is presented in Fig. 14(b) where the stall flutter boundary in terms of pitch and rotational velocity is shown for the positive thrust cases. An increasing freestream velocity pushed the stall flutter boundary further, i.e. a greater pitch and rotational velocity was required before encountering stall flutter.

In terms of the negative thrust study, similar trends and values were found for all tested freestream velocities. Unlike the positive thrust cases, no coupling between the pitch angle and rotational velocity was found at a given freestream, with increases in stress observed at similar pitch angles. With the sharp rise in shear stress occurring around $0^{\circ}$ , the stress values remained fairly stable before slowly decreasing at the lowest pitch angles.

One of the final investigations was the effect of axis yaw on the stall flutter boundary. This effect is shown in Fig. 14(c) where the highest freestream velocity results are presented for a yaw angle of $15^{\circ}$ over a range of rotational velocities. The introduction of a yaw angle had little effect on the shear stress results. To understand, the thrust and torque coefficients were obtained along the blade radius from the rake data. The thrust values indicated only the outer section of the blade (20%) to be stalled with a small magnitude found in the change of angle-of-attack. With the change in angle-of-attack, driven from the axis yaw, too small, the detached flow dominated resulting in similar stress results.

In order to understand the trend between generated thrust and the stall flutter boundary, Rogallo & Yaggy extracted the thrust coefficient results along the blade radius and compared them to the values indicated from the shear stress plots when above $\pm 1,000\ psi$ . In Fig. 14(d) are the thrust coefficient results for the $90\%R$ radial station at the highest freestream velocity in positive thrust. It was shown that the stall flutter boundary appears to coalesce with the drop-off in thrust coefficient, with this observation clearest at $90\%R$ compared to inner stations. This gives further indication to the dominance of the separated flow.

A similar investigation was conducted by Dowty Propellers in 1979 at the spinning tower of the Royal Aircraft Establishment, Farnborough [37]. The aim was to determine the torsional stress levels, and subsequently the flutter characteristics, of the Aero Twin Commander three-bladed propeller when spun at fixed pitch over a range of rotational velocities. The torsional stress levels were measured via strain gauges placed along the propeller radius.

During the first stage of the test, high levels of torsional stress were seen for a given pitch angle, with the test terminated due to excessive oscillations. Upon examination of the propeller apparatus, the propeller crosshead had failed and hence a redesign was required. Following, overspeed and 30 minute power run tests were completed at lower blade pitch angles, without further attempts to probe the stall flutter boundary and therefore further trends in terms of the stall flutter boundary could not be obtained. From the torsional stress results, a clear stall flutter boundary was observed with a sharp increase in stress seen around 1,650rpm, in a similar manner to the plots of Rogallo & Yaggy (Fig. 14(a)).

Smith conducted a static experimental investigation into three prop-fan model propellers [38]. These models were designated SR-2, SR-3 and SR-5, with the blades featuring increasing levels of sweep from the unswept, SR-2, design to the highly swept, SR-5. Additionally, different blade numbers are found across these test cases with the SR-2 and SR-3 having eight blades and ten for the SR-5. No attempt was made to understand the effect of these changes on the flutter response. As in the Dowty tests, for a fixed pitch angle, the propeller rotational velocity was increased to maximum before returning to its baseline value. The stress levels were measured via strain gauges. The positioning of the strain gauges and the blade profiles are presented in Fig. 15. This study was motivated by the transition to thin and highly swept blade designs, conclusions that were drawn from past studies, with the aim to ensure the structural stability of the blades.

Figure 15. Strain gauge positioning and SR propeller blade definitions for the experiments of Smith [38].

Each propeller was installed on an isolated nacelle in the UTRC large subsonic wind tunnel where previous in-house tests had been conducted [39]. The SR-2 propeller was an 8-bladed rotor constructed of steel. The SR-3 and SR-5 were both made using aluminum with the SR-5 having an additional two blades. All configurations were derived from the full scale designs at an intended operating Mach ${\rm M}_a = 0.8$ . A reduction to $1/8^{th}$ scale was required due to the wind tunnel constraints.

The total vibratory stress results were compared against previous experiments of Smith in the NASA-Lewis Research Center [40]. The additional tests of Smith [40] not only looked into static conditions, but also introduced a forward velocity and yaw angle to the freestream to determine its effect. Focusing on the static results, plots of the total vibratory shear stress are presented in Fig. 16 for each blade. Outputs were also presented for the mid and tip blade bending, with shear only presented here due to stall flutter typically being active in the torsional mode.

Figure 16. Total vibratory shear stress results for the two sets of static experiments by Smith [38, 40].

For each propeller blade, results are presented for a pitch angle of $32^{\circ}$ and a subsequent higher value. This higher value was selected as $40^{\circ}$ , $60^{\circ}$ and $50^{\circ}$ for the SR-2, SR-3, and SR-5 propellers, respectively. Focusing on the comparison between the experiments (Fig. 16(a–c)), it was shown by Smith [38] that good correlation was found between the UTRC and NASA-Lewis investigations. The UTRC results at the critical speeds were found to be higher in all blades with this theorised to be the effect of turbulence.

Figure 16(d) presents the results for the $32^{\circ}$ pitch runs only for each propeller. It was observed that during the lower velocity range, similar trend and gradients were seen between the blades. The results start to deviate around 6,000rpm. At this stage, the SR-2 propeller was stopped due to excessive vibration with the maximum speed of 9,000rpm not reached. At this point, the SR-2 propeller was seen to have the highest stress values. As the stress levels are lower during this stage for the SR-3 and SR-5 propellers, the runs continue towards the maximum velocity. Upon reaching the maximum velocity, the SR-3 propeller was seen to have a greater stress value, with the SR-5 results dropping following 7,000rpm. This overall trend mirrors previous experiments of Smith [40] in that an increase in blade sweep results in a decrease in the stall flutter boundary. However, a trend which was not previously captured was the excitation of a bending mode during a stall flutter response and this was found, particularly, for the SR-3 propeller. It was also concluded that the SR-5 oscillated in the second bending mode, however, the excitement is found to be limited and therefore any such conclusion may be difficult to correlate. This response is different from expected responses in that it is usually a torsional mode that is triggered and this was found for the SR-2. Upon examination of the results, Smith found that the high stress regions were found around a pitch angle of $50^{\circ}$ with a rotational velocity of $6000 $ rpm ( $195.4{\rm m/s}$ tip velocity) with these values way beyond stall. Hence, the response at these angles may be that of a buffeting response.

Along with the static experiments, Smith conducted a similar range of tests on the SR-2, SR-3 and SR-5 propellers in low forward-flight conditions in axial and yawed conditions [40]. A maximum forward-flight ${\rm M}_a = 0.35$ was selected with a yaw range of $-2^{\circ}$ to $+15^{\circ}$ . In Fig. 17 is the total stress results for the Mach number and yawed inflow effects on the SR-3 propeller. In Fig. 17(a) the total vibratory stress results are presented with Fig. 17(b) showing the total stress per excitation factor (termed stress sensitivity). The excitation factor was used by Smith due to its linear dependence on yaw and was defined as:

(1) \begin{equation}EF = \psi_i \left(V_{eq}/348\right)^2\end{equation}

where $\psi_i$ is the inflow angle and $V_{eq}$ the equivalent airspeed. Focusing on the inflow effect (Fig. 17(a)) and the results presented are for a pitch angle of $35^{\circ}$ over a range of inflow velocities. A significant decrease in the stress values was seen with the introduction of an inflow velocity. The reduction in stress was mainly seen for the lowest forward-flight speed ( ${\rm M}_a = 0.1$ ), with similar trends observed during the higher velocities. A comparable profile was also found for the SR-2 and SR-5 blades with these profiles presented in the experimental article. Along with these trends, an increase in stress was observed with an increase in blade pitch. This result pushed the propeller further towards detached flow, and was observed across all tested Mach numbers.

Figure 17. Stress results for the forward-flight and yawed experiments by Smith [40].

The secondary factor analysed by Smith was the introduction of an angle to the inflow velocity. Furthermore, Smith presented the results for the $15^{\circ}$ yawed inflow for all blades, with the SR-3 at ${\rm M}_a = 0.1$ , across a range of blade pitch angles, presented in Fig. 17(b). The introduction of a yawed inflow resulted in a 1/revolution fluctuation in the stress. This was seen as the main component and therefore resulted in a linear increase in stress with yawed inflow angle. A resonance-like effect was found when the 1/revolution frequency equaled the rotational velocity, resulting in a spike in stress. Although not directly related to stall flutter, this highlights the influence of the yaw effect on the blade vibrations. This spike was mitigated with the introduction of greater sweep. As observed from Fig. 17(b), the stress sensitivity increased with increasing blade pitch and rotational velocity.

Due to the amount of data available and the publicly available blade definitions, they remain a viable option for future stall flutter analysis. The structural properties would need to be defined with a clear flutter boundary not present within the test range. Smith also acknowledged the effect of the highly-stiff titanium model propellers due to being of single material construction used for this study. The full discrepancy between model- and full-scale would need to be examined.

Presented in Table 4 is a summary of the experimental datasets currently available. The table includes some of the pros and cons highlighted within this review along with some of the conditions examined. It is clear from this section, and table, that no ideal test case is currently available and a factor which may have to be revisited, especially with the introduction of the eVTOL market.

Table 4. Summary of propeller stall flutter experimental test cases.

3.2 Numerical studies

In addition to the experiments, Smith [38, 40] conducted numerical simulations of the SR propellers with a comparison made to the wind tunnel results. Both investigations used the F203 stability analysis code developed by Turnberg [41]. Although initially derived for classical flutter, through the aerodynamic update it was used to analyse for stall flutter. This code required the supply of the blade mode shapes and frequencies, along with the aerodynamic loads, to determine an eigenvalue solution.

The mode shapes and frequencies were derived via a finite element analysis, with a beam based method used for the SR-2 propeller and NASTRAN for the SR-3 and SR-5 blades. The beam method was limited to long, slender, isotropic blades. Due to the sweep, offset and large chord lengths of the SR-3 and SR-5 blades, the beam models would fail and, hence, finite element grids were derived in NASTRAN to utilise its non-linear static analysis method. The grids for these blades used a combination of quadrilateral and triangular elements. These results were then transformed into the required coordinate system. The transformation splits the displacements into radial sections with two translational and one rotational degree of freedom. This was done in order to linearise the aerodynamics during the stability analysis. For the unsteady aerodynamics, a Goldstein-type [42] performance strip analysis was used to determine the steady-static sectional lift and moment curves.

Figure 18(a) shows the estimated stall flutter predictions of the static experiments of Smith [38], with the prediction assumed to be at the point where the damping goes negative for the third mode of type torsion, where the only instability was found. As shown, the stall flutter predictions mirror the experiments in that the introduction of greater sweep was seen to increase the stall flutter boundary or reduce torsional stress. At this current condition (7,000rpm), the highly swept SR-5 propeller does not suffer from stall flutter and requires the velocity to be increased towards the maximum value to induce the instability.

Figure 18. Estimated stall flutter predictions for the static and forward-flight experiments of Smith [38, 40].

The forward-flight results (Fig. 18(b)) predicted the same trend as the experiments [40] where the difference of Mach number indicates a change to the freestream velocity with constant speed of sound. With increasing freestream velocity, higher blade pitch angles and velocities were required to induce stall flutter. Although the trend was captured through the numerical method, a poorer correlation was found in terms of the exact boundary. This will be primarily a result of the aerodynamic modelling and, hence, greater effort is required to model the unsteady flow.

Based upon the experimental results of Smith [38], Reddy & Kaza conducted a numerical investigation into the SR-2 propeller [43]. They sought to improve on the numerical methods of Smith [38, 40], in particular the aerodynamic modelling. As was described in Section 2.3, the instability of stall flutter can be derived from the understanding of dynamic stall and such a study was conducted by Reddy & Kaza in 1987 [44]. This study compared three semi-empirical dynamic stall models in terms of their lift and moment hysteresis loop prediction for a sinusoidally oscillating aerofoil.

The three semi-empirical models include the Gormont model [45, 46], the Gangwani model [47, 48] and the ONERA model [49]. The Gormont model required only one empirical parameter, with a number of stall parameters required for the Gangwani and ONERA models. A short description is given for each modelling method:

Gormont model : This model uses a similar technique as per the Leishman-Beddoes model [50] in the form of an effective angle-of-attack. The lift and moment coefficient are then obtained from static aerofoil data. This correction for the effective angle-of-attack is a function of the rate and is based on oscillating aerofoil test data. The correlation between the true and effective angles is given by:

(2) \begin{equation}\alpha_E = \alpha - K_1 \Delta \alpha_{DS},\end{equation}

where $\Delta \alpha_{DS}$ is an incremental dynamic stall angle given by:

(3) \begin{equation}\Delta \alpha_{DS} = \gamma \sqrt{ \left| \frac{\dot{\alpha} c}{2 U} \right|}.\end{equation}

Here U is the resultant velocity, c the aerofoil chord, $\dot{\alpha}$ the rate of angle-of-attack and $\gamma_g$ the non-dimensional empirical function derived from experimental tests. This empirical parameter is determined by:

(4) \begin{equation}\gamma_g = \sqrt{\mid Ar \mid},\end{equation}

where, Ar is the non-dimensional rate of angle-of-attack. This parameter is a function of Mach number and thickness/chord. The final parameter $K_1$ is given by:

(5) \begin{equation}K_1 = 0.75 + 0.25 \sin (\dot{\alpha})\end{equation}

Gangwani model : This synthesised approach was development by Gangwani [47, 48]. The model utilises semi-empirical expressions representing the qualitative physical features of dynamic stall. The model is seen to capture the formation and streamwise travel of the dynamic stall vortex.

In this model, three stages of dynamic stall are defined: the stall onset; the trailing edge vortex; and the reattachment. Each of these stages are related empirically by the static stall angle, $\alpha_s$ , the non-dimensional rate of angle-of-attack, Ar, and the unsteady decay parameter, $\alpha_W$ . This decay parameter was previously defined as part of the findings of Beddoes [51] using a compressibility correction for the approximated Wagner function [52] (Equation (6)).

(6) \begin{equation}\alpha_W = 0.5 \alpha_{QS} - \int_0^s \alpha_{QS} (\sigma) \beta \phi^{\prime}_c (s-\sigma) d \sigma\end{equation}

Stall onset: The angle at which the dynamic stall vortex is initiated is given based upon the instance of moment stall (since moment stall occurs before lift) and this is represented as:

(7) \begin{equation}\alpha_{Dm} = \left( 1 + \epsilon + C_{Am} A_{m} + C_{Wm} \alpha_{Wm} \right) \alpha_s,\end{equation}

where the empirical parameters $\epsilon$ , $C_{Am}$ and $C_{Wm}$ are obtained from curve fitted experimental data.

Trailing edge vortex: Following the stall onset, the moment coefficient increases significantly in magnitude and this increase is a result of the vortex travel across the aerofoil. The point at which the moment coefficient reaches its maximum value is seen as the point at which the vortex reaches the trailing edge of the aerofoil. The time at which the vortex detaches from the trailing edge is given by:

(8) \begin{equation}s_{mt} = \frac{1}{C_{At} A_{m} + C_{\alpha t} \alpha_{Dm}},\end{equation}

where $C_{At}$ and $C_{\alpha t}$ are the empirical parameters defined using experiments.

Reattachment: The point at which the flow reattaches to the aerofoil depends on the freestream Mach number. For ${\rm M}_a < 0.4$ , the reattachment occurs at an angle which is less than the static stall angle. When ${\rm M}_a > 0.4$ , the reattachment angle occurs at an angle greater than the static stall. Using the empirical parameters $\epsilon$ , $C_{AR}$ and $C_{wR}$ , this reattachment angle ( $\alpha_{re}$ ) is defined as:

(9) \begin{equation}\alpha_{re} = \left( 1 - \epsilon + C_{AR} A_{m} + C_{wR} \alpha_{Wm} \right) \alpha_s\end{equation}

Using these equations for the dynamic stall events, the lift and moment coefficients can be derived alongside additional empirical parameters.

ONERA model : The ONERA model defines the lift and moment coefficients as a set of ordinary differential equations [49]. When attached, the equations are first order with this increasing to third order for the detached flow regime. A total of 11 empirical parameters are required for the differential equations (6 for lift and 5 for moment) with both lift and moment containing a linear ( $C_1$ ) and non-linear ( $C_2$ ) expression resulting in the following simple equations:

(10) \begin{equation}\begin{split}c_{\ell} &= C_{L1} + C_{L2}, \\c_{m,0.25c} &= C_{M1} + C_{M2},\end{split}\end{equation}

where each linear and non-linear expression is given as follows:

(11) \begin{equation}\begin{split}\dot{C_{L1}} + \lambda_{on} C_{L1} &= \lambda_{on} C_{L\ell} + \lambda_{on} s_{on} \dot{\theta} + \sigma_{on} \dot{\alpha} + s \ddot{\theta} \\\ddot{C_{L2}} + a_{on} \dot{C_{L2}} + r_{on} C_{L2} &= - \left( r_{on} \Delta C_L + E_{on} \dot{\alpha} \right) \\C_{M1} &= C_{M\ell} + s_{on} \dot{\theta} + \sigma_{on} \dot{\alpha} + s_{on} \ddot{\theta} \\\ddot{C_{M2}} + a_{on} \dot{C_{M2}} + r_{on} C_{M2} &= - \left( r_{on} \Delta C_M + E_{on} \dot{\alpha} \right)\end{split}\end{equation}

The total aerodynamic angle-of-attack ( $\alpha$ ) is a summation of the angles due to the pitching ( $\theta$ ) and plunging ( $\dot{h}/b$ ) motions. $C_{L\ell}$ and $C_{M\ell}$ represent the static lift and moment coefficients in the linear region, with $\Delta C_L$ and $\Delta C_M$ the difference between the extended linear curve and the true static curve. The operators $\dot{()}$ and $\ddot{()}$ represent derivatives with respect to the non-dimensional time, $\tau = U t / b$ . This leaves the empirical parameters $\lambda_{on}$ , $s_{on}$ , $a_{on}$ , $\sigma_{on}$ , $r_{on}$ and $E_{on}$ which are derived from wind tunnel data.

In addition to the pitching aerofoil study of Reddy & Kaza [44], this model has been used to obtain the dynamic response of a typical helicopter blade section [53] and the entire helicopter blade [54]. For the dynamic stall study [44], all models were used across a range of test cases for the NACA 0012 aerofoil. For the lift coefficient, all models correlate well across the range of test conditions to the experimental data. The moment coefficient was not presented for the NACA 0012 test cases using the ONERA model, however, it was seen to predict the moment coefficient well for a reference OA212 aerofoil section. The moment coefficient was not predicted well for the Gormont model, with very good agreement for the Gangwani. The results from the Gangwani model were expected due to the amount of empirical data taken from the experimental results.

Reddy & Kaza [44] concluded from this initial dynamic stall investigation that the Gormont and ONERA models can be used for propeller flows in typical operating conditions due to the blade being primarily within the light stall regime. The ONERA model involves fewer experimental parameters than the Gangwani, with the Gormont only requiring one. This could potentially be a factor in the choosing of an aerodynamic model as it depends on the availability of static aerofoil data of the blade sections. This situation was found in the propeller stall flutter study [43] and, hence, the Gormont model was used due to the lack of blade sectional data. This model was used to provide the stall effects with quasi-steady, incompressible two-dimensional strip theory used to calculate the generalised aerodynamic forces across a range of sections along the blade.

The structural modelling used by Reddy & Kaza [43] involved the derivation of the static-displacement and normal modes and frequencies from a non-linear finite element analysis which was then supplied to the standard modal equations of motion for a propeller blade [55]. The final equations were derived as:

(12) \begin{equation}\ddot{q_i} + 2 \zeta_i \omega \dot{q_i} + \omega_i^2 q_i = \frac{Q_i}{m_i} \qquad i=1,2\ldots n,\end{equation}

where $m_i$ was the generalised mass, $\omega_i$ the natural frequency, $\zeta_i$ the critical damping, $Q_i$ the generalised force terms and $q_i$ the normal coordinate for the $i^{th}$ mode where n modes were supplied. These equations were solved within the time-domain using the Wilson- $\theta$ method [56]. This method is an implicit integration method that assumes a linear variation in acceleration between time-steps. The static and modal displacements were derived using NASTRAN, with four normal modes supplied. The aeroelastic flowchart used by Reddy & Kaza is presented in Fig. 19.

Figure 19. Aeroelastic flow-chart for the numerical study of the SR-2 propeller by Reddy & Kaza [43].

The experimental results of Smith [38] presented the stall flutter boundaries in terms of stress across the blade. Using the Gormont dynamic stall model, qualitative correlation was found between the experiment and numerical study, however, the quantitative comparison of the damping plots showed a conservative estimation of the boundary using this model. This conservative estimation may be acceptable from a design point of view if the boundary is outside the normal operating range of the propeller.

Due to the fact that only a qualitative comparison could be made to the experiments, a numerical comparison of the modal amplitude results was conducted against the study of Bielawa [57]. The investigation of Bielawa used the Gangwani dynamic stall model with a non-linear beam model for the structural modelling. The use of a non-linear beam model for the unswept SR-2 is acceptable due to the design of the blade. However, for a complex structure, such as the swept SR-5, the general finite element analysis of Reddy & Kaza [43] is more suitable. The modal amplitude predictions of both the simulations agreed well. A slight improvement on the damping-boundary plots was found in the Bielawa study whilst using the Gangwani model with the predictions shown in Fig. 20. This slight improvement is shown via the closer approximation of the zero damping pitch angle to the single experimental point.

Figure 20. Damping stall flutter boundary plots for the SR-2 propeller for the experimental [38] and numerical studies [43, 57].

In 2007, Delamore-Sutcliffe [58] conducted a numerical investigation into propeller stall flutter with a comparison made to the experiments of Baker [9]. The investigation involved the development of an unsteady aerodynamic model, coupled with the Brookes and Houbolt equations of motion for a rotating beam, to determine the stall flutter boundaries. The main objective of this work was to transition from the empirical models still used by industry to a more advanced formulation without penalising the fast prediction characteristics of such methods.

Delamore-Sutcliffe built-up the aerodynamic model from the linear to non-linear regime in two-dimensions, with experimental tests conducted to ensure a validated model was derived. Three fundamental models were tested for the attached flow (linear) regime, with this involving two established models and a third derived in the time-domain. The experiments used to validate the attached flow model were those of Piziali in 1994 [59]. Details such as the pressure distribution over a semi-span wing and oscillating aerofoil were presented, thus ensuring an extensive range of data from which to compare.

The first model was defined as the indicial response model and considers the instantaneous changes between two steady states. A review of this method was conducted by Beddoes [60]. With an indicial response function defined for the circulatory and non-circulatory terms as a function of semi-empirical coefficients, the total force coefficient is seen as a sum of the subsequent parts. These semi-empirical coefficients have been discussed by Leishman in which a method for obtaining [61] and validating [62] was presented. Using a finite difference approximation to Duhamel’s integral [63], it was possible to determine the response to an arbitrary time-history. A full derivation was presented by Leishman [64].

The second model used was the Theodorsen model [65, 66] where a pitching and plunging aerofoil is modelled in static conditions. Taking the fixed-wing expressions for lift (L) and moment (M) on an aerofoil with unit span:

(13) \begin{equation}\begin{split}L &= \pi \rho b^2 \left[ \ddot{h} + U \dot{\theta} - b a \ddot{\theta} \right] + 2 \pi \rho U b C(k) \left[ \dot{h} + U \theta + b (0.5 - a) \dot{\theta} \right], \\[4pt]M &= \pi \rho b^3 \left[ a \ddot{h} - U (0.5-a) \dot{\theta} - b \left(0.125 + a^2\right) \ddot{\theta} \right] + 2 \pi \rho U b^2 (0.5+a) C(k) \left[ \dot{h} + U \theta + b(0.5-a) \dot{\theta} \right],\end{split}\end{equation}

where C(k) is the Theodorsen function. This function defines the unsteady aerodynamic component. It is complex and typically defined as:

(14) \begin{equation}C(k) = F(k) + i G(k) = \frac{H_1^{(2)}(k)}{H_1^{(2)}(k) + i H_0^{(2)}(k)},\end{equation}

where $H_n^{(2)}$ are Hankel functions of the second kind, which are Bessel functions of the first and second kinds. With assumptions made for the rotorcraft equivalent, taking the differential across the blade radius allows the equivalent lift and moment equations are to be derived (Equation (15)). Although such a theory would appear restrictive to the Theodorsen function, by substituting different lift-deficiency functions the unsteady airloads for other assumed conditions can be expressed.

(15) \begin{equation}\begin{split}\frac{dL}{dr} &= - \frac{\pi}{4} \rho c^2 \Omega^2 \left[ \frac{\ddot{z}}{\Omega^2} - \frac{r \dot{\theta}}{\Omega} + \frac{(x_a - c/4)\ddot{\theta}}{\Omega^2} \right] + \frac{1}{2} \rho c_{l_{\alpha}} c r \Omega^2 C(k) \left[ -\frac{\dot{z}}{\Omega} + r \theta + \frac{\left(c/2 - x_a\right)\dot{\theta}}{\Omega} \right] \\[4pt]\frac{dM}{dr} &= - \frac{\pi}{4} \rho c^2 \Omega^2 \left[ \frac{\left(x_a - c/4\right)\ddot{z}}{\Omega^2} + \frac{r\left(c/2 - x_a\right)\dot{\theta}}{\Omega} + \left(\frac{3}{32} c^2 - \frac{c}{4} x_a + x_a^2\right) \ddot{\theta} / \Omega^2 \right] \\[4pt] & \quad + \frac{1}{2} \rho c_{l_{\alpha}} c x_a r \Omega^2 C(k) \left[ -\dot{z}/\Omega + r \theta + \frac{\left(c/2 - x_a\right)\dot{\theta}}{\Omega} \right]\end{split}\end{equation}

Due to the fact that the Theodorsen model was constructed within the frequency domain, an ideal coupling to a time-marching aeroelastic model remains tricky. This leads on to the third model, which was a time-domain approximation of the Theodorsen model derived by Dowell [67]. Here, it was found that a link between the frequency and time-domain can be obtained via a Fourier transform pair. This results in an updated Theodorsen function with coefficients $\lambda_1$ and $\lambda_2$ , which can be substituted into the full lift and moment equations (Equation (13)) before converting to the final time-domain differentials by using the inverse of the differential rule for Laplace transforms.

(16) \begin{equation}C(k) = \frac{\lambda_2 i k + \lambda_1}{i k + \lambda_1}\end{equation}

All attached flow models were tested and provided satisfactory correlation with the experimental results [59]. The third model was chosen as the basis for the full aerodynamic model due to the lack of empirical data and formulation within the time-domain.

Following the derivation of the linear aerodynamic model, an extension was required to include the dynamic stall effects. Using a state-space approach, a non-dimensional parameter was introduced in order to provide the detached flow mechanism. This parameter was the static upper surface separation point ( $x_s$ ) and this was combined with the Kirchoff and Helmholtz theory [68] to determine equations for the normal force and moment coefficients. The separation point was defined as $x_s \in [0,1]$ , where $x_s = 0$ corresponds to detached flow at the aerofoil leading edge and $x_s = 1$ to a fully attached flow. The normal force coefficient ( $C_N$ ) was given by:

(17) \begin{equation}C_N = c_{l_{\alpha}} \left( \frac{1+\sqrt{x_s}}{2} \right)^2 \alpha,\end{equation}

where $c_{l_{\alpha}}$ is the lift curve slope and $\alpha$ the aerofoil angle-of-attack. For the moment coefficient, a reference to the work of Khrabrov and Ol was made as they produced a theory which combined thin aerofoil theory and Kirchoff concepts to determine a closed solution for the pitching moment coefficient at high angles of attack [69]. However, to overcome the small angle assumption, the work of Leishman and Beddoes is cited [50] for the ratio equation between the normal force and moment coefficients ( $C_M$ ), due to the fact it used the separation point. The ratio was given as:

(18) \begin{equation}\frac{C_M}{C_N} = k_0 + k_1 \left(1-x_s\right) + k_2 \sin \left(\pi x_s^2\right),\end{equation}

where $k_0 = (0.25-x_{ac})$ is the non-dimensional aerodynamic centre offset from the quarter chord. The remaining two parameters, $k_1$ and $k_2$ , are semi-empirical coefficients derived from experiments. In order to estimate the separation point, further semi-empirical coefficients were required depending on the current state of the flow and the static stall angle. Leishman and Beddoes stated that the static stall angle ( $\alpha_s$ ) corresponds closely to the point where $x_s=0.7$ . This fact was used alongside experimental data to determine the static stall angle. The separation point was then given by:

(19) \begin{equation}x_s = \begin{cases} 1 - 0.3 e^{\left[ s_1 (\alpha - \alpha_s) \right]} & \alpha \leq \alpha_s \\ 0.04 + 0.66 e^{\left[ s_2 (\alpha_s - \alpha) \right]} & \alpha > \alpha_s \end{cases}\end{equation}

where $s_1$ and $s_2$ were semi-empirical coefficients chosen to define the shape of the state variable curve pre- and post-stall. In order to model the unsteady effects of detached flow, an effective angle-of-attack was introduced:

(20) \begin{equation}\alpha_E = \alpha - \tau_2 \dot{\alpha},\end{equation}

where $\tau_2$ defines the time delay associated with quasi-steady effects. The model was completed by modelling the transient aerodynamics associated with the relaxation process following detached flow. These effects are determined via a first order differential equation for the effective separation point $x_{s_E}$ :

(21) \begin{equation}\tau_1 \frac{d x_{s_E}}{d t} + x_{s_E} = x_s,\end{equation}

where $\tau_1$ was the semi-empirical relaxation time constant.

This dynamic stall model was validated against the experiments of Jumper et al. [70]. Jumper et al. conducted wind tunnel tests of the pitching NACA 0015 aerofoil over a range of constant pitch rates. The aerodynamic model was tested over a range of the experimental rates with the semi-empirical parameters $s_1$ , $s_2$ and $\alpha_s$ determined from the data. Initial estimates of the two dynamic coefficients ( $\tau_1$ and $\tau_2$ ) were sourced from the studies of Goman and Khrabrov [71] and Khrabrov and Ol [69]. These estimates were then adjusted to provide improved correlation. Excellent correlation was found across all pitch rates.

In order to combine the attached and detached flow effects, a similar summation as the linear model was used. This involved the addition of the zero-angle, circulatory and non-circulatory lift coefficients, with the high angle-of-attack effects included via an unsteady factor (LF).

(22) \begin{equation}C_L = C_{L_0} + ( C_{L_c} \times LF ) + C_{L_{nc}}\end{equation}

This unsteady factor was multiplied by the circulatory lift and was derived in the Kirchoff/Helmholtz theory [68] as:

(23) \begin{equation}LF = \left( \frac{1 + \sqrt{x_{s_E}}}{2} \right)^2\end{equation}

A similar process was found for the pitching moment term. An additional unsteady pitching moment factor (PF) was introduced and this was multiplied by the unsteady lift factor and circulatory lift. This term was then included in the summation of the zero-angle, circulatory and non-circulatory moment coefficients.

(24) \begin{equation}\begin{split}C_M &= C_{M_0} + \left( C_{L_c} \times LF \times PF \right) + C_{M_c} + C_{M_{nc}} \\PF &= k_0 + k_1 \left( 1 - x_{s_E} \right) + k_2 \sin \left(\pi x_{s_E}^2\right)\end{split}\end{equation}

The full aerodynamic model was verified against the pitching aerofoil experiments of Piziali [59], with pre-requisite testing used to determine the semi-empirical coefficients. Good correlation was found in static and dynamic tests, with some small deficiencies at high static angles of attack and the dynamic overshoot of the lift and pitching moments.

In order to verify the aerodynamic model for the stall flutter phenomena, Delamore-Sutcliffe conducted a range of wind tunnel experiments of a two-dimensional aerofoil with a single degree of freedom in pitch. This verification was conducted via a time-marching and eigenvalue stability analysis through the coupling of the aerodynamic and structural models. The structural model was derived using Lagrange’s equations of motion [72]. Due to the motion of the aerofoil, the potential remains for the aerofoil to stall at a negative angle-of-attack. As a result, the equations for the separation point needed updating to include the negative range with respect to the negative static-stall angle ( $\alpha_{s_n}$ ):

(25) \begin{equation}x_s = \begin{cases} 1 - 0.3 exp \left[ s_{1_n} \left(\alpha_{s_n} - \alpha\right) \right] & 0 > \alpha > \alpha_{s_n} \\[6pt] 0.04 + 0.66 exp \left[ s_{2_n} \left(\alpha - \alpha_{s_n}\right) \right] & \alpha \leq \alpha_{s_n} \end{cases}\end{equation}

where semi-empirical parameters $s_{1_n}$ and $s_{2_n}$ were updated for the negative range. The equations were coupled in a second order matrix equation including the mass, damping and stiffness matrices. A Runge-Kutta scheme was used for the time-marching formulation and a linearised solution in the form of an exponential was assumed in order to obtain eigenvalues for the stability analysis. Further details were presented by Delamore-Sutcliffe et al. [73]. A comparison of the experimental and numerical results showed good correlation, with the numerical results from both the stability and time-marching methods slightly below the experiment.

In preparation for the transition to the propeller test cases, the aeroelastic model developed for the two-dimensional aerofoil was extended to look at a three-dimensional cantilever wing. The wing was structurally modelled using strip theory and assumed modes with two degrees of freedom in pitch and heave. A key benefit to using this method was the limitation of two-dimensional aerofoil characteristics and the computational solution time. Disadvantages included the lack of three-dimensional crossflow modelling and the requirement for empirical data, a common theme in rotor aerodynamic modelling.

As per the structural model, the aerodynamic model was updated to compute the spanwise lift and moment coefficients in both attached and detached conditions. With the use of strip theory, the aerodynamic model was applied at each spanwise location with a station-by-station update required for the semi-empirical coefficients and lift curve slope. Experiments [59] were again used for the coefficients, along with curve fitting, to determine the spanwise distribution of the values. Theoretical expressions used to convert the two-dimensional slope to three-dimensions with a quintic power distribution over the span was assumed for the lift curve slope with a value of zero applied at the tip. Furthermore, the effect of cantilever wing twist was taken into account via lifting-line theory and Prandtl’s integral equation for circulation [74].

These expressions for aerodynamics and structure were combined in a similar manner to the two-dimensional aerofoil equations. The structural equations were generalised with respect to the assumed modes and then integrated over the span. Two bending and a single torsional mode were included, with higher modes having a minimal effect on the accuracy. A similar process was conducted for the aerodynamic terms including the circulatory and non-circulatory coefficients. The solution to these equations were obtained via a time-marching and eigenvalue analysis. Due to the increase in matrix terms, the time-marching calculations were limited due to the increase in computational cost. These were subsequently compared to the eigenvalue analysis which was conducted across the entire test case range.

This three-dimensional wing method was compared against the experiment of Dunn and Dugundji [19]. Good correlation between the experimental and numerical calculations were found. For the first test case examined, the time-marching method appeared to produce an improved correlation of the flutter frequency in comparison to the stability analysis. However, the flutter velocity was under-predicted at high initial pitch angles. Only the stability analysis was used for the second test case with an under-prediction of the flutter velocity at low and high initial angles, and under-prediction of the frequency across the full range.

Following the verification of the three-dimensional wing aeroelastic model, the aerodynamic model was extended and used with the Brookes and Houbolt equations of motion for a rotating beam [75] to conduct a propeller stall flutter investigation. For the structural model derived by Brookes and Houbolt, several assumptions were made to reduce the complexity of the equations without losing accuracy. These were:

1. (1) The chordwise bending was assumed to be very small, and hence negligible, in comparison to the flapwise bending. This was due to the fact that the chordwise bending stiffness is much greater than the flapwise stiffness when the thickness/chord ratio is less than 15%. This is acceptable with modern propellers due to their need to utilise thin aerofoil sections for performance. This assumption may fail for tiltrotor blades where thicker sections are used.

2. (2) The blade was assumed to have no set twist distribution in the spanwise direction. This was applicable within this study due to the fact that the experimental blades have zero twist. For more complex, highly twisted blades, this assumption would fail and the terms associated would need to be included.

3. (3) The final assumption was that the tensile and rotational axis of the blade align with the elastic axis. This was, again, driven by the experimental data were no information was supplied regarding these offsets from the elastic/mass axis.

These assumptions resulted in the following structural equations for the propeller blade:

(26) \begin{equation}\begin{split}L &= \frac{\partial^2}{\partial y^2} \left(E I_1 \frac{\partial^2 w}{\partial y^2}\right) - \frac{\partial}{\partial y} \left(Ta \frac{\partial w}{\partial y}\right) - \frac{\partial}{\partial y} \left( \Omega^2 m (y+R_{hub}) e \alpha \right) + m \left( \frac{\partial^2 w}{\partial t^2} + e \frac{\partial^2 \alpha}{\partial t^2} \right) \\[4pt]M &= - \frac{\partial}{\partial y} \left( \left(GJ + Ta k_y^2\right) \frac{\partial \alpha}{\partial y} \right) + \Omega^2 m \left(y+R_{hub}\right) e \frac{\partial w}{\partial y} \\[4pt] &\quad + \Omega^2 m \left(k_{zz}^2 - k_{xx}^2\right) \alpha + I_y \frac{\partial^2 \alpha}{\partial t^2} + m e \frac{\partial^2 w}{\partial y^2}\end{split}\end{equation}

Where the axial tension was defined as $Ta = \int_y^L \Omega^2 m (y + R_{hub}) dy$ and e the distance between the elastic and mass axes was given by $e = (x_{ea}-x_{cg})c$ . An important point to note is the reduction of the above equations from propeller to wing configurations is possible via the setting of the rotational velocity ( $\Omega$ ) and hub radius ( $R_{hub}$ ) to zero. This provides continuity in the aeroelastic modelling.

In terms of the aerodynamic modelling, the strip theory used in the cantilever wing was updated to include the spanwise variation in induced velocity. One of the key physics of rotorcraft aerodynamic modelling is the influence of the wake on the preceding blade. This effect can be captured by methods such as Loewy [76] who took into account this effect via the modification to the lift-deficiency function from Theodorsen’s theory. Although such a method aims to capture a greater amount of the physics, it was assumed that such an effect is minimal in comparison to the overall load results obtained via Theodorsen’s method. This is true in high thrust/inflow condition were the propeller wake is induced far away from the disc.

In order to update the aerodynamic model for the rotorcraft environment, the effect of the induced incidence must be taken into account. The induced incidence is the effect of the blade rotation and inflow velocity, and was modelled via Blade Element Momentum Theory (BEMT) [77]. This model was used for each section with an iterative process conducted to obtain the induced velocity effect.

Both the structural and aerodynamic terms were combined as the previous versions, in state-space form of stiffness, damping and mass matrices. In the cantilever wing model, the structural damping was set to zero. This was not the case for the propeller model due to the fact it was taken into account in the experimental studies [9] and, hence, must be included within the equations. To include structural damping, the theory of Theodorsen [66] was applied where the structural damping was seen as a function of the generalised stiffness matrix, i.e.:

(27) \begin{equation}[C] = i g [K],\end{equation}

where, g is the structural damping coefficient. In order to solve the equations of motion, the stability analysis conducted for the aerofoil and cantilever wing was used. Due to the increased computational cost associated with the time-marching method, it was subsequently not used for this method. Further details regarding the aeroelastic modelling was presented by Delamore-Sutcliffe and Greenwell [78].

A comparison to the experimental data of Baker [9] was made with three test cases selected and presented in Fig. 21. A good correlation between the experiment and numerical analysis was found. The trend of lowering velocity with increased pitch angle was seen in both results with the minimum flutter velocity well captured. Differences at higher angles of attack were seen within the comparison and these were determined to be an effect of the high angle-of-attack detached flow modelling. However, the pitch angles where this was seen are outside the blades operational regime.

Figure 21. Stall flutter boundary comparison for the model propellers of Baker [9] compared against the numerical aeroelastic model of Delamore-Sutcliffe [58].

Ognev in 2011 investigated different unsteady aerodynamic models and their influence on the determined flutter boundary [79]. These unsteady aerodynamic models included Theodorsen’s model [80], three cascade models [81, 82, 83] and a detailed three-dimensional model [84]. The cascade models took into account the effect of the returning wake (something which is neglected in the Theodorsen model). However, they were limited by their two-dimensionality with the effect only applied from equivalent cross sections. This was mitigated by the use of the three-dimensional model by Isoilevskii et al. [84] where bound, trailing and shed vortices were included. For all models, the effect of viscosity is neglected and this resulted in the loss of detached flow modelling, not ideal for a stall flutter analysis, hence the comparison between the simulations and experiments [33] were only fair within the classical bending-torsion flutter regime.

In 2019, Higgins et al. conducted a numerical investigation into the Aero Twin Commander propeller blade [85]. This investigation used the HMB3 CFD solver and compared the results of URANS and SAS [21] calculations against experimental data. The use of the SAS method for aerodynamic modelling had proven beneficial in a typical two-dimensional form via the study of pitching aerofoils [20]. Such a method allows for a balance between the resolved and statistical content to be found, and improves the accuracy without significantly increasing the computational cost.

For the structural modelling, a loosely coupled time-marching modal method was used. This method, again, allows for a balance to be made between modelling accuracy and the heavy computational cost associated with coupled methods. For the modal method, the mode shapes and frequencies of the Commander propeller blade were required, and for this, a NASTRAN non-linear beam model was built. A non-linear beam model was selected due to its ability in previous helicopter rotor aeroelastic studies [86], with good validation found between the experimental and numerical mode shapes and frequencies for the propeller blade.

As in the experimental investigation [37], the blade was accelerated from 1,400 to 1,750rpm at a fixed pitch angle of $27^{\circ} 40^{\prime}$ . Good correlation was found in terms of the trend of torsional stress across the blade radius between test and simulations (Fig. 22(b)). The SAS results were found to be slightly higher than the URANS, indicating the need for scale-resolving. This was confirmed with the modal amplitude results where a qualitative agreement to the experiment was found in terms of the third mode amplitudes. As the blade accelerates, a significant increase in amplitude was found which mirrors the trend in stress from the test results (Fig. 22(a)).

Figure 22. Time history of the modal amplitude results and torsional stress trends across the blade for the Commander investigation [85].

This investigation highlighted the ability and requirement for scale-resolving CFD methods when studying stall flutter. Advances in computing power has enabled the use of such methods, thus allowing for closer validation which can lead to greater advances to be made in design. Expanding the boundary which the blade can operate.

Presented in Table 5 is a summary of the numerical investigations conducted for propeller stall flutter. A list of the techniques used, the range of conditions and the date is presented. From the table, and this section, it can be found that a limited number of aeroelastic numerical investigations have been conducted. This is likely down to the small range of potential trigger points during a typical flight path. Thus, allowing for fast, eigenvalue type studies to be dominant within the field as the ultimate aim of avoidance could be achieved from conservative boundaries. However, again with the introduction of the eVTOL market, propeller stall flutter is becoming a more active area of research due to shift in propeller flight path and the re-design of the blade for greater efficiency.

Table 5. Summary of propeller stall flutter numerical test cases.

4.0 Concluding remarks

4.1 Discussion

A range of experimental investigations have been carried out over the years. Due to the heavy use of propeller driven aircraft between 1945 and 1960, a high number of experiments were conducted to understand blade flutter. As these studies were not specifically focused on stall flutter, merely propeller flutter in general, correlation with more recent studies is difficult as only a few fundamental parameters were examined. This included rotational velocity, pitch and number of blades. However, these studies by Sterne [29] and Theodorsen & Regier [30] provided some of the initial trends and influences stall flutter can have on propellers, with this being the significantly reduced boundary as a result of a propeller entering stall flutter.

One of the additional factors studied by Sterne [29] was the effect of number of blades. A different number of blades was studied by Theodorsen & Regier, but due to the fact that one of the two propellers was single-bladed, a trend with blades numbers would be invalid due to lack of influence on the single blade from a proceeding blade wake. However, the study by Ewing et al. [32] did examine different blade numbers. Ewing et al. found that the effect of number of blades had an influence on the vibratory phase between blades and, as a result, any comparison in terms of stress levels between propellers of different blade numbers would also have to take into account the phase. Such a comparison remains difficult due to the fact that a peak stress level for one propeller may aline with another which induces a phase change and thus the stress levels are lower. This is potentially one of the influences that was found in the study by Sterne in that a different phasing was introduced by the 3-bladed propeller that resulted in higher stress levels, and thus lower stall flutter boundaries, than the 2- and 4-bladed alternatives.

In the 1950s, a range of model blades were designed for wind tunnel investigations [9, 32, 33, 34, 36]. Such blades were made of a modern alloy to allow for a reduction in the variability of the manufacturing process. A range of parameters could be investigated and the greatest insight into propeller stall flutter was found during this period. One of the parameter investigated during this period was the effect of blade twist, with this parameter being the only factor examined across the different studies. This was primarily conducted by Baker [9] and Hubbard et al. [33]. However, during the tests, only two data points were examined. In the case of Baker [9], an untwisted and $17^{\circ}$ twisted blade were compared with very little difference found in terms of minimum flutter speed. For the investigation by Hubbard et al. [33], again an untwisted and a $31^{\circ}$ twisted blade were tested but the untwisted blades were conducted at different conditions. Hence, a true side-by-side comparison could not be conducted.

With the transition from propeller-driven aircraft to turbofan engines, interest was lost in propeller stall flutter; therefore, very few studies were conducted after the 1950s. A private-sector test was conducted by Dowty Propellers in 1979 using an in-service blade [37]. Although the primary aim of the test was not specifically for stall flutter understanding, stall flutter was encountered due to the pitch angle range. Furthermore, interest within propfans grew in the 1980s and this resulted in new studies to understand the blade flutter boundaries [38, 40]. These investigations into the propfans by Smith looked into a range of factors previously examined. This included blade sweep previously studied by Baker [9] and freestream velocity and yaw which was previously examined by Rogallo & Yaggy [36].

Focusing on the sweep trends produced by Smith [38, 40] and Baker [9], Baker initially found the positive benefit of sweep. An unswept blade was compared to three swept blades, two of which had a linear profile and the other a non-linear. Although the non-linear profile had the greatest impact across the stall pitch region, both the linear swept blades produced a benefit in terms of the flutter velocity. The comparisons of sweep by Smith looked at the unswept SR-2 against the non-linear swept SR-3 and SR-5. Both the SR-3 and SR-5 reduced the levels of stress acting on the blade, confirming the trend found by Baker. One of the benefits to the introduction of sweep is the delay in separation onset with minimal performance losses. This correlates with the stall flutter results as increasing sweep is seen to delay the onset of high stress.

The majority of the experimental studies have been conducted in static/hover conditions. This is mainly due to the fact that propeller stall flutter is likely to occur during low inflow conditions, such as take-off. However, the experiments of Rogallo & Yaggy [36] and Smith [40] looked into the effect of added freestream velocity and found that the overall stress acting on the blade was significantly reduced and thus an improvement in the stall flutter boundary was achieved. The improvement being the requirement for higher rotational velocity and pitch angles before stall flutter is induced compared to the hover trigger values. Such a response is expected when the effective angle-of-attack is taken into account. Although the geometric pitch remains the same, the additional inflow vector reduces the overall effective angle and thus, stall is delayed. Such a trend is not only observed in propellers but also wings and aerofoils.

The final cross-over factor examined experimentally was the effect of yaw and this was again studied by Rogallo & Yaggy [36] and Smith [40]. In order to examine for the effect of yaw, a freestream velocity is required. For the Rogallo & Yaggy test, only the highest freestream condition was reported and therefore, any yaw influence on the stall flutter boundary is minimised due to the dominating reduction in boundary of the freestream. This is essentially confirmed due to the fact that almost no change in stress was found as a result of yaw. Smith did present results for the yawed cases across a range of freestream velocities, however, redefined the stress values presented by taking into a yaw factor. As a result, no direct comparison to a non-yawed case could be made with only observations presented by Smith available. At the lowest freestream velocities, Smith presented a linear trend of total maximum stress with yaw angle and this is induced due to the one-per-revolution fluctuation in loads due to the variation of velocity across the blade disc. This additional fluctuation could have a significant impact on the stall flutter boundary, however, the dominating effect of the freestream must be minimised in order to gain a true picture of its effects.

With the minimal development in aeroelastic models during the 1950s, interest in numerical stall flutter investigations only began with the development/study of the SR propfans. These studies, conducted in the 1980s, involved both eigenvalue and time-marching modal methods with semi-empirical aerodynamic models. The use of the eigenvalue analysis was selected to ensure a fast calculation and this trend continued into the modern studies. With the availability of faster computing power, more complex finite element models could be used and a greater amount of time-marching analysis was conducted. However, the focus still remained on the use of the semi-empirical dynamic stall models and therefore conservative boundaries were achieved.

4.2 Conclusion

From the initial studies of Sterne [29] and Theodorsen & Regier [30], wooden propellers are used. One of the benefits of using modern alloys is the ability to refine the manufacturing process to reduce any inconsistencies in structural properties. Modal analysis can indicate frequency differences in blades of the same design. Therefore, to ensure no additional inconsistencies are applied on top of numerical errors, metal alloy propellers are preferable, making wooden blades non-ideal.

It is clear from the studies in the 1950s [9, 32, 33, 34, 36] that a great deal of data and insight was gleaned. However, due to the purpose-built nature for each of the studies, the blades were not fully defined within the reports. As a result, such investigations remain non-ideal for numerical investigation. Additionally, no significant cross-over of parameters was conducted and therefore any trends, regardless of propeller type, remain sought after.

For each experimental investigation described, the structural properties are unknown and unreported. These properties, therefore, have to be estimated in order to conduct numerical investigations.

From the propeller stall flutter review it is concluded that there is a lack of definitive propeller stall flutter data. The experimental investigations are historic and often lacking in clear boundary assessment. This minimises the range of validation that can be carried out with modern numerical methods.

4.3 Recommendations

For any future study, the influence of the freestream velocity must remain minimal with ideal test conditions involve zero to low speed conditions. The exact boundary of the low speed condition is dependent upon the specific design. However, from the studies of Rogallo & Yaggy [36] and Smith [40], even the lowest freestream increment had an effect during low rotational velocity conditions. In the Rogallo & Yaggy study, the lowest ratio of freestream to tip speed was $\sim5\%$ and this $\sim5\%$ introduction of freestream found an increase in pitch of $5^{\circ}$ , or 22%, with respect to the hover boundary. This highlights the significant influence freestream velocity has and hence the recommendation.

In order to improve the numerical techniques, and additionally obtain a greater insight into the phenomena via non-intrusive tools, a new set of experiments is required to extract greater datasets. Not only with the current datasets being minimal and historic, it is further emphasised that modern blades for novel aircraft configurations are likely to be considerably different from wooden propellers of 80 years ago, and operate in completely different regimes (from eVTOL to tilt-rotor aircraft), and are made of completely new materials, such as composites, which have different aeroelastic response. For these new systems, there are virtually no data.

A range of numerical methods are available from which to study propeller stall flutter, with each model’s application dependent upon the overall goals. Should a significant range of calculations be required, an eigenvalue analysis with semi-empirical aerodynamic model would be of use as this would allow for a full conservative boundary to be formed. Should greater detail about a specific condition be required, a time-marching method with high-fidelity CFD aerodynamics could be used as this would allow for an assessment of the structural and aerodynamic response in time. With the latter more computationally expensive, improvements in computing power are continually being made and therefore the likelihood of using a time-marching method is increasing.