1. Introduction
Aeroelastic flutter is a critical instability in lightweight and flexible structures that interact with unsteady flows, characterised by self-excited limit-cycle oscillations (LCOs) that can compromise structural integrity (Dowell Reference Dowell2021). Controlling these oscillations is essential to improve aerodynamic performance and prevent undesirable responses (Jonsson et al. Reference Jonsson, Riso, Lupp, Cesnik, Martins and B.2019). However, predicting stability transitions remains challenging due to the nonlinear nature of fluid–structure interactions and their sensitivity to flow separation and vortex shedding (Schuster, Liu & Huttsell Reference Schuster, Liu and Huttsell2003).
Classical theories assume thin airfoils and inviscid, attached flows (Theodorsen Reference Theodorsen1935; Bisplinghoff, Ashley & Halfman Reference Bisplinghoff, Ashley and Halfman2013), limiting their applicability in regimes dominated by unsteady separation. In such conditions, laminar separation flutter and stall-induced LCOs can emerge from nonlinear interactions between separated flow and structural motion (Lee, Price & Wong Reference Lee, Price and Wong1999; Barnes & Visbal Reference Barnes and Visbal2019). Recent advances in flutter prediction include multi-fidelity models (Thelen, Leifsson & Beran Reference Thelen, Leifsson and Beran2020), frequency-domain identification (Simiriotis & Palacios Reference Simiriotis and Palacios2023) and data-driven techniques (Hickner et al. Reference Hickner, Fasel, Nair, Brunton and Brunton2023; Guo et al. Reference Guo, Li, Zhou, Ma and Wang2025). Complementary efforts in active flutter suppression employ high-bandwidth actuators and adaptive controllers to stabilise flutter-prone modes and mitigate LCOs (Li, Xiang & Guo Reference Li, Xiang and Guo2011; Livne Reference Livne2018; Chai et al. Reference Chai, Gao, Ankay, Li and Zhang2021).
However, the impact of localised structural variations remains poorly understood, particularly in systems exhibiting subcritical bifurcations, where small changes can trigger instability (Castravete & Ibrahim Reference Castravete and Ibrahim2008; Khodaparast, Mottershead & Badcock Reference Khodaparast, Mottershead and Badcock2010; Beran, Stanford & Schrock Reference Beran, Stanford and Schrock2017). Energy maps have proven useful for understanding the nonlinear aeroelastic response to gusts (Menon & Mittal Reference Menon and Mittal2019,Reference Menon and Mittal2020); however, their dependence on steady-state responses from forced simulations across a range of parameters limits their utility for real-time flutter characterisation.
To address these limitations, phase reduction offers a low-dimensional tractable framework by modelling unsteady periodic flows as nonlinear oscillators (Taira & Nakao Reference Taira and Nakao2018; Monga et al. Reference Monga, Wilson, Matchen and Moehlis2019). Recent advances include timing control of vortex-shedding dynamics (Nair et al. Reference Nair, Taira, Brunton and Brunton2021), adjoint-based formulations for rapid wake synchronisation (Godavarthi, Kawamura & Taira Reference Godavarthi, Kawamura and Taira2023) and data-driven latent manifold approaches for suppressing transient gust responses (Fukami, Nakao & Taira Reference Fukami, Nakao and Taira2024).
Unlike canonical oscillators, aeroelastic systems can exhibit abrupt transient instabilities that lead to flutter, often before reaching a stable periodic orbit. In this study, we apply phase-based methods to characterise and control the onset of such instabilities. Using two-dimensional incompressible flow over a NACA0015 airfoil at a low Reynolds number, we introduce time-localised perturbations during the free-response oscillation cycle and show that their impact is highly phase-dependent. We demonstrate how this sensitivity can be exploited to suppress flutter through targeted, phase-based control (PBC) strategies. Notably, despite the fluid and structural modes exhibiting distinct frequencies, a single-oscillator phase-based model grounded in structural observables was sufficient to capture and control the instability.
While the underlying methods such as phase reduction, energy-transfer metrics and dynamic mode decomposition (DMD) are well established, our study uniquely applies them to uncover and exploit phase-dependent flutter sensitivity in an aeroelastic system. In particular, we integrate stiffness-triggered instability diagnostics with heaving-based actuation to design an energy-optimal control strategy that disrupts phase synchronisation and suppresses aeroelastic flutter. An overview of the present work is shown in figure 1. Through direct numerical simulations, phase-sensitivity analysis and targeted perturbation tests, we design a PBC strategy that suppresses the flutter instability.
Figure 1. Overview of the present work. (a) Direct numerical simulations (DNS) of the baseline free-pitching dynamics of a NACA0015 airfoil. (b) Phase-localised stiffness perturbations are applied, and the resulting pitch displacement trajectories are recorded. (c) The phase-sensitivity function is computed, revealing critical phases associated with flutter. (d) A PBC strategy is designed based on heave phase-sensitivity function and its gradient. (e) The control input is implemented in DNS, demonstrating suppression of flutter.
2. Methodology
We perform two-dimensional simulations of the aeroelastic response of a NACA0015 airfoil immersed in an incompressible flow using the immersed boundary projection method (Taira & Colonius Reference Taira and Colonius2007; Goza & Colonius Reference Goza and Colonius2017). The fluid motion is governed by the incompressible Navier–Stokes equations in non-dimensional form:
\begin{equation} \frac {\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla}\!\boldsymbol{u} = -\boldsymbol{\nabla}\!p + \frac {1}{Re} {\nabla} ^2 \boldsymbol{u}, \quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0, \end{equation}
where
$\boldsymbol{u}$
is the velocity field,
$p$
is pressure and
$Re = \rho U_\infty c / \mu$
is the Reynolds number, with
$U_\infty$
,
$\rho$
,
$c$
and
$\mu$
denoting the free-stream velocity, density, chord length and dynamic viscosity, respectively. The immersed boundary projection method discretises these equations on a staggered Cartesian grid, employing a Crank–Nicolson scheme for viscous terms and a second-order Adams–Bashforth scheme for convective terms. No-slip boundary conditions are enforced via Lagrange multipliers applied at Lagrangian boundary points, with interpolation and spreading handled via regularised delta functions. A nested multi-domain grid with a fast Poisson solver (Colonius & Taira Reference Colonius and Taira2008) is used to efficiently resolve near- and far-field interactions.
Figure 2. Free-pitching dynamics of a NACA0015 airfoil. (a) Schematic of the computational set-up showing the airfoil mounted at an elastic axis located at
$x_e = 0.33c$
from the leading edge. (b) Nested multi-domain mesh layout used for simulations, with finest resolution near the airfoil. (c) Validation of gust-induced LCOs in pitch angle
$\theta (t)$
, reproducing aeroelastic responses consistent with Menon & Mittal (Reference Menon and Mittal2020).
The airfoil is free to pitch about an elastic axis located at
$x_e = 0.33c$
, corresponding to 33 % of the chord length from the leading edge (figure 2
a). Its structural dynamics is modelled as a single-degree-of-freedom linear torsional oscillator governed by
\begin{equation} I^* \ddot {\theta } + k^*(\theta - \theta _0) = C_M, \end{equation}
where
$\theta$
is the pitch angle,
$\theta _0$
is the equilibrium angle,
$C_M$
is the dimensionless aerodynamic moment and
$I^*$
and
$k^*$
are the non-dimensional moment of inertia and stiffness, respectively.
We adopt parameter values consistent with those used by Menon & Mittal (Reference Menon and Mittal2020), who identified an aeroelastic regime exhibiting strong nonlinear fluid–structure coupling and the onset of gust-induced flutter. These choices are further supported by the extensive parametric studies conducted in Menon & Mittal (Reference Menon and Mittal2019). The selected parameter values and their rationale are summarised in table 1.
Table 1. Summary of selected parameters and rationale.
The structural natural frequency of the aeroelastic system is given by
$ f_s = (1/2\pi )\sqrt {k^*/I^*} \approx 0.19$
, with a corresponding time period
$ T$
and reduced velocity
$ U^* = 1/f_s \approx 5.26$
. Prior studies (Menon & Mittal Reference Menon and Mittal2019) have shown that a lower structural natural frequency, along with a more aft location of the elastic axis, increases the system’s susceptibility to flutter onset by enhancing fluid–structure coupling.
The computational domain spans
$x/c \in [-16, 16]$
and
$y/c \in [-16, 16]$
. The finest grid near the airfoil has a uniform spacing
$\Delta x = \Delta y = 0.0055c$
, and a fixed time step
$\Delta t = 0.001$
is used for time integration. The near-body mesh and nested multi-domain layout are shown in figures 2(a) and 2(b), respectively. The set-up is validated by reproducing gust-induced LCO responses as seen in figure 2(c), consistent with Menon & Mittal (Reference Menon and Mittal2020), confirming the nonlinear amplification behaviour near the onset of flutter.
To investigate the timing-dependent sensitivity of the aeroelastic system, we adopt the observable-based phase-reduction framework of Taira & Nakao (Reference Taira and Nakao2018), which projects the high-dimensional fluid–structure dynamics onto a scalar phase variable defined along the periodic orbit. Specifically, we define the phase
$ \phi$
on the
$ \dot {\theta }$
–
$ \theta$
plane, where
$ \theta$
denotes the airfoil pitch angle and
$ \dot {\theta }$
is the angular velocity. This observable-based representation traces the aeroelastic limit cycle with period
$ T$
and provides a smooth, monotonically increasing phase coordinate over the oscillation cycle, as shown in the polar plot of figure 3(a). The bottom panel of figure 3(a) shows the corresponding time series of the aerodynamic moment coefficient
$ C_M$
, with vertical red lines marking the bounds of one period of the pitch oscillation
$T$
. Within this interval,
$ C_M$
exhibits higher-frequency oscillations compared with
$ \theta$
, with a dominant frequency of
$ f \approx 0.67$
, which is also observed when the airfoil is held static. This frequency mismatch between
$ \theta$
and
$ C_M$
in the free response indicates a lack of synchronisation between structural motion and aerodynamic forcing, thereby preventing sustained energy transfer and averting flutter in the baseline configuration. As a result, although
$ C_M(t)$
exhibits high-frequency fluctuations and is not strictly periodic, the phase definition based on structural observables
$ (\theta , \dot {\theta })$
remains consistent for characterising phase sensitivity. Figure 3(b) presents representative vorticity fields at two characteristic phases, revealing the evolution of coherent vortical structures around the airfoil as it traverses the limit cycle.
Figure 3. Phase-based analysis of aeroelastic response to impulsive stiffness perturbations. (a) Definition of phase
$\phi$
on the
$\dot {\theta }$
–
$\theta$
plane (top) and time series of the aerodynamic moment coefficient
$C_M$
(bottom), with red vertical lines indicating one period of pitch oscillation. (b) Instantaneous vorticity fields at representative phases. (c,d) Phase-sensitivity function
$Z(\phi )$
estimated from asymptotic pitch-angle shifts following negative (decrease in
$k^*$
) and positive (increase in
$k^*$
) stiffness impulses, with impulses of different strengths applied in the positive stiffness case. Dashed lines mark phases where perturbations lead to divergent flutter-like behaviour. (e) Normalised energy extraction
$E$
computed over 18 post-perturbation cycles, showing increased energy transfer at critical phases. (f,g) Exponential growth rate
$\lambda$
of the envelope of
$C_M(t)$
, with positive values indicating instability. Also shown are the Arnold tongue lock-on regimes (highlighted in pink) corresponding to negative and positive stiffness phase-sensitivity functions in insets of (f) and (g), respectively.
Under an external perturbation applied at time
$t_0$
, the phase dynamics evolves as
\begin{equation} \dot {\phi }(t) = 2\pi f_s + \epsilon Z(\phi )\delta (t - t_0), \end{equation}
where
$Z(\phi )$
is the phase-sensitivity function,
$\epsilon$
is the perturbation magnitude and
$\delta$
is the Dirac delta function centred at the impulse time
$t_0$
(Nakao Reference Nakao2016). Following the impulse-response method of Taira & Nakao (Reference Taira and Nakao2018), we apply Gaussian-shaped impulses at uniformly spaced phases across the oscillation cycle as
\begin{equation} \epsilon \delta (t-t_0) \approx \frac {\beta }{\sqrt {2\pi \sigma }} \exp \left [-\frac {1}{2} \left ( \frac {t - t_0}{\sigma } \right )^2 \right ], \end{equation}
where
$\beta$
and
$\sigma = 0.0019T$
are the strength and width of the impulse introduced at time
$ t_0$
, respectively.
These impulses are introduced by perturbing structural parameters (e.g. stiffness) or by imposing surging/heaving velocity impulses on the airfoil. In this study, stiffness perturbations are used to probe phase-based flutter sensitivity due to their direct relevance to structural degradation under cyclic loading. Stiffness variations, arising naturally from material fatigue or geometric non-uniformities, directly affect the system’s natural frequency
$f_s$
, making them effective triggers for phase-dependent stability transitions. The impulse strength
$ \beta$
for stiffness perturbations is selected such that the relative elastic energy input (REI),
\begin{equation} \text{REI} = \frac {\int _0^T \frac {1}{2} \Delta K(t) \boldsymbol{\cdot }\theta (t)^2 \, {\textrm{d}}t}{\int _0^T \frac {1}{2} K_0 \boldsymbol{\cdot }\theta (t)^2 \, {\textrm{d}}t}, \end{equation}
remains sufficiently small to preserve linearity and uphold the assumptions of the phase-reduction framework. For the perturbation amplitude used in this study,
$\text{REI} \approx \pm 0.00175$
.
The phase change induced by each perturbation is evaluated relative to the unperturbed baseline by tracking the observable
$\theta (t)$
. The phase-sensitivity function is then approximated as
$Z(\phi ) \approx \Delta \phi /\beta $
, where
$\Delta \phi$
denotes the asymptotic phase difference between the baseline and the perturbed trajectories. A positive value of
$Z(\phi )$
implies a phase delay (i.e. the perturbed trajectory lags the baseline), while a negative value indicates a phase advance. This observable-based formulation bypasses the need for full-state measurements and provides an experimentally feasible way to quantify phase sensitivity in high-dimensional aeroelastic systems.
Once
$Z(\phi )$
is known, we design an optimal open-loop control input
$u(t)$
to advance or delay the oscillation phase over a finite time horizon. This strategy leverages the observable-based phase-reduction model to implement timing-sensitive actuation without requiring full-state feedback. The controlled phase dynamics is
\begin{equation} \dot {\phi }(t) = 2\pi f_s + u(t) Z(\phi ), \end{equation}
with the objective of achieving a prescribed phase shift
$\Delta \phi$
while minimising control effort. This leads to a two-point boundary-value problem, which can be solved using standard calculus of variations techniques (Monga et al. Reference Monga, Wilson, Matchen and Moehlis2019; Nair et al. Reference Nair, Taira, Brunton and Brunton2021). We assess the effectiveness of this strategy in suppressing the onset of flutter-like instability.
3. Results
We first examine the effect of impulsive stiffness perturbations across the free-response oscillation cycle. Gaussian-shaped impulses are applied to the stiffness parameter
$k^*$
in (2.2) at 32 uniformly spaced phases, and the resulting system evolution is monitored. Figures 3(c) and 3(d) show the estimated phase-sensitivity function
$Z(\phi )$
for negative (decrease in
$k^*$
making it less stiff impulsively) and positive (increase in
$k^*$
) stiffness perturbations, respectively, computed from asymptotic shifts (after 19 fundamental oscillation cycles) in pitch-angle trajectories as detailed in § 2.
The sensitivity functions reveal distinct regions of phase advance (
$Z(\phi ) \lt 0$
) and phase delay (
$Z(\phi ) \gt 0$
), reflecting how timing of perturbation affects the flow–structure coupling. For negative stiffness perturbations (figure 3
c), four specific phases,
$\phi \approx 0.353\pi$
,
$1.410\pi$
,
$1.472\pi$
and
$1.526\pi$
, lead to divergent responses and trigger flutter onset. Similarly, in the case of positive perturbations (figure 3
d), flutter is initiated when impulses are applied at
$\phi \approx 1.472\pi$
,
$1.656\pi$
,
$1.712\pi$
and
$1.768\pi$
. These critical phases are marked with vertical dashed lines and are associated with sharp variations in
$Z(\phi )$
. For the positive stiffness perturbation case (figure 3
d), we also assess robustness by varying the perturbation magnitude. For REI = 0.00084, all phases remain stable, while the primary case (REI = 0.00175) exhibits instability only at the identified phases above. Error bars for these two cases represent the relatively small standard deviation in the phase-sensitivity function values. Increasing to REI = 0.0026 introduces additional instability points, indicating a critical perturbation threshold between 0.00084 and 0.00175 beyond which flutter onset becomes more widespread.
The observed asymmetry between sensitivity functions for positive and negative perturbations stems from stiffness being a material property. Unlike velocity-based actuation, stiffness perturbations directly alter the natural frequency
$ f_s$
, leading to inherently nonlinear and directionally biased responses. This contrasts with prior studies (Taira & Nakao Reference Taira and Nakao2018), where momentum-based impulses yielded symmetric
$ Z(\phi )$
. This directional asymmetry leads to biased synchronisation behaviour, as seen in the Arnold tongue structures in the insets of figure 3( f,g). The Arnold tongue is constructed by evaluating the phase-coupling function
$ \varGamma (\phi )$
as the convolution of the phase-sensitivity function
$ Z(\phi )$
with a sinusoidal forcing input over one period, following the averaging method described in Taira & Nakao (Reference Taira and Nakao2018). Specifically, lock-on regions (highlighted in pink) appear more pronounced for forcing frequencies higher than the structural natural frequency
$ f_s$
. This suggests that the aeroelastic system is more susceptible to phase synchronisation under faster periodic inputs.
The structure of
$Z(\phi )$
further reveals distinct phase-dependent regimes that align with the dynamics in the
$\dot {\theta }$
–
$\theta$
phase plane (figure 3
a). In the case of negative stiffness perturbations (figure 3
c), four alternating intervals are observed:
$Z(\phi ) \lt 0$
from
$0$
to
$\pi /2$
and again from
$\pi$
to
$3\pi /2$
, indicating phase advance during upstroke and downstroke, respectively; and
$Z(\phi ) \gt 0$
from
$\pi /2$
to
$\pi$
and from
$3\pi /2$
to
$2\pi$
, indicating phase delay near the turning points of pitch motion. These transitions reflect how structural input interacts with unsteady aerodynamic loads depending on the direction of motion. In contrast, the positive perturbation case (figure 3
d) displays two broad lobes of phase delay, peaking near
$\phi = \pi /2$
and
$3\pi /2$
, where pitch angle
$\theta$
reaches its maxima. Although the precise locations of instability differ, both cases demonstrate that flutter onset is strongly localised in phase. For both perturbations,
$Z \rightarrow 0$
at
$\phi = 0$
and
$\pi$
as
$\theta = 0$
at these phases and thus the stiffness perturbation does not modify the structural dynamics of (2.2) as much. We note that the flutter-sensitive phases generally occur during the portion of the oscillation cycle where the airfoil is returning from its peak angle of attack towards equilibrium.
To evaluate the response of the aeroelastic system following impulsive stiffness perturbations, we compute the normalised energy extraction coefficient
$E$
over a transient window starting immediately after the period in which the impulse is applied and spanning 18 fundamental oscillation cycles. This energy metric, adapted from Menon & Mittal (Reference Menon and Mittal2020), is
\begin{equation} E = \frac {1}{(T_{f}-T_i)} \int _{T_i}^{T_{f}} C_M(t) \, \dot {\theta }(t) \, {\textrm{d}}t, \end{equation}
where
$C_M(t)$
is the aerodynamic moment coefficient,
$\dot {\theta }(t)$
is the pitch rate,
$T_i = T$
is the initial time and
$T_f = 19T$
is the final time. The sign and magnitude of
$E$
indicate whether the fluid imparts energy into the structure (positive
$E$
) or extracts energy from it (negative
$E$
) over the transient evolution. As shown in figure 3(e), negative stiffness perturbations (circles) and positive perturbations (triangles) yield large positive
$E$
values at specific phases, aligning with those identified in figure 3(c,d) as prone to flutter onset. These high-energy-extraction phases mark the system’s transition from stable periodic motion to divergent, flutter-like growth. To assess the sensitivity of
$E$
to the chosen evaluation duration, we repeated the calculation with shorter windows with final times of
$12T$
,
$14T$
,
$16T$
and
$18T$
. The standard deviation across these cases is shown as error bars in figure 3(e). The results confirm that the identified high-energy-extraction phases are robust to changes in window length, with variations in
$E$
remaining small.
To quantify the instability, we fit an exponential envelope to the pitch moment coefficient
$C_M(t)$
using
$\max (C_M)(t) \sim C_0 e^{\lambda t}$
, where
$\lambda$
is the instantaneous growth rate. Figures 3(f) and 3(g) show
$\lambda$
for negative and positive stiffness perturbations, respectively. Positive values of
$\lambda$
indicate exponential growth and align with phase regions exhibiting high energy
$E$
and sharp variations in
$Z(\phi )$
, while negative values reflect damping and return to periodic behaviour. The consistency across
$\lambda$
,
$Z(\phi )$
and
$E$
confirms that flutter onset is phase-dependent and can be predicted using phase-based metrics. To assess robustness, the growth rate estimation was repeated with varying fitting windows; the resulting standard deviations, shown as error bars in figure 3(f,g), confirm the relative insensitivity of
$\lambda$
to window selection.
Figure 4. Transient response to impulsive stiffness perturbations at selected phases. (a,d,g,j) Baseline (unperturbed) response; (b,e,h,k) stable response to perturbation at
$\phi = 0.29\pi$
; (c,f,i,l) unstable (flutter) response to perturbation at
$\phi = 0.353\pi$
. (a–c) Phase portraits in the
$\dot {\theta }$
–
$\theta$
plane with unperturbed (blue) and perturbed (red) trajectories, overlaid with instantaneous vorticity fields. (d–f) Time series of pitch angle
$\theta (t)$
and aerodynamic moment coefficient
$C_M(t)$
. Time–frequency spectrograms (g–i) of
$\theta (t)$
and (j–l) of
$C_M(t)$
, computed using continuous wavelet transforms to capture joint temporal and spectral features.
To visualise the transient response to stiffness perturbations, we compare three representative cases in figure 4: (i) the unperturbed baseline (figure 4
a,d,g,j), (ii) a stable response to a perturbation applied at
$\phi = 0.29\pi$
(figure 4
b,e,h,k) and (iii) an unstable (flutter) response resulting from a perturbation at
$\phi = 0.353\pi$
(figure 4
c,f,i,l). Figure 4(a–c) shows phase portraits in the
$\dot {\theta }$
–
$\theta$
plane, with the baseline trajectory in blue and the perturbed trajectory in red. For the stable case, the perturbed trajectory initially deviates but gradually converges back to the baseline limit cycle, demonstrating recovery to the original periodic motion. In contrast, the flutter case exhibits sustained divergence and expanding amplitude, indicating exponential growth and onset of instability. Overlaid instantaneous vorticity fields at representative instants capture the evolution of the flow structures. In the baseline and stable cases, the flow remains largely periodic with limited fluctuation in vortex strength. However, in the flutter case, the larger structural amplitude leads to intensified vortex shedding.
Figure 4(d–f) presents the corresponding time histories of pitch angle
$\theta (t)$
and aerodynamic moment coefficient
$C_M(t)$
. The stable case reveals bounded oscillations that settle after the transient, whereas the flutter case displays an exponential increase in both
$\theta$
and
$C_M$
, consistent with the positive growth rate
$\lambda$
observed in figure 3(f).
The bottom two rows present time–frequency spectrograms of
$\theta (t)$
(figure 4
g–i) and
$C_M(t)$
(figure 4
j–l), obtained using continuous wavelet transforms to capture both temporal and spectral dynamics. All three cases exhibit a consistent dominance of structural natural frequency
$f_s$
in
$\theta$
, validating the robustness of the structural response. However, the spectrograms of
$C_M$
provide more nuanced insights into the flow–structure interaction. For the baseline case,
$C_M$
exhibits a persistent high-frequency content around
$f_{C_M} \approx 0.67$
, associated with vortex-shedding dynamics. In the stable perturbed case, a transient amplification of low-frequency components at
$f_{C_M} \approx 0.19$
(matching
$f_s$
) and its harmonic at
$f_{C_M} \approx 0.38$
is observed, but these decay over time. In contrast, the flutter case shows a continuous and growing amplification of these low-frequency components post-perturbation, indicating sustained energy transfer from the flow to the structure. This progressive emergence and amplification of the dominant instability
$f_{C_M} \approx 0.38$
and its subharmonic
$f_{C_M} \approx 0.19$
(matching
$f_s$
) in
$C_M$
reflect a phase lock-on phenomenon between the aerodynamic forcing and the structural motion. This lock-on effectively channels energy into the structural mode, resulting in flutter.
To further examine the underlying flow structures contributing to stability and instability, we apply DMD (Schmid Reference Schmid2022) to the unsteady flow fields following impulsive stiffness perturbations. Specifically, we compare the flutter-prone unstable case (
$\phi = 0.353\pi$
) and the stable case (
$\phi = 0.29\pi$
) discussed in figure 4. Dynamic mode decomposition is performed in a body-fixed frame to isolate flow oscillations associated with the airfoil’s pitching motion. The exact DMD analysis uses 1872 snapshots collected over 18 oscillation cycles with a time interval between snapshots of
$\Delta t/T = 0.0095$
.
Dynamic mode decomposition approximates the system’s evolution using a linear operator
$A$
such that
${\boldsymbol{\omega}}_{k+1} \approx A {\boldsymbol{\omega}}_k$
, where
${\boldsymbol{\omega}}_k$
denotes the vorticity field at time
$t_k$
. The discrete-time eigenvalues
$\lambda _j$
of
$A$
encode the growth rate and frequency of the corresponding DMD modes
$\varPhi _j$
. The continuous time spectra are obtained as
$\varLambda _j = {\log (\lambda _j)}/{\Delta t}$
, where
${\textrm{Re}} (\varLambda _j)$
gives the modal growth rate and
$f_\varPhi = {\textrm{Im}} (\varLambda _j)/2\pi$
gives the oscillation frequency. The mode amplitudes
$b_j$
are computed by projecting the initial state
${\boldsymbol{\omega}}_0$
onto the DMD modes via
$\boldsymbol{b} = \varPhi ^\dagger \boldsymbol{x}_0$
, where
$\dagger$
denotes the Moore–Penrose pseudoinverse.
Figure 5. The DMD of the flow field for stiffness perturbations at (a–c)
$\phi = 0.353\pi$
(flutter-prone) and (d–f)
$\phi = 0.29\pi$
(stable). (a,d) The DMD eigenvalues plotted in terms of growth rate versus frequency, with colours indicating mode amplitudes
$|b|$
. (b,e) The DMD spatial mode at
$f_\varPhi \approx 0.19$
(structural frequency). (c,f) The DMD spatial mode at
$f_\varPhi \approx 0.38$
(first harmonic) for flutter-prone case and spatial mode at
$f \approx 0.67$
for stable case. The DMD is performed in the laboratory frame aligned with airfoil pitching motion. In (a,d), the red circle indicates the DMD mean mode.
Figure 5(a–c) shows results for the flutter-prone unstable case. Figure 5(a) presents the DMD eigenvalues projected in the complex plane, with the growth rate
${\textrm{Re}} (\varLambda _j)$
on the
$x$
axis and frequency
$f_\varPhi$
on the
$y$
axis. The markers are coloured by modal amplitude
$|b_j|$
. The dominant modes appear at
$f_\varPhi \approx 0.38$
and its subharmonic
$f_\varPhi \approx 0.19$
. The subharmonic mode with
$f_\varPhi \approx 0.19$
has a positive growth rate, which confirms it as the driver of flutter. The corresponding spatial structures, shown in figures 5(b) and 5(c), exhibit strong vortex shedding from the trailing edge and large-scale separation from the suction surface. Figure 5(d–f) presents results for the stable case. In figure 5(d), all eigenvalues lie in the left-half complex plane, confirming that the dynamics is stable. The dominant modes, shown in figures 5(e) and 5(f), correspond to frequencies
$f_\varPhi \approx 0.19$
and
$f_\varPhi \approx 0.67$
(associated with moment coefficient oscillations). Repeating the analysis with a larger snapshot interval (
$\Delta t/T = 0.019$
) yielded similar dominant modes, confirming insensitivity to chosen parameters for DMD analysis. These findings support the conclusion that flutter is driven by phase lock-on between structural and fluid modes, resulting in selective energy amplification of low-frequency DMD modes. In the stable case, this coupling remains limited, preventing instability.
To suppress the flutter-prone instability triggered by negative stiffness perturbation at
$ \phi = 0.353\pi$
, we introduce a PBC strategy through transient heaving motion of the airfoil. Although the instability is induced via stiffness perturbations, stiffness itself is not easily actuated in practice. Heaving motion, on the other hand, can be implemented more readily through control surfaces, making it a physically realisable control input.
Following the PBC framework described in § 2, we compute the phase-sensitivity function
$ Z_h(\phi )$
and its gradient
$ {\textrm{d}}Z_h/{\textrm{d}}\phi $
by applying Gaussian-shaped (downward) heaving impulses at different phases of the oscillation cycle with impulse strength
$ \beta = -0.01$
. The resulting sensitivity profile, shown in figure 6(a), reveals alternating regions of phase advance and delay, identifying windows of high control effectiveness.
With this phase sensitivity, we formulate a control problem to compute an energy-optimal input
$ u(t)$
that induces a prescribed phase shift
$ \Delta \phi$
over a finite time horizon. The underlying rationale is that by advancing the phase trajectory to reach the same oscillatory position slightly ahead of its natural timing, we disrupt the temporal coherence between structural motion and unsteady aerodynamic forcing. This intentional desynchronisation weakens the feedback loop responsible for energy transfer and flutter amplification. The control input is obtained by minimising the cost function
\begin{equation} \mathcal{C} = \min _{u(t)} \quad \int _0^{T^*} \left [ \frac {1}{2} u^2(t) + \lambda (t) \left ( \dot {\phi }(t) - 2\pi f_s - Z_h(\phi ) u(t) \right ) \right ] {\textrm{d}}t, \end{equation}
where
$ \lambda (t)$
is the Lagrange multiplier enforcing the phase evolution constraint and
$ T^* = 0.9T$
is the reduced target period. The resulting optimal control signal
$ u(t)$
, shown in figure 6(c), is localised in time and requires minimal energy. While the present approach targets phase manipulation, in theory, the same framework could be extended to suppress specific unstable modes by directly shaping the modal energy growth through adjoint-informed or data-driven controllers. To evaluate robustness of the control design, we also introduce additive Gaussian noise to the phase-sensitivity function and its gradient and recompute the PBC signal. The resulting control waveform under a 20 dB signal-to-noise ratio (SNR) is shown as a dashed grey line in figure 6(c). Despite the noisy sensitivity function and its gradient, the control retains its overall structure, indicating robustness of the controller design to moderate levels of modelling noise.
Figure 6. The PBC of flutter instability via heaving actuation. (a) Phase-sensitivity function and its gradient for heaving input; negligible deviations at reduced impulse strength of
$\beta = -0.005$
are shown with error bars for characteristic phases. (b) Schematic of the energy-optimal phase-control framework targeting a prescribed phase shift. (c) Control inputs: optimal phase-based actuation (black), low-frequency sinusoid (
$f_1 = f_s/0.9$
) and high-frequency sinusoid (
$f_2 = f_s/0.18$
). (d) Resulting pitch angle
$\theta (t)$
and moment coefficient
$C_M(t)$
. (e) Spectrograms of
$C_M(t)$
for each actuation. (f) The DMD eigenvalues and dominant modes for the PBC case. Also shown is the Arnold tongue lock-on regime (highlighted in pink) corresponding to heave phase-sensitivity function in the inset of (d).
For comparison, we also consider two sinusoidal reference actuations with identical input power as PBC: one with frequency
$f_1 = f_s / 0.9$
(shown in red) and the other with
$f_2 = f_s/0.18 = 1.05$
(shown in blue). Here,
$f_2$
is obtained by selecting the dominant frequency of PBC. Interestingly, the dominant frequency of PBC lies in the high-frequency regime. This aligns with the asymmetry observed in the Arnold tongue computed for heave phase-sensitivity function in figure 6(d) (inset), where stronger lock-in behaviour occurs at frequencies higher than the structural natural frequency. The heaving actuation is turned off after
$T^*$
, allowing the system to evolve freely. We note that negative control input refers to downward heave motion of the airfoil. The corresponding pitch angle and moment coefficient responses are shown in figure 6(d). Only the PBC successfully stabilises the dynamics and suppresses flutter growth. In contrast, the low-frequency sinusoid accelerates divergence, while the high-frequency sinusoid delays instability but fails to prevent it. This comparison highlights the advantage of PBC in generating waveform-timed actuation aligned with the system’s intrinsic dynamics. However, its effectiveness in practice depends on precise timing. Delays in phase estimation or actuation due to bandwidth limits, mechanical lag or sensor inaccuracies can degrade performance. These constraints underscore the need for delay-compensated control schemes to account for actuator and sensing lags in experimental settings.
Spectrograms of
$C_M(t)$
(figure 6
e) reveal that only the phase-based controller suppresses the subharmonic growth at
$f \approx 0.19$
, while the sinusoidal inputs result in lock-on to the structural frequency. This underscores the importance of actuation timing and supports the role of phase-specific control in interrupting energy transfer to unstable modes. The DMD analysis of the vorticity field under PBC (figure 6
f) shows all modal eigenvalues (with the dominant spatial modes shown below) in the left-half complex plane, indicating fully stabilised dynamics.
4. Conclusion
This study demonstrates a phase-based framework for characterising and controlling flutter instabilities in a freely pitching airfoil. By applying localised stiffness perturbations across oscillation phases, we uncover critical phase regions that trigger instability, supported by transient energy growth, exponential amplification and spectral signatures. Dynamic mode decomposition further reveals low-frequency mode lock-on as the underlying mechanism of flutter onset. Building on these insights, we design an energy-optimal phase-control strategy using transient heaving motion. The proposed controller successfully suppresses flutter by preventing phase lock-on and inhibiting subharmonic amplification. Our findings highlight the critical role of phase in nonlinear aeroelastic interactions and for controlling flow-induced instabilities. This phase-based framework can be extended to investigate a wide range of flow-induced phenomena across large parameter spaces, including gust-driven and transonic flutter and aeroacoustic instabilities. At higher Reynolds numbers, turbulent flows exhibit greater modal complexity and diminished global coherence, making PBC more challenging. Nonetheless, recent work by Godavarthi et al. (Reference Godavarthi, Kawamura, Ukeiley, Cattafesta III and Taira2025) demonstrates that the phase dynamics of dominant coherent structures can still be extracted using modal decomposition techniques combined with phase reduction. This indicates that phase-based strategies may remain effective when applied to these structures through a reduced-order control framework.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10721.
Funding
A.G.N. acknowledges support from the Air Force Office of Scientific Research (Award No. FA9550-23-1-0483; Program Manager: Dr Gregg Abate).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data and code that support the findings are openly available on https://github.com/chathu14/Phase-based-control-for-low-Re-number-flows.
Open access
