2.1. Stall flutter under uniform freestream conditions

Stall flutter is defined as a self-sustained large-amplitude oscillation in which the flow supplies net energy to the elastically mounted aerofoil through dynamic-stall effects. Classical observations of this aeroelastic instability reveal free pitch oscillations about a mean angle of attack close to or within the static stall regime (Ham & Young Reference Ham and Young1966; McCroskey Reference McCroskey1982). The stall flutter investigated in the current study is peculiar in the sense that the self-sustained pitch oscillations occur symmetrically about a zero angle of attack. It has been analysed in some detail both experimentally in one-degree-of-freedom pitch-only and two-degree-of-freedom pitch–heave (Poirel et al. Reference Poirel, Goyaniuk and Benaïssa2018; Goyaniuk et al. Reference Goyaniuk, Poirel and Benaïssa2020), and numerically in one degree of freedom pitch-only (Goyaniuk et al. Reference Goyaniuk, Poirel, Benaïssa and Yuan2017; Zhang et al. Reference Zhang, Poirel and Yuan2021b ). Both two-dimensional unsteady Reynolds-averaged Navier–Stokes simulations and three-dimensional aeroelastic large-eddy simulations (LES) were performed. Although both sets of simulations were able to reproduce qualitatively the experimental observations, as expected, the LES-based results captured reasonably well the oscillation amplitude and frequency.

Similar experimental observations were reported by Dimitriadis & Li (Reference Dimitriadis and Li2009) for a NACA0012 with an EA position at $0.38c$ , and by Zhu et al. (Reference Zhu, Su and Breuer2020) also on a NACA0012 but with the EA at the mid-chord. In both cases, the stall flutter occurred via a subcritical bifurcation with its onset airspeed past the divergence airspeed.

Figure 3 summarises the experimental response of the pitch-only flutter rig, treated in the current work, in a steady uniform freestream. The stall-flutter onset speed is $U_{\textit{SF,onset}}=6\,\mathrm{m\,s^{-1}}$ ( $ \textit{Re}\approx 6.2\times 10^{4}$ ), and the static loss of equilibrium at approximately $7.5\,\mathrm{m\,s^{-1}}$ is due to divergence. Thin-aerofoil theory predicts divergence to occur at $U_D = 7.32\,\mathrm{m\,s^{-1}}$ ; a slightly lower divergence airspeed is obtained if the measured value of the moment at the zero angle of attack is used for the calculations instead of the theoretical value. For $U_{\infty }\lt U_{\textit{SF,onset}}$ , small perturbations decay to the zero-equilibrium state. For $U_{\textit{SF,onset}}\lt U_{\infty }\lt U_{D}$ , two stable attractors coexist (two solid branches), i.e. sufficiently large initial disturbances lead to stall flutter exhibiting large-amplitude oscillation, whereas small perturbations still return to zero equilibrium. When $U_{\infty }\gtrsim U_{D}$ , stall flutter is self-triggered without external perturbations, and the baseline amplitude $\bar {\theta }_0$ approaches a plateau $\sim 41^{\circ }$ near $U_{\infty }=10\,\mathrm{m\,s^{-1}}$ ( $ \textit{Re}\approx 10.3\times 10^{4}$ ). This amplitude plateau of stall flutter was also observed by Dimitriadis & Li (Reference Dimitriadis and Li2009) at a similar level. In addition, previous experiments have shown similar stall-flutter dynamics for both the clean and tripped wings (Goyaniuk et al. Reference Goyaniuk, Poirel and Benaïssa2020).

Regarding figure 3(b), the frequency in quiescent air is the structural natural frequency. As the airspeed is increased to the stall-flutter speed and beyond, its frequency in the no-gust case, $f_{\theta ,0}$ , increases approximately linearly. Taking the point of view of an effective aeroelastic (total) stiffness, which is the addition of the structural and effective aerodynamic stiffness, one can conclude that this effective aerodynamic stiffness grows with dynamic pressure. This interpretation is compatible with the measured nonlinear dependence of the aerodynamic moment on angle of attack to be discussed later in this section. Hence, the concept of effective aerodynamic stiffness is akin to an averaged value of the nonlinear stiffness as it varies within the range of angles encountered during the LCO.

Figure 3.

Baseline stall-flutter (SF) characteristics under uniform freestream (no-gust) conditions. (a) Bifurcation diagram of the LCO amplitude ( $\bar {\theta }_0$ ) against $U_{\infty }$ , where solid blue curves and markers denote the experimental measurements, while the solid black traces provide schematic representations of the attractors observed within each airspeed subrange. (b) Flutter frequency ( $f_{\theta ,0}$ ) trend against $U_{\infty }$ .

Table 1.

Nominal values of structural parameters.

Table 1. Long description

A table comparing structural parameters with their nominal values. The table has one row and six columns. The columns are labeled EA location, I_EA (kg m^2), D_theta (N m s rad^-1), K_theta (N m rad^-1), Chord c (m), and Span, ℓ (m). The row data is as follows: EA location: 0.35c, I_EA: 0.0011 kg m^2, D_theta: 0.0014 N m s rad^-1, K_theta: 0.3 N m rad^-1, Chord c: 0.156 m, Span, ℓ: 0.61 m.

The equation of motion for the pitch degree of freedom is expressed as:

(2.1) \begin{align} I_{EA}\ddot {\theta } + D_{\theta }\dot {\theta } + K_{\theta }\theta = \frac {\rho U^{2}}{2} c^{2}\ell C_{M,\textit{EA}}. \end{align}

The specific structural parameters are given in table 1; note that the structure is assumed to be linear as per (2.1), notably with a viscous damping model and pitch stiffness linearly proportional to pitch angle. Figure 4 confirms this assumption, except at very small angles where dry friction starts to dominate. Note as well that the red dots highlight the local peaks and troughs, which were used to perform a cycle-to-cycle analysis. Both the local frequencies and the damping ratios remained reasonably constant, except at very small angles. Further details can be found in Poirel et al. (Reference Poirel, Goyaniuk and Benaïssa2018) and Goyaniuk et al. (Reference Goyaniuk, Poirel and Benaïssa2020).

Figure 4.

Time history of the pitch free decay response in quiescent air (no flow), overlaid with the exponential decay envelope derived from the parameters in table 1.

A typical stall-flutter time response of the pitch motion is shown in figure 5(a), with the corresponding spectral content shown in figure 5(b). Both plots indicate a nearly simple harmonic motion, with a strong dominant spectral peak defined as the limit-cycle oscillation frequency $f_{\theta ,0}$ . The next strongest peak occurs at the third superharmonic $3f_{\theta ,0}$ , followed by much weaker odd superharmonics. The even harmonics are almost invisible, indicating that the motion is approximately symmetric about the zero pitch angle. The specific airspeed shown is $U_\infty = 10\,\mathrm{m\,s^{-1}}$ , which is the baseline condition used for the subsequent gust-response study.

Figure 5.

Stall-flutter pitch response at $U_\infty =10\,\mathrm{m\,s^{-1}}$ : (a) time history; (b) power spectral density (PSD).

As alluded to previously, the nonlinear dynamics originates exclusively from the aerodynamics. The nonlinear dependence of the aerodynamic moment coefficient on the pitch angle, and hence on the geometric angle of attack, is illustrated in figure 6 for both static and dynamic conditions. The dynamic case follows the pitch-motion schedule over one LCO cycle at $U_\infty =10\,\mathrm{m\,s^{-1}}$ . The dynamic moment coefficient is calculated from the aeroelastic equation of motion (2.1). The measured pitch motion shown in figure 5(a), which is used as the input for the moment calculation, is first filtered to remove all frequency content above $3f_{\theta ,0}$ . A filter cut-off frequency slightly above the third superharmonic provides a good compromise between retaining the essential nonlinear features of the aerodynamic moment and allowing a clear interpretation of the dynamics. Both the filtered pitch response and the resulting aerodynamic moment coefficient are then phase-averaged to preserve the phase dependence between the two quantities.

The arrow in figure 6 indicates the direction of increasing time. The two clockwise loops near the extrema of the LCO correspond to positive work done by the airflow that feeds the stall flutter. Their total enclosed area is larger than that of the anticlockwise loop around the linear regime, indicating a net transfer of energy from the flow to the structure. The change in loop direction occurs at, or slightly above, the static stall angle.

Figure 6.

Phase-averaged aerodynamic moment coefficient about the EA at $0.35c$ under static (dashed red) and dynamic (solid blue) conditions at $U_\infty =10\,\mathrm{m\,s^{-1}}$ .

The static aerodynamic moment coefficient was calculated from direct measurements of the normal and axial forces, as well as the moment about the $0.3c$ position of the static wing using another NACA0012 apparatus with a similar chord length at $U_\infty =10\,\mathrm{m\,s^{-1}}$ . The lift, drag and moment coefficients about any point can then be calculated. For completeness, the lift coefficient and moment coefficient about the quarter-chord point are shown in figure 7 for a Reynolds number compatible with the stall-flutter observations, and for the wing with a boundary-layer strip. Averaged values are shown because some degree of unsteadiness was observed, especially in the high-angle-of-attack regime. Further details can be found in Zhang et al. (Reference Zhang, Poirel and Yuan2025) for both the clean and boundary-layer-tripped wings.

Figure 7.

(a) Averaged lift coefficient and (b) aerodynamic moment coefficient about $c/4$ for the boundary-layer-tripped wing at $Re_c=12\times 10^4$ (where AoA means angle of attack).

The effect of the nonlinear aerodynamic stiffness, as shown in the static moment coefficient, can also be observed in the transient pitch response in cases where the pitch motion is triggered by a small initial condition and develops towards the steady-state LCO. An evolution of the instantaneous oscillation frequency, as measured peak-to-peak, is thus observed with small frequency values at low pitch amplitudes, and larger frequencies at high pitch amplitudes.

Static stall occurs at angle of approximately $10^\circ$ , which is lower than the classical values obtained for the NACA0012 aerofoil at higher Reynolds numbers (Abbott & von Doenhoff Reference Abbott and von Doenhoff1959), indicating that Reynolds number effects remain present even though a boundary-layer strip is used. Reynolds number effects are also observed at lower angles, especially in the moment coefficient, through its nonlinear dependence on angle of attack.

2.2. Gust generation

As shown in figure 2, to excite the flutter rig with controlled gusts, a gust generator is installed upstream of the tripped flutter wing. The gust generator, composed of two NACA0015 vanes, was initially designed for our last-stage study, i.e. gust-induced dynamic stall (Zhang et al. Reference Zhang, Poirel and Yuan2025), with a stationary wing mounted downstream with an internal balance. The two vanes with chord length $c_V = 0.3$ m are controlled by two independent AC motors. A constant-temperature anemometer with an X-hot-wire probe is employed to measure the flow variations in two directions: $u$ and $v$ components. In-phase vane motion generates transverse gusts; anti-phase motion generates longitudinal gusts. The probe was positioned at the mid-span (mid-height) plane of the test section and at the primary streamwise location corresponding to the leading edge of the downstream flutter wing, with the wing removed. This location was selected because it represents the inflow encountered by the flutter wing while avoiding direct interference from the wing motion and wake. This location also ensures the gust amplitude behaves nearly uniformly along the streamwise direction (Yuan et al. Reference Yuan, Zhang and Poirel2024; Zhang et al. Reference Zhang, Poirel, Benaïssa and Yuan2024). With the anti-phase vane motion, the longitudinal gusts generated are depicted in figure 8, behaving in the nearly harmonic pattern as desired, $U(t) = U_\infty +u_{g,\textit{max}}\sin (2 \pi f_g t)$ . The $v$ component variation is negligible compared with the $u$ component. When the vane oscillation is limited to $12^\circ$ (to avoid vane stall), the streamwise gust ratio $\varepsilon _g = {u_{g,\textit{max}}}/{U_\infty }$ depends on factors such as $U_\infty$ and the vane frequency, which is equal to the gust frequency $f_g$ . At $U_\infty = 10\,\mathrm{m\,s^{-1}}$ for the current facility, $\varepsilon _g$ ranges from 0.09 to 0.14, with $f_g$ varying from 2 to 4 Hz.

Figure 8.

Phase-averaged (a) $u$ and (b) $v$ components downstream of the anti-phase oscillating vanes (vane amplitude $12^\circ$ ) with various frequencies at $U_{\infty } = 10\,\mathrm{m\,s^{-1}}$ .

4.1. Aerodynamic stiffness simplification

Before proposing a minimal model, how a longitudinal gust $U(t)=U_{\infty }+u_g(t)$ modulates the system needs to be identified. In the canonical VdP oscillator, the nonlinearity resides in the state-dependent damping (shown later). To capture the parametric action of the streamwise gust, a modulation of the effective pitch stiffness is introduced, and the damping modulation will be discussed later. Classical aeroelasticity predicts that the aerodynamic stiffness felt by a lifting surface in pitch scales quadratically with the freestream speed (Fung Reference Fung1955; Poirel & Price Reference Poirel and Price2003; Zhang et al. Reference Zhang, Kheiri and Xie2020). At the most basic level, aerodynamic pressure scales as $q=(1/2)\rho U^2$ , and the linearised moment response to a small angular displacement contains a stiffness-like term proportional to $q$ (times appropriate geometry/aerodynamic derivatives). Consequently, if a longitudinal (streamwise) gust induces a small relative speed variation, then

(4.1) \begin{align} U(t)=U_\infty \bigl [1+\varepsilon _g \cos (\omega _g t)\bigr ]. \end{align}

Accordingly, the aerodynamic stiffness inherits a parametric dependence through

(4.2) \begin{align} K_a(t) \propto U(t)^2 = U_\infty ^2\bigl [1+\varepsilon _g \cos (\omega _g t)\bigr ]^2. \end{align}

Expanding the square shows what this injects into the pitch equation:

(4.3) \begin{align} \bigl [1+\varepsilon _g \cos (\omega _g t)\bigr ]^2 = 1+2\varepsilon _g\cos (\omega _g t)+\frac {\varepsilon _g^2}{2}\bigl [1+\cos (2\omega _g t)\bigr ]. \end{align}

Three ingredients appear: (i) a small constant stiffening $+( {\varepsilon _g^2}/{2})$ , (ii) a dominant tone at $\omega _g$ of order $\mathcal{O}(\varepsilon _g)$ , and (iii) a weaker tone at $2\omega _g$ of order $\mathcal{O}(\varepsilon _g^2)$ .

Considering small–moderate gust intensities, $\varepsilon _g\ll 1$ , the following practical simplification is adopted.

1. (i) The parametric stiffness with a single-tone modulation is modelled as

   (4.4) \begin{align} \bigl [1+\varepsilon \cos (\omega _g t)\bigr ]\,x, \end{align}
   where $x$ is the normalised pitch angle, and $\varepsilon \approx 2\varepsilon _g$ .

2. (ii) The constant stiffness term $+\varepsilon _g^2/2$ produces only a small upward shift of the carrier frequency $\omega _{0,{eff}}\approx \omega _0(1+\varepsilon _g^2/4)$ , which is higher order and negligible at the gust levels considered ( $\sim 0.1$ ).

3. (iii) The $2\omega _g$ tone has amplitude $\mathcal{O}(\varepsilon _g^2)$ and yields much weaker second-order sidebands.

Therefore, the single-tone parametric form $\bigl [1+\varepsilon \cos (\omega _g t)\bigr ]\,x$ captures the leading-order physics while keeping the derivation compact and transparent.

Stall flutter under longitudinal gusts is modelled by the conventional Mathieu-type van der Pol (M-VdP) equation with stiffness modulation and no external forcing:

(4.5) \begin{align} \fbox{$ \ddot {x}-\mu \omega _0\bigl (1-x^2\bigr )\dot {x}+ \omega _0^2\bigl [1+\varepsilon \cos (\omega _g t)\bigr ]\,x=0, \quad |\varepsilon |\ll 1,\ 0\lt \mu \ll 1 $ }, \end{align}

where $x = {\theta }/{\bar {\theta }_0}$ is the normalised pitch angle, $\omega _0$ is the natural flutter frequency (no-gust), and $\mu$ is the effective aeroelastic damping parameter governing the nonlinearity strength. The term $-\mu \omega _0\bigl (1-x^2\bigr )\dot {x}$ ensures self-excitation with smaller $x$ , and saturation with larger $x$ . A detailed derivation connecting this reduced oscillator to the physical aeroelastic pitch equation is given by Zhang (Reference Zhang2026, ch. 6, § 6.3.1) and is not repeated here. Throughout, $\dot {(\ )}$ denotes time derivatives in dimensional time $t$ .

The present M-VdP model should be interpreted as a minimal phenomenological model for the measured pitch dynamics rather than a general normal-form model for the full fluid–structure system, but its parametric term is motivated by the aerodynamic effect of a longitudinal gust: the streamwise velocity fluctuation modulates the dynamic pressure and therefore the effective aerodynamic stiffness and damping of the aeroelastic system. A broader class of oscillator-flow models is provided by Stuart–Landau-type formulations. For example, Herrmann et al. (Reference Herrmann, Oswald, Semaan and Brunton2020) developed an extended Stuart–Landau model for forced turbulent oscillator flows, and showed that nonlinear mode interactions can produce synchronisation regions, including subharmonic lock-in. The M-VdP model used here has a more restricted purpose: it represents the self-excited stall-flutter LCO through nonlinear damping, and introduces the longitudinal gust as a stiffness modulation, thereby not merely fitting the pitch response, but also isolating the fluid–structure coupling pathway responsible for the observed spectral sidebands and 2 : 1 lock-in.

Modulation of the aerodynamic damping is neglected in (4.5) since the damping force is scaled by the small parameter $\mu$ . Any parametric variation in damping would scale as $\mathcal{O}(\mu \varepsilon )$ , rendering it a significantly weaker coupling mechanism compared to the dominant stiffness modulation of order $\mathcal{O}(\varepsilon )$ , especially considering the complexity of damping modulation that would contribute to the nonlinear state-dependent damping term, $-\mu \omega _0\bigl [1+\varepsilon _{g}\cos (\omega _g t)\bigr ]\bigl (1-x^2\bigr )\dot {x}$ . Thus, instead of a laborious analytical derivation, the damping modulation will be studied using numerical results in § 4.4.

Equations of this form, a self-excited VdP oscillator with a parametrically modulated stiffness, are classical canonical models for 2 : 1 parametric resonance and Mathieu-type synchronisation, and have been analysed extensively in nonlinear dynamics studies, e.g. Pandey et al. (Reference Pandey, Rand and Zehnder2008) and Veerman & Verhulst (Reference Veerman and Verhulst2009), where the emphasis is on a systematic analytical classification of periodic and quasi-periodic responses over broad parameter ranges. In the present study, the same canonical structure is adopted, but with a different emphasis. The derivations developed below are deliberately minimal and tailored to the present aeroelastic problem: only the spectral and averaged amplitude-phase relations are retained to interpret the observed sidebands and 2 : 1 lock-in band in stall flutter under longitudinal gusts. In this sense, the analysis below is intended not as a complete mathematical classification of the M-VdP dynamics, but as a problem-driven reduction tailored to interpret the preceding experimental observations.

4.2. Spectral content analysis away from lock-in

To understand the system’s spectral signature, the autonomous case is first considered by setting the parametric excitation to zero ( $\varepsilon =0$ ). The model reduces to the classical VdP oscillator, which is known to exhibit a stable limit cycle. The carrier frequency of this autonomous oscillation, $\omega _\theta$ , is given by

(4.6) \begin{align} \omega _\theta =\omega _0+\mathcal{O}(\mu ^2). \end{align}

This represents the primary response frequency, which experiences a minor, higher-order shift from the linear natural frequency $\omega _0$ . The nonlinear damping term $-\mu (1-x^2)\dot {x}$ is an odd function of $x$ , resulting in a limit cycle $x_0(t)$ that is, to leading order, half-wave symmetrical. This symmetry suppresses the even Fourier components. Consequently, the autonomous response spectrum consists predominantly of the fundamental frequency $\omega _\theta$ and its odd harmonics:

(4.7) \begin{align} x_0(t)=c_1\cos (\omega _\theta t)+c_3\cos (3\omega _\theta t)+c_5\cos (5\omega _\theta t)+\cdots, \quad c_{2k}\approx 0. \end{align}

Now the parametric excitation is introduced by treating $\varepsilon$ as a small perturbation parameter. A standard perturbation expansion for the solution $x(t)$ is employed:

(4.8) \begin{align} x(t)=x_0(t)+\varepsilon\, x_1(t)+\varepsilon ^2\, x_2(t)+\cdots . \end{align}

Substituting this expansion into (4.5), and collecting terms at the first order $\mathcal{O}(\varepsilon )$ , yields the equation for the first correction $x_1(t)$ . With the linear operator $\mathcal L \equiv ( {{\rm d}^2}/{{\rm d}t^2})+\omega _0^2$ , this equation is

(4.9) \begin{align} \mathcal L\,x_1 = -\,\omega _0^2 \cos (\omega _g t)\,x_0. \end{align}

The right-hand side of (4.9) acts as a multiplicative forcing term. This term causes frequency mixing between the gust, $\cos (\omega _g t)$ , and the autonomous response, $x_0(t)$ . Applying the trigonometric identity $\cos \alpha \cos \beta =( 1/2)[\cos (\alpha +\beta )+\cos (\alpha -\beta )]$ , each odd carrier harmonic $m\omega _\theta$ in $x_0$ interacts with $\omega _g$ to generate a pair of first-order sidebands at

(4.10) \begin{align} \omega ^{(1)}_{m,\pm } = m\omega _\theta \pm \omega _g, \quad m=1,3,5,\ldots . \end{align}

The magnitude of these sidebands scales linearly with the perturbation, $\propto \varepsilon$ .

This recursive process continues at higher orders. At $\mathcal{O}(\varepsilon ^2)$ , the first-order solution $x_1$ , which contains the sidebands from (4.10), is in turn multiplied by the parametric term $\cos (\omega _g t)$ :

(4.11) \begin{align} \mathcal L\,x_2 = -\omega _0^2 \cos (\omega _g t)\,x_1 + \text{(weak VdP cross-terms)}. \end{align}

This secondary mixing generates the second-order sidebands:

(4.12) \begin{align} \omega ^{(2)}_{m,\pm } = m\omega _\theta \pm 2\omega _g, \quad m=1,3,5,\ldots . \end{align}

These components have magnitudes scaling as $\propto \varepsilon ^2$ . Generalising to $\mathcal{O}(\varepsilon ^n)$ , this parametric cascade generates a full lattice of sidebands at $\omega ^{(n)}_{m,\pm } = m\omega _\theta \pm n\omega _g$ .

Thus, for the non-resonant case (away from lock-in), the predicted spectral content is characterised by three key features.

1. (i) Odd carrier harmonics: the fundamental frequency $\omega _\theta$ and its odd multiples ( $m=1,3,5,\ldots$ ), originating from the underlying VdP dynamics.

2. (ii) Parametric sideband lattice: a symmetrical lattice of sidebands, $\omega =m\omega _\theta \pm n\omega _g$ ( $n=1,2,\ldots$ ) surrounding each carrier harmonic. The amplitude of the $n$ th-order sideband is expected to scale as $\propto \varepsilon ^{\,n}$ .

3. (iii) Absence of gust-frequency peak: a critical distinction of this parametric model is the absence of an isolated spectral peak at the gust frequency $\omega _g$ itself. Because the excitation is multiplicative (not additive), the $\omega _g$ component only appears as part of the mixed sideband pairs (e.g. $\omega _\theta \pm \omega _g$ ) unless an external forcing at $\omega _g$ is present.

4.3. Lock-in band

Introducing non-dimensional time $\tau =\omega _0 t$ and frequency ratio $\varOmega =\omega _g/\omega _0$ to (4.5) yields the non-dimensional M-VdP model:

(4.13) \begin{align} x''-\mu \big(1-x^2 \big)\,x'+\bigl [1+\varepsilon \cos (\varOmega \tau )\bigr ]x=0, \end{align}

where $(\,)'$ denotes the derivative with respect to $\tau$ . This non-dimensional framework facilitates the application of the averaging method and the subsequent analysis of lock-in.

4.3.1. Averaging procedure

Considering the principal parametric resonance $\varOmega \approx 2$ , let $\delta =\varOmega -2=\mathcal O(\varepsilon )$ represent the small detuning from exact resonance. To capture the dynamics, we recognise two distinct time scales: the ‘fast time’ $\tau$ , which tracks the natural oscillation of the structure, and the ‘slow time’, over which the amplitude and phase modulation evolve due to the weak parametric perturbation. The method of averaging (Sanders et al. Reference Sanders, Verhulst and Murdock2007) is employed to separate these scales. We adopt the near-harmonic representation with slowly varying amplitude $r(\tau )$ and phase $\phi (\tau )$ :

(4.14) \begin{align} x(\tau )=r(\tau )\cos (\tau +\phi ),\quad \psi (\tau )=\tau +\phi -\frac {\varOmega }{2}\tau , \end{align}

where $\psi$ is the phase difference, or mismatch, between the oscillator’s phase $(\tau +\phi )$ and the effective parametric forcing phase $(\varOmega \tau /2)$ .

Assume $r'=O({\mu },\varepsilon )$ and $\phi '=O({\mu },\varepsilon )$ ; then

(4.15) \begin{align} \begin{aligned} x' &= r'\cos (\tau +\phi ) - r(1+\phi ')\sin (\tau +\phi ), \\ x'' &\approx -r(1+\phi ')^2\cos (\tau +\phi ) - 2r'(1+\phi ')\sin (\tau +\phi ), \end{aligned} \end{align}

where $r''$ and $\phi ''$ are higher-order and are dropped. Substituting (4.14) and (4.15) into (4.13), and applying the standard averaging procedure – projecting onto the basis functions $\cos (\tau +\phi )$ and $\sin (\tau +\phi )$ and integrating over one cycle ( $2\pi$ ) of the fast time $\tau$ – filters out the rapidly oscillating terms, and isolates the secular terms that drive the slow dynamics; see classical nonlinear dynamics textbooks such as Nayfeh & Mook (Reference Nayfeh and Mook1979), Strogatz (Reference Strogatz1994) and Sanders et al. (Reference Sanders, Verhulst and Murdock2007) for more details. Finally, utilising the kinematic relationship $\psi ' = \phi ' - \delta /2$ , the system is transformed into the variables $r$ and $\psi$ . This yields the so-called ‘slow-flow’ equations (a dynamical systems term referring to the flow of trajectories in the reduced phase space, rather than physical fluid motion), which govern the long-term evolution of the system:

(4.16) \begin{eqnarray} r' &=& \frac {1}{2}{\mu }r - \frac {1}{8}{\mu }r^3 + \frac {1}{4}\varepsilon r \sin (2\psi ), \\[-28pt] \nonumber \end{eqnarray}

(4.17) \begin{eqnarray} \psi ' &=& -\frac {1}{2}\delta + \frac {1}{4}\varepsilon \cos (2\psi ). \\[0pt] \nonumber \end{eqnarray}

Equations (4.16) and (4.17) constitute the coupled slow-flow system governing the interaction between stall flutter and longitudinal gusts. Physically, the first two terms in (4.16), $({1}/{2}){\mu }r$ and $- ({1}/{8}){\mu }r^3$ , represent the autonomous flutter mechanism, specifically, the competition between linear negative damping (energy input) and cubic nonlinear damping (saturation). The terms proportional to $\varepsilon$ in both equations capture the impact of the parametric excitation. Synchronisation, or lock-in, occurs when the parametric forcing is strong enough to overcome the frequency detuning $\delta$ , resulting in a stationary amplitude and a constant phase difference.

4.3.2. Lock-in boundary

Steady (phase-locked) solutions satisfy $r'=\psi '=0$ , hence

(4.18) \begin{align} \cos (2\psi )=\frac {2\delta }{\varepsilon },\quad \sin (2\psi )=\frac {1}{\varepsilon }\left(\frac {{\mu }}{2}r^2-2{\mu }\right) \!, \end{align}

and eliminating $\psi$ gives the existence boundary

(4.19) \begin{align} \varepsilon ^2=4\delta ^2+\left(\frac {{\mu }}{2}r^2-2{\mu }\right)^2 \!. \end{align}

For the autonomous VdP LCO, $( {{\mu }}/{2})\,r^2-2{\mu }=0$ so $r_0=2$ ; evaluating at $r=r_0$ removes the amplitude term and yields the half-width in detuning:

(4.20) \begin{align} |\delta |_{max }=\frac {1}{2}\varepsilon \quad \Longleftrightarrow \quad |\varOmega -2|_{max }=\frac {1}{2}\varepsilon . \end{align}

Mapping back to dimensional frequencies ( $\varOmega =\omega _g/\omega _0$ ) gives

(4.21) \begin{align} \fbox{$ |\omega _g-2\omega _0|_{max }=\frac {1}{2}\varepsilon \omega _0 $}. \end{align}

This result indicates that the synchronisation bandwidth scales linearly with the modulation strength $\varepsilon$ and the natural frequency $\omega _0$ . Importantly, to leading order, this width is independent of the autonomous limit-cycle amplitude $r_0$ (and thus the baseline pitch angle $\overline {\theta }_0$ ), a feature that distinguishes this parametric mechanism from the classical forced VdP oscillator. These analytical predictions are next verified against numerical solutions of the full system.

4.4. Analytico-numerical verification

All derivations in §§ 4.1–4.3 assume small $\varepsilon$ . To test the predictions and delineate the validity range of $\varepsilon$ , the M-VdP model in (4.5) is solved numerically. For the numerical simulations, (4.5) is integrated in MATLAB using the ode45 solver with relative and absolute tolerances $10^{-8}$ and $10^{-10}$ , respectively. The maximum solver step is limited to $0.002\,\mathrm{s}$ , and the solution is output at $\Delta t=0.001\,\mathrm{s}$ , corresponding to $f_s=1000\,\mathrm{Hz}$ . Each simulation is run for $180\,\mathrm{s}$ , and the first $30\,\mathrm{s}$ are discarded as transient. The final $150\,\mathrm{s}$ are used to compute the statistically steady-state amplitude, response frequency and PSD. The response frequency is estimated from upward zero-crossings of the de-meaned pitch signal. The PSD values are computed using the same Welch-type averaged FFT periodogram as for the experimental data, yielding frequency-bin spacing $0.00763\,\mathrm{Hz}$ .

The experimental case in figure 9 is reproduced with $f_g=2.5\,\mathrm{Hz}$ , $f_{\theta ,0}=2.93\,\mathrm{Hz}$ , $\bar {\theta }_0=41^{\circ }$ and $\varepsilon =0.2$ . The damping parameter $\mu$ dictates the slow time scale of amplitude evolution (see (4.16)), which has a negligible effect on the steady-state frequency and amplitude (Nayfeh & Mook Reference Nayfeh and Mook1979; Strogatz Reference Strogatz1994). Consequently, $\mu$ mainly affects the transient period, which is set to 0.1 throughout all numerical simulations. The time history in figure 16(a) shows stronger amplitude modulation than the experimental result in figure 9(a). In experiments, there is an amplitude convergence of stall flutter (approximately $41^\circ$ ) shown in figures 3(a) and 12(b), which naturally limits the beating amplitude. The amplitude convergence is a higher-order aeroelastic nonlinearity that is not captured in this minimal model. The PSD in figure 16(b) exhibits primary carriers at $m f_{\theta }$ (odd $m$ ) together with modulation sidebands at $\lvert m f_{\theta } \pm n f_{g}\rvert$ ( $n=1,2,3,\ldots$ ). Selected components are labelled and summarised as $f_i=m f_{\theta }+n f_{g}$ in table 3. These sidebands match the experimental observation in figure 9(b), and support the derivation in § 4.2.

Table 3.

Representative sidebands expressed as $f_i=mf_{\theta }+nf_{g}$ , corresponding to the labels in figure 16(b).

Table 3. Long description

A table titled ‘Coefficient’ compares values for two functions, m(fθ) and n(fg), across seven coefficients: fb, f1, f2, f3, f4, f5, and f6. The table has 3 rows and 8 columns. Row 1: Coefficient, fb, f1, f2, f3, f4, f5, f6. Row 2: m(fθ), +1, -1, +3, +1, +3, +1, +3. Row 3: n(fg), -1, +2, -2, +1, -1, +2, +1.

Figure 16.

Numerical M-VdP results for $f_g=2.5\,\mathrm{Hz}$ , $f_{\theta ,0}=2.93\,\mathrm{Hz}$ , $\bar {\theta }_0=41^{\circ }$ and $\varepsilon =0.2$ : (a) time history; (b) PSD with labelled sidebands (see table 3), reproducing experimental results in figure 9.

A sweep in $f_g$ using the M-VdP model (figure 17) reveals no 1 : 1 lock-in near $f_g\approx f_{\theta ,0}$ , but a clear 2 : 1 synchronisation band, consistent with the classical Mathieu-type parametric resonance and also the experimental observations shown in figure 13. Note that the M-VdP parameters, such as $f_{\theta ,0}$ and $\bar {\theta }_0$ , are adjusted to phenomenologically match the experimental sweep capturing the 2 : 1 lock-in zone. In addition, the modulation strength is physically calibrated to $\varepsilon = 0.26$ , derived from the flow measurement using the X-hot-wire probe that indicated $u_{g,\textit{max}} = 1.3\,\mathrm{m\,s^{-1}}$ at $f_g = 3.9$ Hz (not the ideally central frequency but within the lock-in band). With this physically constrained parameter, the numerical sweep yields predicted lock-in bandwidth 0.48 Hz. This result demonstrates a remarkable three-way quantitative agreement: it matches the analytical prediction 0.486 Hz derived in (4.21), and aligns closely with the experimentally observed bandwidth approximately 0.5 Hz shown in figure 13. This convergence of analytical, numerical and experimental scales provides strong support for the validity of the proposed M-VdP model, confirming that the observed lock-in is indeed driven by the parametric stiffness modulation mechanism.

The $\bar {\theta }$ trend around the 2 : 1 lock-in band is symmetric for the M-VdP model, which is typical resonance response, but asymmetric for the measured results; see figure 13. This discrepancy suggests that subtle physical mechanisms (e.g. higher-harmonic aerodynamic contributions) lie beyond the scope of these minimal models. In particular, Pandey et al. (Reference Pandey, Rand and Zehnder2008) reported that a forced M-VdP–Duffing model can produce asymmetric amplitudes of lock-in bands, similar to figure 13, but the added cubic stiffness (Duffing term) shifts the oscillation frequency from $f_{\theta ,0}$ out of lock-in, contrary to the experimental observation. Nevertheless, such extensions point to a promising pathway for refining the present minimal model in future studies.

Figure 17.

Numerical M-VdP results for a frequency sweep of the parametric excitation (covering both 1 : 1 and 2 : 1 regions), with $f_{\theta , 0} = 1.87$ Hz, $\bar {\theta }_0=29^\circ$ and $\varepsilon = 0.26$ , reproducing experimental results in figure 13.

Figure 17 Long description

Three line graphs are presented. Panel A: The line graph shows the frequency (f_theta) in Hertz on the vertical axis ranging from 1.6 to 2.1 Hz, plotted against the excitation frequency (f_g) in Hertz on the horizontal axis ranging from 1 to 5 Hz. The data points form a line with a noticeable dip around 2 Hz and a peak around 4 Hz. Panel B: The line graph displays the mean angle of attack (theta_bar) in degrees on the vertical axis ranging from 20 to 40 degrees, plotted against the excitation frequency (f_g) in Hertz on the horizontal axis ranging from 1 to 5 Hz. The data points form a line with a peak around 2 Hz and a dip around 4 Hz. Panel C: The line graph illustrates the standard deviation of the angle of attack (sigma_theta) in degrees on the vertical axis ranging from 0 to 0.15 degrees, plotted against the excitation frequency (f_g) in Hertz on the horizontal axis ranging from 1 to 5 Hz. The data points form a line with peaks around 2 Hz and 4 Hz.

For the present M-VdP model, (4.20) predicts that the non-dimensional parametric 2 : 1 bandwidth depends only on the longitudinal-gust strength, represented by $\varepsilon$ , and is essentially independent of the baseline amplitude. These predictions are confirmed by figures 18 and 19, i.e. varying $\bar {\theta }_0$ leaves the band unchanged, whereas increasing $\varepsilon$ expands it. The dimensional widths of the lock-in bands are $0.18$ , $0.36$ and $0.56\,\mathrm{Hz}$ for $\varepsilon =0.1$ , $0.2$ and $0.3$ , respectively, in agreement with (4.21). For $\varepsilon =0.3$ , the response frequency deviates more noticeably from $f_{\theta ,0}$ , reflecting stronger nonlinear effects of the greater modulation.

Figure 18.

Effect of baseline amplitude around the 2 : 1 region: $f_{\theta ,0}=1.87\,\mathrm{Hz}$ , $\varepsilon =0.2$ and $\bar {\theta }_0=19^{\circ },29^{\circ },39^{\circ }$ .

Figure 19.

Effect of stiffness modulation strength around the 2 : 1 region: $f_{\theta ,0}=1.87\,\mathrm{Hz}$ , $\bar {\theta }_0=29^{\circ }$ and $\varepsilon =0.1,0.2,0.3$ .

Revisiting the 1 : 1 region in figures 12(a), 13 (experiment) and 17 (M-VdP), subtle frequency fluctuations near $f_g\approx f_{\theta ,0}$ hint at the presence of a weak interaction. To determine whether stronger parametric modulation could induce 1 : 1 synchronisation in the M-VdP model, numerical sweeps are conducted with modulation strength $\varepsilon$ ranging from 0.05 to 0.6. Figure 20 confirms that while low modulation strengths result in quasi-periodic beating (implied by non-zero $\sigma _\theta$ ), a sufficiently large $\varepsilon$ yields definitive 1 : 1 lock-in. This high-modulation synchronisation is a strictly nonlinear phenomenon not captured by the first-order perturbation analysis, thus outside the scope of §§ 4.1–4.3.

Figure 21 presents the global synchronisation map in the $(f_g,\varepsilon)$ plane. These boundaries are determined numerically using the frequency matching criteria $f_\theta = f_g$ for the 1 : 1 zone, and $f_\theta = ({1}/{2})f_g$ for the 2 : 1 zone. The resulting synchronisation regions, classically referred to as Arnold tongues, after V.I. Arnold (Pikovsky et al. Reference Pikovsky, Rosenblum and Kurths2001), exhibit two qualitatively distinct behaviours. First, the primary 2 : 1 parametric lock-in zone forms a sharp V-shaped wedge originating from $\varepsilon =0$ . The numerical boundaries in this region demonstrate excellent agreement with the analytical prediction derived in (4.21), verifying the validity of the averaging method used to derive the lock-in criteria. Second, the 1 : 1 lock-in zone exhibits a detached tongue structure. Unlike the 2 : 1 resonance, the 1 : 1 synchronisation requires a finite threshold of modulation strength ( $\varepsilon \gt 0.2$ ) to manifest. Below this threshold, the parametric energy input is insufficient to overcome the detuning. This contrast underscores the different physical mechanisms at play: the 2 : 1 resonance is the primary instability of the Mathieu-type term, whereas the 1 : 1 lock-in is a higher-order nonlinear effect. Experimentally, generating the large modulation amplitudes required to observe 1 : 1 lock-in is beyond the current apparatus capabilities, suggesting a compelling direction for future work.

Figure 20.

Numerical frequency sweep around the 1 : 1 region for varying stiffness modulation strengths: $f_{\theta ,0}=1.87\,\mathrm{Hz}$ , $\bar {\theta }_0=29^{\circ }$ and $\varepsilon =0.05, 0.1, 0.2,0.4,0.6$ .

Figure 20 Long description

Three line graphs illustrate the relationship between frequency and various parameters for different stiffness modulation strengths. Panel A: The line graph shows the frequency f_theta (Hz) on the vertical axis and the frequency fg (Hz) on the horizontal axis. Different lines represent varying stiffness modulation strengths (epsilon values of 0.05, 0.1, 0.2, 0.4, and 0.6), indicated by different colors and symbols. Panel B: The line graph displays the average angle theta (degrees) on the vertical axis and the frequency fg (Hz) on the horizontal axis. The same color and symbol scheme is used to represent different epsilon values. Panel C: The line graph shows the standard deviation of the angle sigma_theta (degrees) on the vertical axis and the frequency fg (Hz) on the horizontal axis, again using the same color and symbol scheme for different epsilon values. Each panel reveals how the parameters change with varying stiffness modulation strengths.

Figure 21.

Map of synchronisation regions (Arnold tongues) for the M-VdP model with stiffness modulation (4.13) ( $f_{\theta ,0} = 1.87$ Hz, $\overline {\theta }_0 = 29^\circ$ ): 2 : 1 parametric lock-in (blue) and 1 : 1 lock-in (red, emerging at higher modulation strengths). The V-shaped black dashed lines represent the analytical prediction for the 2 : 1 boundary derived in (4.21).

Finally, to justify the exclusion of damping modulation in the primary model (4.13), numerical simulations are conducted using an alternative formulation where the gust modulates the damping rather than the stiffness. Assuming that the aerodynamic damping inherits the parametric dependence through $D_a(t) \propto U(t) = U_\infty [1+\varepsilon _g \cos (\omega _g t)]$ , (4.13) is modified as

(4.22) \begin{align} x'' - {\mu }\bigl [1+\varepsilon _g\cos (\varOmega \tau )\bigr ](1-x^2)\,x' + x = 0. \end{align}

The resulting synchronisation boundaries, shown in figure 22, reveal a lock-in bandwidth that is significantly narrower than that of the stiffness modulation case (figure 21). Considering the scaling relation $\varepsilon _g = \varepsilon /2$  for the stiffness modulation model, the difference in bandwidth is even more pronounced than it appears. This reduced bandwidth is consistent with the scaling analysis in § 4.1: the parametric damping term scales as $\mathcal{O}({\mu }\varepsilon )$ , rendering it a significantly weaker coupling mechanism than the stiffness modulation of order $\mathcal{O}(\varepsilon )$ . Furthermore, the numerical amplitude trend for damping modulation (figure 23) exhibits a symmetric profile, but with a slight increase at lock-in boundaries, and a distinct dip in $\bar {\theta }$ (lower than $\bar {\theta }_0$ ) near the resonance condition $f_g \approx f_{\theta ,0}$ , which are contrary to the experimental observations (slight drop at the lower boundary and peak at the upper boundary, with the peaking level higher than $\bar {\theta }_0$ ). These discrepancies indicate that damping modulation might not capture the correct physics of the gust interaction, and the dominant physics of longitudinal gust excitation on stall flutter originates from the parametric stiffness modulation.

Figure 22.

Map of synchronisation regions (Arnold tongues) for the M-VdP model with damping modulation (4.22) ( $f_{\theta ,0} = 1.87$ Hz, $\overline {\theta }_0 = 29^\circ$ ): 2 : 1 parametric lock-in (blue) and 1 : 1 lock-in (red, fundamentally invisible).

Figure 23.

Effect of damping modulation strength: $f_{\theta ,0}=1.87\,\mathrm{Hz}$ , $\bar {\theta }_0=29^{\circ }$ and $\varepsilon _g = 0.1,0.2,0.4,0.8$ . Unlike the experimental observations (figure 13) and stiffness modulation results (figure 17), the response amplitude $\bar {\theta }$ exhibits a distinct dip inside the lock-in band.