1. Introduction
Stall flutter is an aeroelastic instability phenomenon closely related to flow separation, which manifests itself as a large-amplitude, single degree-of-freedom limit-cycle oscillation of an airfoil (Poirel, Goyaniuk & Benaissa Reference Poirel, Goyaniuk and Benaissa2018). This instability is commonly observed in various aeroelastic systems such as rotor blades, high-manoeuvrability flight vehicles and low-speed aircraft (Dos et al. Reference Dos Santos and Marques2021). Halfman, Johnson & Haley (Reference Halfman, Johnson and Haley1951) found that the occurrence of stall flutter is accompanied by a significant decrease in the critical flutter velocity of the airfoil, which endangers flight safety. Due to its highly destructive nature and strong coupling between structural and aerodynamic responses, stall flutter poses a serious threat to both performance and structural integrity. This motivates the need for effective control strategies to mitigate the stall flutter.
To establish the physical basis of the proposed control strategy, it is essential to understand the flow evolution associated with stall flutter. A crucial aerodynamic feature of stall flutter is dynamic stall, characterised by the eruption of vorticity near the leading edge that rolls up into a large-scale leading-edge vortex (LEV) and induces massive excursions in aerodynamic loads, which can couple with the structural dynamics of the airfoil (Corke & Thomas Reference Corke and Thomas2015). As long as the LEV remains attached, a delayed stall effect enhances lift (McCroskey Reference McCroskey1982); however, its downstream convection and shedding cause a sharp loss of lift and a rapid change in aerodynamic moment, which may result in negative damping and trigger stall flutter (Li & Feng Reference Li and Feng2022). Therefore, the prediction of the vortex behaviour in terms of shedding instant and peak circulation plays an important role in developing flow control strategies.
Two primary mechanisms of LEV shedding have been identified. The first is the chord-limiting mechanism, sometimes referred to as bluff-body shedding (Rival et al. Reference Rival, Kriegseis, Schaub, Widmann and Tropea2014; Widmann & Tropea Reference Widmann and Tropea2015). In this scenario, the LEV grows only up to the limit imposed by the finite chord length and the trailing-edge vortex (TEV) plays a crucial role in the LEV shedding (Eldredge & Jones Reference Eldredge and Jones2019). Widmann & Tropea (Reference Widmann and Tropea2015) showed that, when circulation of a LEV is normalised, considering the effect of an effective angle of attack, the shedding circulation associated with this mechanism collapses to a nearly uniform value. The second mechanism is boundary-layer eruption, where shedding is triggered not by a TEV but by the action of a secondary vortex (sV). Doligalski, Smith & Walker (Reference Doligalski, Smith and Walker1994) explained this as a viscous boundary-layer response to the adverse pressure gradient induced by the LEV itself, leading to flow separation once a critical threshold is exceeded. They proposed a fractional convection rate k
*
$=$
Γ/(4
$\pi$
yU∞), where Γ is the vortex circulation, y is the wall-normal distance and U∞ is the free-stream velocity, suggesting that k
* > 0.25 marks separation onset. However, this criterion is not directly applicable to all situations, and subsequent research has attempted to refine this picture. For instance, Li & Feng (Reference Li and Feng2022) reported the modulated location of peak pressure as an indicator to predict the LEV shedding. Li et al. (Reference Li, Feng, Kissing, Tropea and Wang2020) reported variations in sV structures associated with boundary-layer eruption and developed a simplified model for LEV circulation that incorporates both an effective inflow velocity and effective angle of attack. Besides, other works noted that actual shedding events can involve features of both mechanisms (Huang & Green Reference Huang and Green2015).
In the present study we find that large-amplitude stall flutter corresponds to a boundary-layer-eruption-induced LEV-shedding topology, accompanied by the emergence of a second leading-edge vortex (sLEV). For this regime, we introduce a new scaling of the shedding circulation that characterises the boundary-layer-eruption mechanism, thereby providing a refined framework for interpreting and controlling stall flutter.
With this flow-topology-based understanding, control strategies can be more precisely tailored. Conventional approaches to stall flutter mitigation generally aim to weaken, delay or suppress the LEV or TEV (Chen et al. Reference Chen, Shi, Chen and Yao2022; Wu et al. Reference Wu, Dai, Yang, Hu and Huang2023; Chen, Wu & Chen Reference Chen, Wu and Chen2025), thereby alleviating the associated moment fluctuations. In contrast, the present work identifies a counter-intuitive mechanism: for large-amplitude pitching oscillations governed by the boundary-layer-eruption topology, stall flutter can be effectively suppressed by enhancing the sLEV. To the authors’ knowledge, the role of the sLEV in stall flutter suppression has not been reported in prior studies. This differentiates our approach from earlier efforts, such as Ju et al. (Reference Ju, Yu, Wang and Meng2025), where plasma-based actuation was used primarily to suppress the formation of a large-scale LEV during the upstroke. They observed that once the LEV was fully developed, the effectiveness of their method diminished because momentum mixing could no longer prevent trailing-edge separation during the downstroke. By contrast, our strategy exploits the post-LEV dynamics: amplifying the sLEV and strengthens its interaction with the TEV, which in turn disrupts the feedback loop that sustains stall flutter.
To enable precisely targeted vortex control, it is essential to analyse the mechanisms by which unsteady aerodynamic loads are generated. Quantitatively attributing the aerodynamic force or moment to individual flow structures is a non-trivial task, and several theories have been developed to address this challenge (Wu, Liu & Liu Reference Wu, Liu and Liu2018). Building on Howe’s integral force formulation (Howe Reference Howe1995), Li & Wu (Reference Li and Wu2018) proposed the vortex force map method to quantitatively extract body forces from flow fields. Gao & Wu (Reference Gao and Wu2019) introduced a Galilean-invariant force decomposition, ensuring that the decomposed force components remain unaffected by translational motion of the reference frame. Gao, Xie & Lu (Reference Gao, Xie and Lu2025) recently established a weighted integral framework, significantly extending the applicability of force decomposition and enabling its application to surface-stress distribution. Menon & Mittal (Reference Menon and Mittal2021a ) extended the force partitioning method to moments, establishing a general force and moment partitioning framework applicable to unsteady aerodynamic flows. In the present study, we utilise a clustering-based vortex-induced load partitioning method (VILPM) to reveal a previously unquantified chain mechanism: sLEV–TEV interaction pathway.
In addition to quantitative load partitioning tools, the realisation of precisely targeted vortex control requires flow control approaches capable of directly manipulating boundary-layer and vortex dynamics. Among the available strategies, active surface morphing has emerged as a particularly promising technique. By introducing localised surface deformations, it can delay flow separation, enhance lift and reduce drag under both steady and unsteady aerodynamic conditions (Huang et al. Reference Huang, Dai, Yang and Xia2022a ; He et al. Reference He, Guo, Xu, Feng and Wang2023). Building on this foundation, our previous work (Xia et al. Reference Xia, Dai, Huang and Yang2024) demonstrated that surface morphing can also mitigate stall flutter by reducing the net energy transfer from the unsteady flow to the structure by as much as 46 %. Since this energy transfer is the primary mechanism sustaining limit-cycle oscillations (LCO) as explained by Goyaniuk, Poirel & Benaissa (Reference Goyaniuk, Poirel and Benaissa2020), altering it through controlled morphing offers a viable pathway for stall flutter suppression.
Active surface morphing can be realised through several mature actuation technologies, including shape memory alloy (SMA) and piezoelectric macro fibre composite (MFC) actuators. Shape memory alloy-based morphing offers advantages such as high actuation force and structural simplicity, and has been widely explored in bionic wing concepts (Jia et al. Reference Jia, Chang, Zhou, Zhang, Wang and Zhang2025). Dong, Boming & Jun (Reference Dong, Boming and Jun2008) and Georges et al. (Reference Georges, Brailovski, Morellon, Coutu and Terriault2009) installed SMA springs between the wing skin and the supporting wing box to achieve controllable surface morphing. Piezoelectric MFC actuators are characterised by their lightweight design, flexibility and ability to provide high-frequency actuation with low energy consumption (Wilkie et al. Reference Wilkie, Bryant, High, Fox, Hellbaum, Jalink, Little and Mirick2000). Using MFC actuators, Debiasi, Chan & Jadhav (Reference Debiasi, Chan and Jadhav2016) fabricated a symmetric wing with actively morphed surface and experimentally verified structural feasibility. More recently, He et al. (Reference He, Guo, Xu, Feng and Wang2023) showed that local surface morphing driven by MFC actuators can directly modify instantaneous vortex structures, thereby improving aerodynamic performance, highlighting the suitability of MFC-based morphing for flow control applications.
Effective flow control requires precise timing of actuation. As shown in Xia et al. (Reference Xia, Dai, Huang and Yang2024), the phase difference between the airfoil’s pitching motion and surface morphing largely determines control performance. In practice, however, once control is applied the system develops a time-varying response frequency and phase, which makes a fixed actuation phase ineffective. Even time-delayed feedback strategies (Wu, Dai & Yang Reference Wu, Dai and Yang2022) often fail to deliver stable performance because phase drift is inherent in nonlinear, unsteady aeroelastic systems. These challenges motivate the use of phase-locked-loop (PLL) control. Phase-locked-loop control (Hsieh et al. Reference Hsieh and Hung1996) has been widely adopted in oscillatory systems to track and regulate nonlinear dynamics (Peter & Leine Reference Peter and Leine2017; Abeloos et al. Reference Abeloos, Müller, Ferhatoglu, Scheel, Collette, Kerschen, Brake, Tiso, Renson and Krack2022). It maintains a specified phase lag between system response and actuation through a feedback loop, without requiring detailed system modelling (Denis et al. Reference Denis, Jossic, Giraud-Audine, Chomette, Renault and Thomas2018). For stall flutter suppression via active surface morphing, PLL offers a practical way to sustain a consistent phase lag between morphing deformation and airfoil motion, even when frequency and amplitude vary in time. The key is real-time estimation of instantaneous phase (Zhou et al. Reference Zhou, Raze, Kosova and Kerschen2025), which provides the input for continuously adjusting the morphing signal to remain locked to the target phase. Accurate phase detection is therefore central to PLL control. Two common approaches are the short-time Fourier transform (Griffin & Lim Reference Griffin and Lim1984) and variational mode decomposition (Chen et al. Reference Chen, Yang, Peng, Dong, Zhang and Meng2019). Short-time Fourier transform is limited in stall flutter applications by the low-frequency nature of the response, where frequency resolution is comparable to the signal itself, making transient tracking unreliable. Variational mode decomposition requires the amplitude’s time derivative to be much smaller than the response frequency, a condition that breaks down under control scenarios with strong amplitude attenuation. As a result, neither method is suitable for stall flutter control based on optimal phase morphing.
A further challenge in implementing PLL control lies in determining the optimal actuation phase. Conventional approaches rely on exhaustive parametric sweeps (Hu et al. Reference Hu, Dai, Wu and Yang2022; Xia et al. Reference Xia, Dai, Huang and Yang2024), which are computationally expensive and offer limited predictive value for new configurations. The underlying principle of active surface morphing is to locally perturb the boundary layer and thereby alter the dynamics of vortex shedding associated with stall flutter (Jaiman, Li & Chizfahm Reference Jaiman, Li and Chizfahm2023). Recent studies suggest that optimal control parameters in such unsteady flow control may correlate with certain baseline flow features, such as vortex shedding characteristics (Hirato et al. Reference Hirato, Shen, Gopalarathnam and Edwards2021; Li & Feng Reference Li and Feng2022; Zhang et al. Reference Zhang, Calderon-Sanchez, Almeida-Medina, Medina-Manuel and Souto-Iglesias2025) or energy transfer metrics (Huang, Dai & Yang Reference Huang, Dai and Yang2022b ). Motivated by this, we propose a flow-feature-based indicator that links the pre-control vortex dynamics directly to the optimal morphing phase, eliminating the need for costly trial-and-error searches. This indicator is formulated through a framework that combines VILPM with an energy-based map (EBM). While previous studies have highlighted the potential of load partitioning techniques in flow control design (Wu et al. Reference Wu, Liu and Liu2018), to the authors’ knowledge, the present work is the first to explicitly employ such methods to determine optimal control parameters.
Based on the above considerations, this study develops a PLL control framework for stall flutter suppression via phase-optimised surface morphing. Two distinct dynamic stall topologies governing stall flutter are first identified, and for the flow topology of boundary-layer-eruption-induced shedding at large pitching amplitudes, a new scaling of the shedding circulation is proposed. Within this topology, the sLEV emerges as a critical flow feature. A counter-intuitive control mechanism is then uncovered: by enhancing the sLEV rather than suppressing dynamic stall vortices, stall flutter is mitigated through the sLEV–TEV interaction. By combining VILPM with an EBM, a flow-feature-based indicator is further established that quantitatively links pre-control vortex dynamics to the optimal morphing phase, marking the first use of load partitioning techniques to guide control design. Finally, the PLL framework demonstrates complete suppression of stall flutter in the large-amplitude pitching regime dominated by boundary-layer-eruption topology, thereby delineating the effective operational boundary of sLEV-based control.
The remainder of the paper is organised as follows. Section 2 introduces the numerical set-up, the construction of the EBM, the clustering-based VILPM and the PLL control framework. Section 3 presents the flow topologies and circulation scaling, the sLEV–TEV control mechanism, the optimal phase migration mechanism and its critical indicator, followed by the performance of the PLL strategy. The effect of morphing amplitude is also examined, highlighting the robustness and generality of the proposed approach. Section 4 concludes with the main findings.
2. Computational model and numerical methods
2.1. Problem set-up
The two-dimensional computational model studied in this work employs a NACA0012 airfoil in an incompressible fluid with a low Reynolds number Re
$=$
U
∞
c/ν
$=$
1000, where c is the chord, U
∞
is the free-stream velocity and ν is the kinematic coefficient of viscosity, as shown in figure 1. The unsteady, incompressible Navier–Stokes equation is taken as the governing equations,
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}&=0 ,\end{align}
\begin{align} \frac{\partial \boldsymbol{u}}{\partial t}+\left(\left(\boldsymbol{u}-\boldsymbol{u}_{g}\right)\boldsymbol{\cdot }\boldsymbol{\nabla }\right)\boldsymbol{u}&=-\frac{1}{\rho }\boldsymbol{\nabla }\!p+\nu {\nabla} ^{2}\boldsymbol{u} ,\end{align}
where ρ, u and p denote the density, velocity and pressure of the fluid, respectively, and u g denotes the velocity of the grid.
Schematic of pitch elastically mounted airfoil with active surface morphing.
Figure 1. Long description
The diagram illustrates a pitch elastically mounted airfoil with active surface morphing. It features an elastic axis, a maximum morphing shape depicted by a blue dashed line, and an active morphing region indicated by a blue solid line. The airfoil is shown with two key points labeled as xm,1 and xm,2, and an angle alpha between the elastic axis and the chord line. The diagram also includes a representation of the critical flutter velocity U infinity, a damper labeled D, and a spring labeled K. The active morphing region is highlighted, and the overall structure demonstrates the interaction between aerodynamic and structural responses.
The torsional elasticity of the airfoil is modelled using a linear torsional spring with the pitch stiffness K and the pitch damping D, attached to the elastic axis X
ea
that is normalised by chord c and valued as X
ea
$=$
0.5. The equilibrium angular position of the airfoil is
$\alpha$
e
$=$
15°. The structural system of the airfoil is governed by the equation
\begin{equation} I_{\textit{ea}}\ddot{\alpha }+D\dot{\alpha }+K\left(\alpha -\alpha _{e}\right)=M_{\textit{ea}}\left(t\right)\!,\end{equation}
where
$\alpha$
is the angle of attack, t is the time, I
ea
is the moment of inertia about the elastic axis and M
ea
is the aerodynamic moment about the elastic axis. The directions of
$\alpha$
and M
ea
is defined as positive in the clockwise direction. Here
$C_{M}=M_{\textit{ea}}/(0.5\rho {U}_{\infty }^{2}c^{2})$
is the coefficient of aerodynamic moment. The dimensionless moment of inertia is
$I^{*}=I/(0.5\rho c^{4})=4.5$
. The damping ratio is
$\xi ^{*}=D/D_{cr}=0.03$
, where
$D_{cr}=2\sqrt{KI_{\textit{ea}}}$
is the critical damping for the harmonic oscillator. The spring stiffness is commonly expressed in terms of a reduced velocity
$U^{*}=U_{\infty }/f_{s}c=6$
, where
$f_{s}=(1/2\pi )\sqrt{K/I_{\textit{ea}}}$
is the natural frequency of the spring.
The unsteady incompressible Navier–Stokes equations with an arbitrary Lagrangian–Eulerian reference frame are solved by the finite volume method. The velocity–pressure coupling is solved with the pressure-implicit with splitting of operators (PISO) algorithm. A second-order backward implicit scheme is selected for time discretisation. The second-order scheme based on a linear upwind interpolation is used for the convection term. The diffusion term is discretised with a second-order accurate scheme based on linear interpolation, using an explicit non-orthogonal-corrected surface normal gradient scheme. Therefore, the numerical discretisation scheme has a second-order accuracy in temporal and spatial domains.
Computational domain and mesh details.
Figure 2. Long description
The diagram illustrates a computational domain with dimensions of 40c by 40c. It features a fixed mesh in the outer region and a moving mesh within a circular boundary of radius 14c. A non-conformal sliding interface separates the two meshes. A zoomed-in view highlights a smaller rectangular area within the moving mesh, showing detailed contour lines and an arrow indicating direction.
A coupled dynamic mesh technique, integrating the arbitrary mesh interface (AMI) method and the Laplacian smoothing method (Salehi & Nilsson Reference Salehi and Nilsson2023) has been employed to handle mesh updates induced by rigid-body pitching and surface morphing, respectively. As shown in figure 2, the computational domain is split into two subdomains: a circular moving mesh zone and a fixed mesh zone. The centre of a moving mesh zone is the elastic axis position of the pitching airfoil. The moving mesh zone can rotate with the airfoil while the fixed mesh zone remains stationary. The diameter of the moving mesh zone is 14c. The fixed mesh zone is a square with its side length as 40c. The moving mesh zone is linked together with the fixed mesh zone using the AMI technique through the non-conformal sliding interface. The free-stream velocity and pressure are given on the inlet and outlet boundaries, respectively. On the airfoil surface, the no-slip condition is applied, requiring that the fluid velocity at the airfoil surface matches the surface velocity. The code has been extensively validated in our previous studies to deal with the moving-boundary flow problems (Huang et al. Reference Huang, Dai, Yang, Wu and Xia2021, Reference Huang, Dai and Yang2022b
). The structural governing equation (2.3) is solved analytically at each time step. Taking the angle of attack
$\alpha ^{n}$
and pitch velocity
$\dot{\alpha }^{n}$
as an initial condition, the analytical solution of
$\alpha ^{n+1}$
and
$\dot{\alpha }^{n+1}$
at time step n + 1 can be calculated as
\begin{equation} \begin{array}{l} \alpha ^{n+1}={\rm e}^{\lambda {\rm d}t}\left(c_{1}\cos \mu {\rm d}t+c_{2}\sin \mu {\rm d}t\right)+{M}_{\textit{ea}}^{n+1}/K+\alpha _{e},\\ \dot{\alpha }^{n+1}={\rm e}^{\lambda {\rm d}t}\left[\left(c_{1}\lambda +c_{2}\mu \right)\cos \mu {\rm d}t+\left(c_{2}\lambda -c_{1}\mu \right)\sin \mu {\rm d}t\right] \!,\end{array} \end{equation}
where
\begin{align} \lambda & =-D/2I_{\textit{ea}},\nonumber\\ \mu & =\sqrt{4KI_{\textit{ea}}-D^{2}}/2I_{\textit{ea}},\nonumber\\ c_{1} & =\alpha ^{n}-{M}_{\textit{ea}}^{n+1}/K-\alpha _{e},\nonumber\\ c_{2} & =\left(\dot{\alpha }^{n}-\lambda c_{1}\right)/\mu ,\nonumber\\ {d}t & =t^{n+1}-t^{n}. \end{align}
The fluid-structure coupling employs a loosely coupled approach wherein the flow and structural equations are solved sequentially, and the fluid-structure interaction (FSI) code has been validated in our previous work (Huang et al. Reference Huang, Dai, Yang and Xia2022a ; Xia et al. Reference Xia, Dai, Huang and Yang2024). Further details about numerical validation can be found in Appendix A.
Uncontrolled FSI case (onset of stall flutter, baseline): (a) time-domain evolution of the airfoil pitching amplitude and dimensionless time
$t^* = t U_\infty / c$
; (b) C
M
hysteresis curve (black solid line) with flow field snapshots of six representative moments (t
1–t
6) coloured by spanwise vorticity contours.
Figure 3. Long description
The image contains two graphs. The first graph shows the time-domain evolution of the airfoil pitching amplitude in degrees against dimensionless time, tU/c. The pitching amplitude starts at zero, oscillates with increasing amplitude, and eventually stabilizes around 40 degrees. The second graph presents a hysteresis curve of the aerodynamic moment coefficient, CM, plotted against the angle of attack, alpha, in degrees. The black solid line represents the CM hysteresis curve. Six flow field snapshots at different moments (t1 to t6) are included, colored by spanwise vorticity contours. These snapshots illustrate the flow evolution during the stall flutter process. The snapshots are placed at specific points along the hysteresis curve, showing the dynamic stall and vortex shedding behavior. The flow field snapshots provide visual insights into the aerodynamic phenomena associated with stall flutter.
Figure 3 shows the development of the airfoil stall flutter under the uncontrolled FSI case, as the baseline of this work. The pitching amplitude grows rapidly from an initial state of rest and reaches a fully developed LCO before a non-dimensional time of
$t^*=tU_{\infty}/c=150$
, with a pitching amplitude of
$A_\alpha$
$=$
24.4°. The resulting stall flutter response exhibits a highly harmonic behaviour and can be well approximated by the expression
\begin{equation} \alpha =\alpha _{0}+A_{\alpha }\sin \left(2\pi f_{\alpha }t+\phi _{\alpha }\right)\!,\end{equation}
where
$\alpha$
0 is the mean angle of attack,
$f_{\alpha}$
and
$\phi_{\alpha}$
are the frequency and phase of the pitching airfoil.
To mitigate the stall flutter shown in figure 3, active surface morphing is applied to reduce the pitching amplitude A
$\alpha$
. The active morphing motion of the airfoil’s surface is considered as a harmonic motion of a single mode varying both in temporal and spatial domains, as shown in figure 1, which is given by the equation
\begin{align} y_{m}\left(t,\overline{x}\right) & =\frac{A_{m}}{2}\left[\sin \left(2\pi \overline{x}-\frac{\pi }{2}\right)+1\right]\sin \left(2\pi f_{m}t+\phi _{m}\right)\nonumber\\ & =A\left(\overline{x}\right)\sin \left(2\pi f_{\alpha }t+\phi _{\alpha }+{\unicode[Arial]{x0394}} \phi \right)\!, \end{align}
where y
m
is the surface morphing in the direction perpendicular to the chord and scaled by chord c; A
m
, f
m
and
$\phi_{m}$
are the amplitude, frequency and phase of the local surface oscillation, respectively, and f
m
$f_{m}=f_{\alpha}$
, Δ
$\phi$
is the phase difference between airfoil pitching and surface morphing. The coordinate
$\overline{x}=(x_{m}-x_{m,1})/(x_{m,2}-x_{m,1})$
, where x
m,1 and x
m,2 are the end points of the active morphing region, and scaled by chord c. In this study, x
m,1
$=$
0, x
m,2
$=$
0.6.
2.2. Energy-based map of stall flutter boundary
Previous studies have shown that the complex amplitude response of aeroelastic flutter can be systematically interpreted using energy-based approaches. Early works applied potential-flow or reduced-order models to estimate the stability boundary (Dowell Reference Dowell1966; Holmes & Marsden Reference Holmes and Marsden1978; Jones & Platzer Reference Jones and Platzer1996; Ramesh, Murua & Gopalarathnam Reference Ramesh, Murua and Gopalarathnam2015), while more recent efforts demonstrated that the energy extracted from the flow during one oscillation cycle provides a reliable indicator for predicting flutter bifurcations (Morse & Williamson Reference Morse and Williamson2009; Bhat & Govardhan Reference Bhat and Govardhan2013; Menon & Mittal Reference Menon and Mittal2019). The central idea is that the balance between the fluid energy input and the structural energy dissipation determines whether oscillations decay, persist as limit cycles or diverge.
We extend the energy-based idea to predict control performance efficiently and to track the evolution of the real-time aeroelastic response under control. For a pitching airfoil, the net energy transfer over one cycle is defined as
\begin{equation} E^{*}={E}_{f}^{*}-{E}_{d}^{*}, \end{equation}
where
$E_{f}^{\ast}$
is the energy extracted from the flow and E
d
* is the energy lost to structural damping,
\begin{align} {E}_{f}^{*}&={\int }_{t_{n}}^{t_{n}+T}C_{M}\dot{\alpha }{\rm d}t ,\end{align}
\begin{align} {E}_{d}^{*}&={\int }_{t_{n}}^{t_{n}+T}C_{\textit{damping}}\dot{\alpha }{\rm d}t={\int }_{t_{n}}^{t_{n}+T}\frac{D}{0.5\rho {U}_{\infty }^{2}c^{2}}\dot{\alpha }^{2}{\rm d}t ,\end{align}
where
$\dot{\alpha }$
is angular velocity and T is the oscillation period. The
$E^{\ast}$
can be used to analyse the changes in pitch oscillation amplitude. Equilibrium of the FSI system is attained when
$E^{\ast}$
$=$
0. Here
$E^{\ast}$
< 0 indicates the decay of oscillation amplitude, while
$E^{\ast}$
> 0 denotes the growth of oscillation amplitude. The stability of equilibrium is evaluated by
${\rm d}E^{*}/{\rm d}A_{\alpha }$
. The stable equilibrium is given by states with
$E^{\ast}$
$=$
0 and
${\rm d}E^{*}/{\rm d}A_{\alpha }\lt 0$
. More explanation about
$E^{\ast}$
can be found in Appendix B. For example, in the baseline uncontrolled case (figure 3), the airfoil develops a large-amplitude stable LCO, where the conditions
$E^{\ast}$
$=$
0 and
${\rm d}E^{*}/{\rm d}A_{\alpha }\lt 0$
are satisfied simultaneously. The associated vortex dynamics sustain the oscillation by replenishing energy extracted from the flow. Under control, the reduction of E
* below zero reflects the suppression of stall flutter, with the degree of mitigation increasing as
$E^{\ast}$
decreases.
For control analysis,
$E^{\ast}$
is mapped onto the parameter space spanned by the control parameter (morphing phase, here) and the response of the vibration system (pitching amplitude, here), as shown in figure 11. The contour
$E^{\ast}$
$=$
0 delineates the stall flutter boundary, enabling prediction of the controlled response without fully coupled simulations. The validity of the EBM is confirmed in § 3.4.
Compared with conventional aeroelastic simulations or experiments, the EBM offers distinct advantages: it allows full control over the prescribed kinematics, enables systematic parameter sweeps to delineate stall flutter boundaries and provides direct visual insight into LCO behaviour and control parameter dependence. Moreover, it captures subcritical instabilities that are often missed in flow-driven simulations, reduces the computational cost by avoiding long convergence times in time-resolved FSI computations and incorporates structural damping into the energy balance by (2.10) without added complexity (Turner, Seo & Mittal Reference Turner, Seo and Mittal2024). In this paper, it further serves as a foundation for constructing flow-feature-based indicators when combined with VILPM.
2.3. Clustering-based VILPM
In this study we integrate the load partitioning method with the adaptive variable density-based spatial clustering (AVDSC) algorithm developed by Xia et al. (Reference Xia, Dai, Huang, Yang and Ferrer2026) to systematically identify and evaluate the contributions of distinct vortical structures to the total aerodynamic load. A brief description of the combined methodology is presented below.
The key load we are concerned with in the stall flutter problem is the aerodynamic moment M ea on the airfoil with surface B, immersed in a flow domain V f bounded by surface Σ, consisting of a pressure component and friction component, i.e.
\begin{equation} M_{\textit{ea}}=\boldsymbol{M}\boldsymbol{\cdot }\boldsymbol{e}_{k}=\int _{B}p\left[\left(\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}}\right)\times \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{k}\right]{\rm d}S+\int _{B}\nu \rho \left[\left(\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}}\right)\times \left(\boldsymbol{n}\times \boldsymbol{\omega }\right)\boldsymbol{\cdot }\boldsymbol{e}_{k}\right]{\rm d}S ,\end{equation}
where p is pressure, ω is the vorticity vector, X is the position vector of a point on the airfoil surface, n is the normal unit vector of the airfoil surface and e k is the Cartesian basis vector in the k direction. Due to the direction of M ea is defined as positive in the clockwise direction, e k is set as (0, 0, –1) in this study.
The first step in the VILPM is the calculation of an auxiliary potential
$\psi_{k}$
, which is defined as
\begin{align} {\nabla} ^{2}\psi _{k}=0\text{ with }\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}=\begin{cases} \left(\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}}\right)\times \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{k} \quad\text{ on }B,\\ 0 \quad\text{ on }\sum. \end{cases} \end{align}
The next step is to project the Navier–Stokes equation onto the gradient of the auxiliary potential and integrate over the volume of the fluid domain. This results in
\begin{align} M_{\textit{ea}}&=M_{k}+M_{Q}+M_{vp}+M_{v}+M_{\sum }\nonumber\\ &=\underset{M_{k}}{\underbrace{-\rho \int _{B}\boldsymbol{n}\boldsymbol{\cdot }\left(\frac{{\rm d}\boldsymbol{u}}{{\rm d}t}\psi _{k}\right){\rm d}S}}+\underset{M_{Q}}{\underbrace{\rho \int _{V_{f}}2\left[\frac{1}{2}\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\right)\right]\psi _{k}{\rm d}V}}\nonumber\\&\quad +\underset{M_{vp}}{\underbrace{\nu \rho \int _{B}\left(\boldsymbol{\omega }\times \boldsymbol{n}\right)\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}S}}+\underset{M_{v}}{\underbrace{\nu \rho \int _{B}\left(\boldsymbol{n}\times \boldsymbol{\omega }\right)\boldsymbol{\cdot }\left[\boldsymbol{e}_{k}\times \left(\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}}\right)\right]{\rm d}S}}\nonumber\\&\quad \underset{M_{\sum }}{\underbrace{-\rho \int _{\sum }\boldsymbol{n}\boldsymbol{\cdot }\left(\frac{{\rm d}\boldsymbol{u}}{{\rm d}t}\psi _{k}\right){\rm d}S+\nu \rho \int _{\sum }\left(\boldsymbol{\omega }\times \boldsymbol{n}\right)\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}S}}, \end{align}
where the term
$M_{\Sigma}$
is negligible as demonstrated by Zhang (Reference Zhang2015) and
\begin{equation} M_{Q}=\int _{V_{f}}\left(-2\rho Q\psi _{k}\right){\rm d}V=\int _{V_{f}}m{\rm d}V .\end{equation}
More details of the derivation can be found in Appendix C. When the integrand m is integrated over the region occupied by an individual vortex, it gives the contribution of an individual vortex to the aerodynamic moment. This specific aspect is a key focus of the VILPM.
To calculate the moment contributions from distinct vortical structures, vortex detection is first performed using the Q-criterion (Jeong & Hussain Reference Jeong and Hussain1995). A threshold value Q th is applied to identify the dominant vortical regions (Q > Q th ), effectively filtering out weak vortices with negligible influence on the aerodynamic moment. These identified regions are then segmented into individual vortex structures using the AVDSC algorithm. A subsequent tracking procedure assigns consistent labels to each vortex across time steps, enabling temporal continuity. By isolating the corresponding spatial volumes of these structures throughout their evolution, the aerodynamic moment induced by each vortex can be quantitatively evaluated following (2.14).
2.4. Phase-locked-loop control framework
The effects of targeted vortex control are validated in a closed-loop control framework. According to Xia et al. (Reference Xia, Dai, Huang and Yang2024), the phase difference Δ
$\phi$
is a crucial factor in surface morphing control. Even slight variations in Δ
$\phi$
can significantly impact control effectiveness. Once control is applied, the fluid-structure system is a time-varying system, whose amplitude, frequency and phase will change with time, causing the phase of the pitching motion to shift over time gradually. To maintain Δ
$\phi$
at a specific value, such as the optimal phase difference, real-time estimation of the amplitude, frequency and phase of the pitching motion is required. After initiating control, the pitching angle signal at the current time step,
$\alpha$
n
, is collected and combined with the pitching angle signals from the preceding Ns time steps to form a time series, i.e.
\begin{equation} \alpha=\left[\begin{array}{llll} \alpha _{n-Ns} & \cdots & \alpha _{n-1} & \alpha _{n} \end{array}\right]^{T} \end{equation}
for each current time step, assuming that the corresponding pitching angle time series follows the time-domain response function given by
\begin{equation} \alpha \left(t\right)=\alpha _{0,n}+A_{\alpha ,n}\sin \left(2\pi f_{\alpha }t+\phi _{\alpha ,n}\right) .\end{equation}
According to the least squares method,
\begin{align} \left[\begin{array}{l} A_{\alpha ,n}\cos \phi _{\alpha ,n}\\ A_{\alpha ,n}\sin \phi _{\alpha ,n}\\ \alpha _{0,n} \end{array}\right]&=\left[\begin{array}{lll} \sin \left(2\pi f_{\alpha }t_{n-Ns}\right) & \cos \left(2\pi f_{\alpha }t_{n-Ns}\right) & 1\\ \vdots & \vdots & \vdots \\ \sin \left(2\pi f_{\alpha }t_{n}\right) & \cos \left(2\pi f_{\alpha }t_{n}\right) & 1 \end{array}\right]^{-1}\left[\begin{array}{l} \alpha _{n-Ns}\\ \vdots \\ \alpha _{n} \end{array}\right] \!,\end{align}
\begin{align} A_{\alpha ,n}&=\sqrt{\left(A_{\alpha ,n}\cos \phi _{\alpha ,n}\right)^{2}+\left(A_{\alpha ,n}\sin \phi _{\alpha ,n}\right)^{2}} ,\end{align}
\begin{align} \phi _{\alpha ,n}&=\arctan \left(A_{\alpha ,n}\cos \phi _{\alpha ,n}/A_{\alpha ,n}\sin \phi _{\alpha ,n}\right)\!,\end{align}
the amplitude
$A_{\alpha,n}$
, phase
$phi_{\alpha,n}$
and mean value
$\alpha_{0,n}$
of the pitching oscillation at the current time step are computed. These computed values are then appended to their respective time series
$A_{\alpha,n}$
,
$phi_{\alpha,n}$
and
$\alpha_{0,n}$
. In general, the instantaneous amplitude, frequency and phase is a smooth (or low-pass) function. To alleviate the influence of noise, a low-pass filtering function is applied to these signal sequences to obtain the filtered values given by
\begin{align} {\boldsymbol{A}}_{\alpha ,n}^{F}&=\boldsymbol{\varPhi }^{-1}\boldsymbol{A}_{\alpha ,n},\nonumber\\ {\boldsymbol{\phi }}_{\alpha ,n}^{F}&=\boldsymbol{\varPhi }^{-1}\boldsymbol{\phi }_{\alpha ,n},\nonumber\\ {\boldsymbol{\alpha }}_{\alpha ,n}^{F}&=\boldsymbol{\varPhi }^{-1}\boldsymbol{\alpha }_{\alpha ,n}, \end{align}
where
\begin{align} \boldsymbol{A}_{\alpha ,n} & =[\!\!\begin{array}{llll} A_{\alpha ,n-Nt+1} & \cdots & A_{\alpha ,n-1} & A_{\alpha ,n} \end{array}\!\!]^{T},\nonumber\\ \boldsymbol{\phi }_{\alpha ,n} & =[\!\!\begin{array}{llll} \phi _{\alpha ,n-Nt+1} & \cdots & \phi _{\alpha ,n-1} & \phi _{\alpha ,n} \end{array}\!\!]^{T},\nonumber\\ \boldsymbol{\alpha }_{\alpha ,n} & =[\!\!\begin{array}{llll} \alpha _{\alpha ,n-Nt+1} & \cdots & \alpha _{\alpha ,n-1} & \alpha _{\alpha ,n} \end{array}\!\!]^{T}, \end{align}
and Nt represents the length of the filtered signal sequence. The low-pass filtering function used in this study is defined as
\begin{equation} \boldsymbol{\varPhi }=\boldsymbol{I}+\frac{1}{\mu }\boldsymbol{D}^{\boldsymbol{T}}\boldsymbol{D} ,\end{equation}
with
\begin{equation} \boldsymbol{D}=\left[\begin{array}{lllllll} 1 & -2 & 1 & 0 & 0 & \cdots & 0\\ 0 & 1 & -2 & 1 & 0 & \cdots & 0\\ \vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & 1 & -2 & 1 & 0\\ 0 & \cdots & 0 & 0 & 1 & -2 & 1 \end{array}\right]\in \mathrm{\mathbb{R}}^{(Nt-2)\times Nt}, \end{equation}
where I is the identity matrix and μ is a regularisation parameter controlling the balance between signal fidelity and smoothness, in which smaller values of μ lead to stronger filtering. The filtered amplitude, phase and mean value at the current time step are given by
\begin{align} {A}_{\alpha ,n}^{F}&={\boldsymbol{A}}_{\alpha ,n}^{F}\left(Nt\right)\!,\nonumber\\ {\phi }_{\alpha ,n}^{F}&\boldsymbol{=}{\boldsymbol{\phi }}_{\alpha ,n}^{F}\left(Nt\right)\!,\nonumber\\ {\alpha }_{\alpha ,n}^{F}&={\boldsymbol{\alpha }}_{\alpha ,n}^{F}\left(Nt\right)\!,\end{align}
and the corresponding instantaneous frequency is expressed as
\begin{equation} {f}_{\alpha ,n}^{F}=f_{\alpha }+\frac{{\boldsymbol{\phi }}_{\alpha ,n}^{F}\left(Nt\right)-{\boldsymbol{\phi }}_{\alpha ,n}^{F}\left(Nt-1\right)}{{d}t} ,\end{equation}
where dt is the time step. At this point, the output of the surface morphing actuator is determined by
\begin{equation} y_{m,n}=A\left(\overline{x}\right)\sin \left(2\pi f_{\alpha }t_{n}+{\phi }_{\alpha ,n}^{F}+{\unicode[Arial]{x0394}} \phi \right) .\end{equation}
The time-domain fluid-structure-control (FSC) coupled simulation of the surface morphing airfoil plays a significant role in the design of the control law and the analysis of closed-loop dynamics. In this study, a time-domain FSC simulator based on a loosely coupled method (Zhang & Hisada Reference Zhang and Hisada2004) is employed, as shown in figure 4.
Phase-locked-loop control framework in the time-domain FSC simulator.
Figure 4. Long description
The flowchart begins with the start node and proceeds to the structural solver, which solves structural equations involving mass, damping, and stiffness matrices. The solver outputs displacement and velocity, which are then used in the flow solver. The flow solver updates mesh node positions induced by rigid-body pitching using the AMI method and solves CFD governing equations. The process then checks if the end time is reached. If not, it updates the time and loops back. The control system decides if control starts and whether to use a phase-locked loop or open-loop control. Depending on the decision, it obtains real-time estimation of the pitching state or uses the pitching state at the initial time. The control system then updates mesh node positions induced by surface morphing using the Laplacian smoothing method.
3. Results and discussions
3.1. Flow topologies and circulation scaling
To establish the physical basis of the proposed control strategy, it is necessary to first characterise the uncontrolled flow evolution that underlies stall flutter. As shown in § 2.1 (figure 3), the baseline FSI case develops a strongly harmonic pitching LCO. We therefore take this baseline case as a reference and examine its instantaneous flow topology and vortex dynamics. Figure 5 presents representative vorticity fields for the baseline; black contours denote vortex boundaries detected by the Q-criterion, which remains consistent with the vortex outlines obtained by the VILPM identification procedure described in § 2.3.
Flow field development in the uncontrolled FSI baseline (
$A_\alpha$
$=$
24.4°): (a) pitching-up phase; (b) pitching-down phase. Fields are coloured by spanwise vorticity and overlaid with VILPM-identified vortex boundaries (black contours).
Figure 5. Long description
The image shows the flow field development in the uncontrolled FSI baseline during pitching-up and pitching-down phases. Fields are colored by spanwise vorticity and overlaid with VILPM-identified vortex boundaries in black contours. The pitching-up phase includes six time points from -0.24T to 0.25T, showing the formation and shedding of the leading-edge vortex (LEV). The pitching-down phase includes six time points from 0.28T to 0.75T, showing the interaction between the LEV, trailing-edge vortex (TEV), and secondary vortex (sV).
Topologically, LEV formation and its subsequent shedding are central to the unsteady aerodynamics. As illustrated in figure 6(a), the nascent LEV is initially bounded by two half-saddles: one at the leading edge where the separated shear layer feeds the LEV and a second at the reattachment point on the suction surface. During the pitching-up phase the LEV grows and the reattachment point moves downstream along the airfoil surface. The subsequent topology of shedding depends on the behaviour of this reattachment motion and the generation of secondary vorticity. Two distinct LEV-shedding mechanisms have been identified in the literature, commonly termed boundary-layer eruption and bluff-body shedding (Rival et al. Reference Rival, Kriegseis, Schaub, Widmann and Tropea2014; Widmann & Tropea Reference Widmann and Tropea2015). A useful operational criterion is whether LEV shedding occurs before the reattachment point reaches the trailing edge or only after the reattachment has reached the trailing edge. In the present baseline the LEV sheds prior to reattachment reaching the trailing edge, placing this case in the boundary-layer-eruption category (hereafter topology I), as schematised in figure 6(b). We examined the flow topologies over a range of pitching amplitudes relevant to stall flutter suppression. Figure 7 displays the uncontrolled flow snapshots for different pitching amplitudes showing the integrand m of the vortex-induced moment defined in (2.14), overlaid with VILPM-identified vortex boundaries (black contours). For amplitudes from
$A_\alpha$
$=$
30° down to 10°, the shedding consistently belongs to the topology I regime, as exemplified by figures 7(b) and 7(c). By contrast, at sufficiently small pitching amplitudes the LEV grows until the reattachment reaches the trailing edge and then detaches in a chord-limited, bluff-body-like manner (topology II), as schematised in figure 6(c). An example of this behaviour is shown for
$A_\alpha$
$=$
5° in figure 7(a).
Schematic flow topologies associated with LEV evolution: (a) LEV formation; (b) topology I, boundary-layer-eruption-induced LEV shedding; (c) topology II, bluff-body-type shedding.
Figure 6. Long description
The diagram consists of three parts labeled (a), (b), and (c), each depicting different stages of LEV evolution. In part (a), the LEV formation is shown with an arrow indicating the direction of flow and the LEV labeled in red. Part (b) illustrates topology I, where the LEV and a secondary LEV (sLEV) are present, along with a secondary vortex (sV). Arrows indicate the flow direction and the shedding of the sLEV. Part (c) shows topology II, depicting the LEV and a trailing-edge vortex (TEV) with arrows indicating the flow direction and the shedding process. The diagram highlights the dynamic stall process and the interaction between the LEV and other vortices, which are crucial for understanding stall flutter and developing flow control strategies.
Uncontrolled flow snapshots for different pitching amplitudes: (a)
$A_\alpha$
$=$
5° (topology II); (b)
$A_\alpha$
$=$
10° and (c)
$A_\alpha$
$=$
20° (topology I).
Figure 7. Long description
A heat map illustrates the flow patterns for different pitching amplitudes. The map is divided into three sections labeled (a), (b), and (c), each representing different amplitudes: 5 degrees, 10 degrees, and 20 degrees respectively. Each section shows the flow patterns at specific times, with color gradients indicating the intensity of the flow. The labels LEV, sLEV, and cV mark specific vortex structures within the flow. The color scale ranges from blue to red, with blue indicating lower intensity and red indicating higher intensity. The flow patterns show distinct vortex formations and interactions, highlighting the changes in flow dynamics as the pitching amplitude increases.
To provide a quantitative corroboration, figure 8 shows the temporal evolution of LEV circulation for several
$A_\alpha$
and the corresponding TEV circulation evaluated at the instant of LEV shedding. Following the classical definition (Fage & Johansen Reference Fage and Johansen1928), we identify LEV shedding as the time when the circulation reaches a peak value, indicating that the vortical structure is no longer fed any more vorticity through the shear layer. At the LEV-shedding instant, the TEV is already strong for the
$A_\alpha$
$=$
5° case, whereas the TEV circulation is negligible for the
$A_\alpha$
$=$
10°–30° cases, which is insufficient to trigger LEV shedding. This distinction reinforces the topological classification above and suggests fundamentally different inter-vortex dynamics in the two regimes.
(a) Temporal development of LEV circulation for several A
$\alpha$
; (b) TEV circulation evaluated at the LEV-shedding instant indicated by red circles and blue circles mark the maximum TEV circulation for comparison.
Figure 8. Long description
The image contains two line graphs. The first graph on the left shows the temporal development of leading-edge vortex (LEV) circulation for various angles of attack (Aα) ranging from 5 degrees to 30 degrees. The x-axis represents time normalized by the period (Time/T), and the y-axis represents the magnitude of circulation (|Γ|) in square meters per second. Different colored lines represent different angles of attack, with labels indicating specific angles. Key points such as LEV shedding are marked with circles. The second graph on the right shows trailing-edge vortex (TEV) circulation evaluated at the LEV-shedding instant, indicated by red circles, and the maximum TEV circulation, marked by blue circles. The x-axis represents the angle of attack (Aα) in degrees, and the y-axis represents the magnitude of circulation (|Γ|) in square meters per second. The graph is divided into two topological regions, labeled Topology II and Topology I. The red circles indicate TEV circulation at the LEV shedding instant, while the blue circles mark the maximum TEV circulation for comparison. All values are approximated.
Given the severe aerodynamic excursions and serious structural failure associated with large-amplitude stall flutter, the remainder of this study focuses on the topology I regime. Within this topology, the LEV induces an adverse pressure gradient along the airfoil surface, which eventually drives boundary-layer separation once a critical threshold is reached. This eruption process ejects vorticity from the boundary layer, forming an sV accompanied by a sLEV (Foss Reference Foss2004). Vorticity from the leading-edge shear layer is then entrained into the sLEV, which in turn starts to grow. At the early stage after the formation of the sV and sLEV, the LEV region identified by the Q-criterion includes some connection region between the LEV and leading-edge boundary layer. The vorticity flux at the boundary of this connection region has a significant negative value. Note that the circulation of the LEV has a negative value too. This illustrates that the formation of the sV and sLEV does not cut off the vorticity transport from the leading-edge boundary layer into the LEV immediately, which results in the further increase of LEV circulation. After the LEV circulation reaches the peak value, the connection between the LEV and this connection region become weak rapidly, and this connection region is identified as an individual vortex independence on the LEV. At this moment, the vorticity transport from the leading-edge boundary layer into the LEV is completely cut off. Meanwhile, the separation of this connection region from the LEV results in the sharp decrease in LEV circulation immediately after LEV shedding.
To account for the kinematic dependence, a new circulation scaling is proposed. The reference velocity is taken as the chord-parallel component U
ref
$=$
U
∞cos
$\alpha$
and the reference length is taken as the chord projection associated with the maximum effective angle, L
ref
$=$
csin
$\alpha$
max. Peak circulation values are therefore non-dimensionalised as Γ/L
ref
U
ref
. Following Widmann & Tropea (Reference Widmann and Tropea2015), these values are further scaled by the dimensionless mean peak circulation Γ
ref
* computed over topology I cases. Figure 9 compares the conventional normalisation Γ/cU
∞ with the proposed normalisation Γ/ L
ref
U
ref
. For topology I (boundary-layer eruption), the new normalisation at the instant of shedding produces more uniform values, whereas the conventional normalisation exhibits strong deviation. For topology II (bluff-body shedding), the peak circulation is relatively small under both normalisations. These results indicate that, when LEV shedding is governed by boundary-layer eruption, the shedding circulation is principally determined by kinematic parameters that can be captured by the proposed scaling; conversely, simple chord-and-free-stream scaling fails to incorporate the role of pitching kinematics.
Comparison of LEV-shedding circulation normalisations.
Figure 9. Long description
A line graph compares LEV-shedding circulation normalizations. The x-axis represents the angle of attack (A_alpha) in degrees, ranging from 0 to 35 degrees. The y-axis represents the normalized circulation (Gamma/L_ref/U_ref) divided by Gamma_ref, ranging from 0 to 1.5. The graph is divided into two regions: Topology II (shaded in gray) and Topology I (unshaded). Two data series are plotted: one with open red circles and another with filled red circles. The open red circles represent data where U_ref equals U_infinity and L_ref equals c. The filled red circles represent data where U_ref equals U_infinity times the cosine of alpha and L_ref equals c times the sine of alpha_max. The data points show variations in normalized circulation across different angles of attack, with some points falling above or below the value of 1.0 on the y-axis. All values are approximated.
3.2. Vortex-interaction mechanism: sLEV–TEV pathway
The large-amplitude stall flutter of interest here corresponds to topology I, which is dominated by a large-scale LEV with a characteristic scale of approximately 0.5c. Achieving control using local surface morphing of amplitude an order of magnitude smaller (here A
m
$=$
0.05c) is therefore a non-trivial challenge. To evaluate the effectiveness of such a control approach, we first impose the same prescribed pitching motion as in the uncontrolled baseline (§ 2.1, figure 3) and apply active surface morphing with A
m
$=$
0.05c. Figure 10(a) shows the net energy gain
$E^{\ast}$
as a function of morphing phase difference Δ
$\phi$
. The minimum value occurs at Δ
$\phi$
$=$
0.5
$\pi$
, indicating this configuration yields the most effective mitigation under the given kinematics. The corresponding hysteresis loop of pitching moment C
M
is presented in figure 10(b). Compared with the baseline, the controlled case (red curve) exhibits two key modifications: (i) a reduction of the moment peak, and (ii) a substantial increase in the moment trough. Together, these changes demonstrate that active morphing can strongly alter the flow-induced moment generation process. Further details of the test procedure are reported in Xia et al. (Reference Xia, Dai, Huang and Yang2024).
(a) Net energy gain
$E^{\ast}$
value under active control with morphing amplitude A
m
$=$
0.05c and prescribed airfoil pitch motion with
$A_{\alpha}$
$=$
24.4°; (b) C
M
hysteresis curve for the active control case with Δ
$\phi$
$=$
0.5
$\pi$
(red solid line), and flow field snapshots of six representative moments (t
1–t
6) coloured by spanwise vorticity contours, compared with the baseline uncontrolled FSI case (black dotted line).
Figure 10. Long description
The image contains two graphs labeled (a) and (b) with flow field snapshots. Graph (a) is a line graph showing the net energy gain E* value under active control with morphing amplitude Am = 0.05c and prescribed airfoil pitch motion with A = 24.4 degrees. The x-axis represents the phase angle Δφ/π, and the y-axis represents the net energy gain E* = Ej* - Ed*. The graph shows a sinusoidal pattern with a peak at approximately 0.04 and a trough at approximately -0.08. Graph (b) is a line graph displaying the CM hysteresis curve for the active control case with Δφ = 0.5 (red solid line), compared with the baseline uncontrolled FSI case (black dotted line). The x-axis represents the angle of attack α in degrees, ranging from -10 to 40 degrees, and the y-axis represents the moment coefficient CM, ranging from -0.2 to 0.5. The graph shows the active control case following a different path compared to the uncontrolled case. Six flow field snapshots (t1 to t6) are included, colored by spanwise vorticity contours, showing different moments in the flow field.
The control performance prediction of EBM in the
$A_\alpha$
–Δ
$\phi$
parameter space for A
m
$=$
0.05c control, and the dashed line denotes the stall flutter boundary where
$E^{\ast}$
$=$
0.
Figure 11. Long description
A heat map displays the control performance prediction of EBM in the parameter space for Am equals 0.05c control. The x-axis represents the phase difference delta phi over pi, ranging from 0 to 2. The y-axis represents the angle of attack A alpha in degrees, ranging from 0 to 35. The color scale on the right indicates the value of E star, with blue representing lower values and red representing higher values. The dashed line denotes the stall flutter boundary where E star equals 0. A red arrow points downwards from a red dot labeled Prediction at delta phi equals 0.3 pi and A alpha equals 25 degrees. The heat map is divided into two regions labeled Topology I and Topology II.
The above analysis considers a fixed pitching amplitude (
$A_\alpha$
$=$
24.4°), corresponding to the baseline LCO. When control is applied to the fully coupled aeroelastic system, however, the pitching amplitude evolves dynamically. To predict this evolution, an EBM introduced in § 2.2 is employed. The EBM, constructed using bilinear interpolation of 216 discrete simulations over the
$A_\alpha$
– Δ
$\phi$
parameter space, is shown in figure 11. Blue regions indicate parameter combinations with a mitigating effect. At different values of
$A_\alpha$
, the control effectiveness of a given phase difference Δ
$\phi$
may vary. Nevertheless, certain phases consistently produce suppression. In particular, the EBM predicts that applying Δ
$\phi$
$=$
0.3
$\pi$
achieves complete suppression of stall flutter within topology I, as highlighted by the red trajectory.
These results raise a central question: How can a relatively small-scale surface morphing actuation (5 %c) disrupt the dynamics of a much larger-scale flow structure (LEV of order 50 %c) to mitigate and eventually suppress the aeroelastic oscillation? To address this, the VILPM introduced in § 2.3 is employed, revealing the sLEV–TEV interaction pathway, as follows.
To systematically investigate the effectiveness across the relevant
$A_\alpha$
range, prescribed pitching motions with amplitudes of 10°, 15°, 20° and 25° are analysed. Figure 12 presents the corresponding C
M
hysteresis loops, comparing the uncontrolled case (blue dotted lines) and the active control case with Δ
$\phi$
$=$
0.3
$\pi$
(red solid lines). For reference, results using the optimal phase determined for each
$A_\alpha$
based on energy map predictions are also included (yellow solid lines).
The
$C_M$
hysteresis curves as the pitching amplitude
$A_{\alpha}$
decreases from 25° to 10°, comparing the uncontrolled case (blue dotted line), the active control case with
$\Delta\phi$
$=$
0.3
$\pi$
(red solid line) and the active control case with the optimal phase for each
$A_\alpha$
(yellow solid line).
Figure 12. Long description
The image contains four graphs labeled (a), (b), (c), and (d), each depicting CM hysteresis curves as the pitching amplitude A decreases from 25 to 10 degrees. Each graph compares three cases: the uncontrolled case (blue dotted line), the active control case with a phase shift of 0.3 (red solid line), and the active control case with the optimal phase for each amplitude (yellow solid line). The x-axis represents the angle of attack (alpha) in degrees, ranging from -10 to 40 degrees, while the y-axis represents the moment coefficient (CM), ranging from -0.2 to 0.5. The graphs show how the CM hysteresis curves change with different pitching amplitudes and control strategies. The arrows indicate the direction of the hysteresis loops. The uncontrolled case shows larger hysteresis loops compared to the actively controlled cases, demonstrating the effectiveness of active control in reducing the hysteresis effect. The optimal phase control further minimizes the hysteresis loops, indicating the most efficient control strategy.
Across the range of
$A_\alpha$
where effective control is observed, the fundamental surface morphing mechanisms remain consistent and can be summarised as identified in figure 10: (i) a reduction in the peak aerodynamic moment, and (ii) a substantial enhancement of the moment trough. Between these two mechanisms, mechanism (ii) is found to be the dominant contributor to control effectiveness. The superiority of the optimal phase settings for each
$A_\alpha$
primarily arises from their greater influence on enhancing the trough compared with non-optimal phases. A detailed analysis of the optimal phase and its underlying mechanism will be provided in § 3.3.
Following the identification of overall control effects, the flow mechanisms under control at Δ
$\phi$
$=$
0.3
$\pi$
are further analysed using the VILPM. A pitching amplitude of
$A_\alpha$
$=$
20° is first selected for detailed study. Figure 13(a–c) shows the results for the uncontrolled case, while figure 13(d–f) presents the corresponding analysis for the active control case. Figure 13(a,d) compares the total aerodynamic moment and its decomposition into four physical contributions: kinematic, vortex-induced, viscous pressure and friction moments, calculated according to (2.13). It is evident that the vortex-induced moment dominates the total aerodynamic moment throughout the pitching cycle in both cases. To further dissect the vortex-induced contribution, the moment is separated into components arising from rotation-dominated (Q > 0) and strain-dominated (Q ≤ 0) flow regions. As shown in figure 13(b,e), while the strain-induced moment follows a similar temporal trend to the rotation-induced moment, its magnitude remains consistently lower across the cycle. Similar findings are also found by Menon & Mittal (Reference Menon and Mittal2021b
) and Xia et al. (Reference Xia, Dai, Huang, Yang and Ferrer2026). This shows that the key contributions are associated with vortical structures in rotation-dominated regions.
The VILPM analysis for the (a–c) uncontrolled case and (d–f) active control case with
$\Delta\phi$
$=$
0.3
$\pi$
under a prescribed pitching motion with
$A_{\alpha}$
$=$
20°. (a,d) Total aerodynamic moment and decomposition into kinematic, vortex-induced, viscous pressure and friction contributions. (b,e) Comparison of vortex-induced moment contributions from regions with Q > Q
th
, Q > 0 and Q ≤ 0. (c, f) Contributions of key vortex structures.
Figure 13. Long description
The image contains six graphs analyzing the VILPM for both uncontrolled and active control cases under a prescribed pitching motion with an amplitude of 20 degrees. The first row shows total aerodynamic moment and its decomposition into kinematic, vortex-induced, viscous pressure, and friction contributions for both cases. The second row compares vortex-induced moment contributions from regions with different Q values. The third row highlights the contributions of key vortex structures. Each graph uses a time over T axis and a C M axis to represent the data. The graphs illustrate how different factors contribute to the aerodynamic moment and how active control affects these contributions.
To focus on the most influential structures, weak vortices are filtered out by setting a threshold Q
th
$=$
20, which remains well below 0.5 % of Q
max
. The close agreement between the Q > Q
th
and Q > 0 curves confirm the appropriateness of this threshold. Flow snapshots in figure 14 show the distributions of the integrand m of vortex-induced moment, overlaid with detected vortex structures (black contours), and they clearly reveal the presence of five key vortices: LEV, sLEV, sV, shear-layer vortex (SLV) and TEV. Here, the SLV is defined as a shear-layer separation vortex induced by the surface bulge, originating from the roll-up of the separated shear layer downstream of the bulge, as illustrated in figure 14(e). Figure 13(c, f) illustrates the contributions of these structures to the vortex-induced moment. The discontinuities observed in the curves are attributed to vortex shedding, splitting and merging. Comparing the uncontrolled and controlled cases, the following two primary control pathways emerge.
Flow snapshots comparing the uncontrolled case (a,b) and the active control case (c,d,e) with
$\Delta\phi$
$=$
0.3
$\pi$
at
$A_{\alpha}$
$=$
20°. Black contours indicate detected vortex structures.
Figure 14. Long description
Flow snapshots comparing the uncontrolled case (a,b) and the active control case (c,d,e) with delta phi = 0.3 at A = 20 degrees. Black contours indicate detected vortex structures. In the uncontrolled case, the leading edge vortex (LEV) and trailing edge vortex (TEV) are visible at different times, showing the evolution and shedding of vortices. In the actively controlled case, the snapshots show the LEV and secondary leading edge vortex (sLEV) at various times, with a closer view in panel (e) highlighting the structures at an earlier time point. The comparison illustrates the differences in vortex dynamics between the uncontrolled and controlled scenarios.
-
(i) During the pitching-up phase, surface bulging delays LEV detachment and promotes the formation of a SLV downstream of the bulge, effectively reducing the peak aerodynamic moment.
-
(ii) During the pitching-down phase, surface concave morphing strengthens the sLEV and simultaneously suppresses the TEV, thereby significantly enhancing the moment trough.
Figure 15(b) quantifies the contributions of these five key vortices to
$E_{f}^{\ast}$
at
$A_\alpha$
$=$
20°. It is evident that the suppression of the TEV is the dominant mechanism, accounting for the largest reduction in
$E_{f}^{\ast}$
. To further understand the mechanism behind the suppression of the TEV-induced moment, figure 16(a) provides a detailed quantitative analysis of TEV evolution. Based on (2.14), the significant weakening of the TEV-induced moment under control is primarily attributed to the reduction in the peak averaged Q value. As shown in figure 16(b), the timing of the Q peak reduction coincides with the period of strengthened sLEV-induced moment. Therefore, the enhancement of the sLEV arises from the increased space available due to surface concave morphing, leading to a larger vortex area. In a word, the principal sLEV–TEV control pathway of surface morphing is clarified: concave morphing enhances the sLEV, which triggers a chain effect that suppresses the TEV, resulting in a substantial increase in the aerodynamic moment trough. Unlike conventional strategies that focus on suppressing dynamic stall vortices, this counter-intuitive approach leverages the enhancement of a smaller-scale vortex (sLEV) whose scale is comparable to that of the surface morphing. By targeting the sLEV rather than the large-scale LEV, the control exploits a more accessible and effective pathway for flow control.
Contributions of five key vortex structures to E
f
* as
$A_{\alpha}$
decreases from 25° to 10°, with blue bars representing the uncontrolled cases, red bars representing the active control cases with
$\Delta\phi=0.3\pi$
and yellow bars representing the active control case with optimal phases that correspond to
$E_\textit{min}^{\ast}$
for each
$A_{\alpha}$
.
Figure 15. Long description
The bar graph compares the contributions of five key vortex structures to Ef* as A decreases from 25 to 10. The x-axis represents the vortex structures: LEV, sLEV, sV, SLV, and TEV. The y-axis represents the values of Ef*. There are three sets of bars for each vortex structure, colored blue, red, and yellow. Blue bars represent the uncontrolled cases, red bars represent the active control cases with a phase shift of 0.3pi, and yellow bars represent the active control case with optimal phases corresponding to Emin* for each A. The graph shows how the contributions of each vortex structure change as A decreases. All values are approximated.
Detailed comparison of the (a) TEV and (b) sLEV evolution at
$A_{\alpha}=20^{\circ}$
among the uncontrolled case (blue lines), active control case with
$\Delta\phi$
$=$
0.3
$\pi$
(red lines) and with optimal phases
$\Delta\phi$
$=$
0.42
$\pi$
(yellow lines): vortex-induced moments, average Q values, average potential
$\psi_{k}$
and areas of TEV over time.
Figure 16. Long description
The image contains two sets of graphs labeled (a) and (b), each with four subplots. The graphs compare the evolution of trailing edge vortex (TEV) and stable leading edge vortex (sLEV) at A equals 20 among three cases: uncontrolled (blue lines), active control with delta phi equals 0.3 (red lines), and optimal phases delta phi equals 0.42 (yellow lines). Each set of graphs includes plots for vortex-induced moments, average Q values, average potential k, and areas of TEV over time. The x-axis represents time normalized by T, and the y-axes represent different metrics for each subplot. The graphs show how these metrics vary over time for each case, highlighting differences in vortex-induced moments, average Q values, average potential k, and areas of TEV. All values are approximated.
Detailed comparison of the sLEV evolution at
$A_{\alpha}=20^{\circ}$
among the uncontrolled case (black dotted lines), active control case with
$\Delta\phi$
$=$
0.42
$\pi$
(yellow lines), shape effect case (red lines) and velocity effect case (blue lines): (a) vortex-induced moments, (b) average Q values, (c) average potential
$\psi_{k}$
, (d) areas, (e, f) centroid coordinates of sLEV over time.
Figure 17. Long description
The image contains six line graphs comparing the evolution of sLEV at A equals 20 among different cases: uncontrolled (black dotted lines), fully controlled (yellow lines), shape effect (red lines), and velocity effect (blue lines). Graph (a) shows vortex-induced moments over time, with the y-axis representing the moment coefficient and the x-axis representing normalized time. Graph (b) displays average Q values, with a notable threshold labeled Q_th. Graph (c) presents average potential k values. Graph (d) illustrates the areas of sLEV, with the y-axis representing the maximum y-coordinate. Graphs (e) and (f) show the centroid coordinates of sLEV over time, with the y-axis representing the average x and y coordinates, respectively. Each graph highlights the differences in sLEV evolution among the four cases, with trends, peaks, and variations clearly depicted. All values are approximated.
Finally, as shown in figure 15, the same dominant mechanism is observed for the other pitching amplitudes (
$A_\alpha$
$=$
10°, 15° and 25°), where surface morphing control consistently regulates sLEV development to induce a chain effect that suppresses the TEV and weakens the TEV-induced moment. This confirms that the sLEV–TEV interaction pathway is not limited to a specific case but represents a robust and universal control mechanism within topology I.
To further reveal the physical mechanism by which the unsteady surface morphing enhances the sLEV, two cases (shape effect case and velocity effect case) are implemented to investigate the effect of concave surface shape or the unsteady boundary motion on the enhancement of the sLEV, respectively. In the shape effect case, the shape of the morphing airfoil surface is updated according to the prescribed morphing motion, while the velocity of the airfoil boundary induced by the surface morphing is discarded. On the contrary, the velocity of the airfoil boundary induced by the prescribed surface morphing is considered in the velocity effect case, while the airfoil geometry remains undeformed. In the fully controlled case, the effect of the surface shape and the unsteady boundary motion are considered simultaneously. As shown in figure 17(a), both effects can enhance the sLEV after time/T
$=$
0.43, but the increment of the moment induced by the individual effect is significantly less than that of the combined effects. As shown in figure 17(d), the concave shape in the shape effect case provides sufficient space for the development of the sLEV, which results in the enlargement of the sLEV area. When the morphing surface is concave, the downstream boundary of the sLEV is located near the chordwise centre of the morphing surface. Meanwhile, the sLEV is clockwise vortex. The downwash velocity induced by the unsteady surface morphing in the velocity effect case can strengthen the intensity of the sLEV, which results in the increment of the average Q value and area of the sLEV as shown in figures 17(b) and 17(d). Besides, the downwash velocity can also makes the sLEV closer to the airofoil surface as shown in figure 17(f), which results in the magnitude increment of the auxiliary potential as shown in figure 17(c). Therefore, both the shape and velocity effect of the unsteady surface morphing can enhance the sLEV, but the physical mechanism of two effects is different. The enhancement of the sLEV in the fully controlled case is the result of the combination of the shape effect and velocity effect, which is significantly greater than that of the individual effects.
The enhanced sLEV subsequently leads to the suppression of the TEV. The formation mechanism of the TEV is that the pressure difference between the upper and lower surface near the trailing edge adsorbs the airflow from the lower surface to the upper surface, and the shear layer of trailing edge rolled up to form the TEV. Therefore, the magnitude of the pressure difference between the upper and lower surface near the trailing edge governs the intensity of the TEV. As shown in figure 18(a), under the effect of unsteady surface morphing, the pressure coefficient on the upper surface increases significantly compared with the uncontrolled case. The chordwise velocity u
m
profiles at x
m
/c
$=$
0.8 and 0.9 are shown in figure 18(b). Under the effect of active control, the velocity magnitude of the reverse flow around the trailing edge decreases, especially at wall-normal locations around y
m
/c
$=$
0.1. This leads to an increase in the pressure coefficient on the upper surface of the trailing edge, which in turn suppresses the TEV.
(a) Pressure coefficient distributions on the upper and lower surfaces of the airfoil, and (b) the corresponding velocity profiles at x
m
/c
$=$
0.8 and 0.9 for the uncontrolled case (blue curves) and active control case with
$\Delta\phi$
$=$
0.42
$\pi$
(yellow curves) at
$\textit{time}/T = 0.47$
.
Figure 18. Long description
The image contains two line graphs. The first graph (a) shows pressure coefficient distributions on the upper and lower surfaces of an airfoil. The x-axis represents the normalized chord length (xm/c), ranging from 0 to 1, and the y-axis represents the pressure coefficient (Cp), ranging from -4 to 4. Four lines are plotted: solid blue for the uncontrolled upper surface, dashed blue for the uncontrolled lower surface, solid yellow for the actively controlled upper surface with a phase shift of 0.42, and dashed yellow for the actively controlled lower surface. The second graph (b) displays velocity profiles at two specific chord lengths, xm/c = 0.8 and 0.9. The x-axis represents the normalized velocity (u/U∞), ranging from -1.5 to 0, and the y-axis represents the normalized height (y/c), ranging from 0 to 0.2. Four lines are plotted: solid blue for the uncontrolled case at xm/c = 0.8, dashed blue for the uncontrolled case at xm/c = 0.9, solid yellow for the actively controlled case at xm/c = 0.8, and dashed yellow for the actively controlled case at xm/c = 0.9. The graphs illustrate the differences in pressure and velocity distributions between the uncontrolled and actively controlled airfoil surfaces.
The decrease in the velocity magnitude of the reverse flow around the trailing edge is a result of the enhancement of the sLEV. In the uncontrolled case, the airflow region between the closed LEV region and the closed SV region originates from the trailing edge, whilst the airflow originating from the leading edge can be neglected as shown in figure 19(a). In the active control case, due to the enhancement of the sLEV, a substantially larger amount of near-wall fluid is entrained to the boundary layer, subsequently transported downstream by the SV, and finally convected away from the airfoil by the LEV, as shown in figure 19(b). This mass transport pathway is markedly intensified compared with the uncontrolled case. Figure 19(c) quantitatively compares the mass transport of the uncontrolled and controlled case at the region between the closed LEV region and the closed SV region. The total mass transport of this region in the two cases is basically the same. However, due to the enhancement of the mass transport pathway from the leading edge, the mass transport of this region from the trailing edge decreases significantly. This led to the decrease in the velocity magnitude of the reverse flow around the trailing edge, which in turn suppressed the TEV.
Streamlines overlaid with vorticity fields for (a) the uncontrolled case and (b) active control case with
$\Delta\phi$
$=$
0.42
$\pi$
, and (c) the corresponding mass flow rates.
Figure 19. Long description
The image contains three panels. The first panel shows streamlines overlaid with vorticity fields for an uncontrolled case, highlighting flow from the leading and trailing edges. The second panel displays streamlines overlaid with vorticity fields for an active control case with a phase difference of 0.42π, also indicating flow from the leading and trailing edges. The third panel presents a bar graph comparing mass flow rates per unit spanwise length for the uncontrolled and active control cases. The graph includes total mass flow rates and contributions from the leading and trailing edges.
Migration of the optimal phase: plot of
$A_{\alpha}$
vs
$E_\textit{min}^{\ast}$
(red dashed line) overlaid on the energy map in the
$A_{\alpha}$
–
$\Delta\phi$
parameter space for
$A_m=0.05$
c control. The yellow marker indicates the optimal PLL control result.
Figure 20. Long description
A heat map displays the migration of the optimal phase, represented by a red dashed line overlaid on an energy map in the A parameter space for Am equals 0.05c control. The yellow marker indicates the optimal phase-locked-loop control result. The x-axis represents the phase difference, while the y-axis represents the angle A alpha in degrees. The color scale ranges from blue to red, indicating varying energy values. The heat map shows regions of different energy intensities, with the red dashed line tracing the path of the optimal phase. The yellow marker highlights a specific point of interest on this path.
3.3. Optimal phase migration mechanism and critical indicator
According to the EBM, the optimal phase corresponding to
$E_\textit{min}^{\ast}$
is plotted against
$A_\alpha$
in figure 20. Within the range of
$A_\alpha$
explored in the topology I regime, the optimal phase systematically decreases as
$A_\alpha$
decreases, illustrating a clear migration trend. From a control design perspective, establishing a direct indicator from pre-control flow features to predict the optimal actuation phase without the need for exhaustive parameter space exploration would be highly valuable (Li & Feng Reference Li and Feng2022). This motivates the following analysis to identify such an indicator.
Vortex-induced moment contributions from five key vortex structures: (a–d) uncontrolled cases, (e–h) active control cases with optimal phases corresponding to
$E_\textit{min}^{\ast}$
for each
$A_\alpha$
.
Figure 21. Long description
The image contains eight line graphs arranged in two columns, each showing vortex-induced moment contributions from five key vortex structures. The left column represents uncontrolled cases, while the right column represents active control cases with optimal phases. Each graph plots the moment coefficient (C_M) on the y-axis against normalized time (Time/T) on the x-axis. The graphs compare different angles of attack (A_alpha) and phase differences (Delta_phi) in both uncontrolled and controlled scenarios. The lines represent various vortex structures, including LEV, sLEV, TEV, SLV, and the next cycle LEV. The uncontrolled cases show higher fluctuations in moment coefficients, indicating instability, while the actively controlled cases demonstrate reduced fluctuations, suggesting effective mitigation of stall flutter. All values are approximated.
Figure 21 presents the vortex-induced moment contributions from the five key vortex structures as
$A_\alpha$
decreases from 25° to 10°. Combined with the energy-based analysis of figure 15, it is observed that, at each
$A_\alpha$
, optimal morphing phases primarily enhance control effectiveness via pathway (ii), i.e. suppression of the TEV by strengthening the sLEV-induced chain effect. Building on the understanding from § 3.2, the suppression of the TEV is attributed to the strengthening of the sLEV through concave morphing. Thus, attention is focused on the temporal characteristics of the sLEV and TEV.
The sLEV-induced aerodynamic moment in (a) the uncontrolled case and (b) the active control case with optimal phase, where peak locations are marked by circles.
Figure 22. Long description
The image contains two line graphs side by side. The left graph (a) shows the sLEV-induced aerodynamic moment under uncontrolled conditions, while the right graph (b) shows the moment under active control with optimal phase. Each graph features multiple lines representing different angles: 10 degrees in red, 15 degrees in orange, 20 degrees in green, 25 degrees in blue, and 30 degrees in purple. The x-axis represents time normalized by T, ranging from 0 to 1. The y-axis represents the sLEV-induced aerodynamic moment, ranging from -0.1 to 0.3. Peak locations are marked by circles in both graphs. The uncontrolled case shows more erratic fluctuations, while the active control case demonstrates smoother peaks and troughs, indicating better stability. All values are approximated.
The TEV-induced aerodynamic moment in (a) the uncontrolled case and (b) the active control case with optimal phase.
Figure 23. Long description
The image contains two line graphs side by side. The first graph (a) shows the TEV-induced aerodynamic moment for the uncontrolled case, while the second graph (b) shows the TEV-induced aerodynamic moment for the active control case with optimal phase. Both graphs plot TEV-induced C_M on the y-axis against Time/T on the x-axis. Each graph includes multiple lines representing different angles: 10 degrees, 15 degrees, 20 degrees, 25 degrees, and 30 degrees. In the uncontrolled case, the lines show varying dips and recoveries over time. In the active control case, the lines also show dips but with different patterns compared to the uncontrolled case, indicating the effect of optimal phase control. The graphs highlight how active control with optimal phase alters the aerodynamic moment compared to the uncontrolled scenario. All values are approximated.
Critical indicator analysis: red bars represent the time of the sLEV-induced C
M
peak in the uncontrolled flow; blue bars represent the time when concave morphing is initiated corresponding to the optimal phase. Error bars account for the sampling resolution of
$\Delta\phi$
at intervals of 1/12
$\pi$
.
Figure 24. Long description
The bar graph compares the time when the surface starts to morph concavely and the time when C M induced by s L E V reaches its peak for different A alpha values. The x-axis represents A alpha in degrees, ranging from 10 to 30. The y-axis represents the normalized time, ranging from 0 to 0.5. The graph features two sets of vertical bars: blue bars represent the time when the surface starts to morph concavely corresponding to the optimal phase, and red bars represent the time when C M induced by s L E V reaches its peak. Error bars are included to account for the sampling resolution at intervals of 1/12. The data shows a general trend where both times decrease as A alpha increases. All values are approximated.
(a) Time-domain response of the airfoil pitching amplitude: uncontrolled case (grey line), open-loop control (red line) and PLL control (blue line). Active control is activated at dimensionless time
$t^*=t U_{\infty}/c=150$
. (b) Real-time evolution of airfoil oscillation frequency and phase difference under open-loop control.
Figure 25. Long description
The image contains two line graphs. The first graph (a) shows the time-domain response of the airfoil pitching amplitude for three different cases: uncontrolled (grey line), open-loop control (red line), and phase-locked loop control (blue line). Active control is activated at a dimensionless time tU/c equal to 150. The second graph (b) illustrates the real-time evolution of the airfoil oscillation frequency and phase difference under open-loop control. The x-axis represents time, while the y-axis on the left shows the pitching amplitude in degrees, and the y-axis on the right shows the normalized oscillation frequency. The red line indicates the phase difference, and the blue line represents the normalized oscillation frequency. The graphs highlight the differences in control strategies and their effects on the airfoil’s pitching amplitude and oscillation frequency.
Figure 22 compares the sLEV-induced C
M
evolution between the uncontrolled and optimally controlled cases. In the controlled case, the peak of the sLEV-induced moment is aligned to occur at time/T
$=$
0.5. This timing closely matches the uncontrolled peak of the TEV-induced moment, as shown in figure 23(a), thereby maximising the suppressive effect of the sLEV on the TEV. Based on this observation, figure 24 establishes the relationship between the critical sLEV feature and the timing of concave morphing initiation. The red bars indicate the timing of the sLEV-induced C
M
peak in the uncontrolled flow, while the blue bars denote the timing of concave morphing initiation corresponding to the optimal phase. The error bars account for the finite resolution of phase sampling (Δ
$\phi$
intervals of 1/12
$\pi$
). Within this error range, a strong agreement between the two timings is observed. Therefore, a critical indicator is identified: the optimal timing for initiating concave morphing corresponds to the instant of the sLEV-induced aerodynamic moment peak, ensuring near-optimal control performance. This indicator provides a direct quantitative link between pre-control flow features and the optimal actuation phase, offering substantial advantages for control law design and practical implementation.
3.4. Stall flutter suppression by PLL control
To validate the effectiveness of the sLEV–TEV control pathway and assess the predictive capability of the EBM, fully coupled simulations are carried out using the time-resolved FSC solver introduced in figure 4.
Open-loop control is first examined by prescribing a fixed phase difference of Δ
$\phi$
$=$
0.5
$\pi$
, corresponding to the optimal phase identified for the baseline pitching amplitude of
$A_\alpha$
$=$
24.4°. The time-domain result is shown in figure 25(a) (red line), with the uncontrolled baseline plotted in grey. Initially, the open-loop control indeed reduces the oscillation amplitude of the airfoil. However, the effect deteriorates over time, and the final amplitude slightly exceeds that of the uncontrolled baseline case. This degradation of control effect arises from phase drift, a known limitation of open-loop strategies. Following the real-time estimation method described in § 2.4, in figure 25(b) we present the actual phase difference evolution under open-loop control. Although initially set to Δ
$\phi$
$=$
0.5
$\pi$
, the phase difference gradually decreases due to minor frequency shifts in the airfoil response, leading to accumulated phase error and loss of control effectiveness.
Combined with the EBM, as shown in figure 26, the red trajectory corresponding to the phase evolution in figure 25(b), closely follows the trend predicted by the energy contours. Initially, the trajectory lies within the
$E^{\ast}$
< 0 region, where the oscillation amplitude decays. As the phase drift continues, it crosses the
$E^{\ast}$
$=$
0 stall flutter boundary and enters the
$E^{\ast}$
> 0 region, where amplitude begins to grow. The response eventually converges to a new limit cycle near the
$E^{\ast}$
$=$
0 isoline with a negative energy gradient ∂
$E^{\ast}$
/∂
$A_\alpha$
< 0, matching the observed final amplitude in figure 25(a).
The EBM in the
$A_\alpha$
–Δ
$\phi$
parameter space for A
m
$=$
0.05c control. The dashed line denotes the stall flutter boundary where
$E^{\ast}$
$=$
0. The red solid and blue solid lines represent the pitching response trajectories under open-loop control and PLL control with Δ
$\phi$
$=$
0.5
$\pi$
, respectively.
Figure 26. Long description
A heat map displays the EBM in the A parameter space for Am = 0.05c control. The heat map features a color gradient ranging from blue to red, indicating different values of E*. The x-axis represents the phase difference (Δφ/π) and the y-axis represents the parameter Aα in degrees. The dashed line denotes the stall flutter boundary where E* equals zero. The red solid line represents the pitching response trajectory under open-loop control, while the blue solid line shows the trajectory under phase-locked loop (PLL) control with a specific parameter value. The heat map highlights regions of varying E* values, with blue indicating lower values and red indicating higher values. The trajectories show how different control strategies affect the system’s response.
To address the limitations of open-loop control, a closed-loop control strategy (PLL control) based on the optimal phase
$\Delta\phi$
$=$
0.5
$\pi$
is implemented within the FSC simulator. The resulting time-domain response (figure 25
a, blue line) shows sustained control effectiveness, with the amplitude converging to 14.3°, which is in excellent agreement with the predicted equilibrium point (
$A_\alpha$
$=$
14.2°, Δ
$\phi$
$=$
0.5
$\pi$
) on the EBM (figure 26, blue line). These results confirm that the
$A_\alpha$
–Δ
$\phi$
energy map serves as a reliable predictive tool for PLL control outcomes. To demonstrate the generality of this approach, we examine the following additional cases with varying phase values.
Because of the optimal phase migration with the variation of pitching amplitude
$A_\alpha$
, the optimal phase Δ
$\phi$
$=$
0.5
$\pi$
, which is determined for
$A_\alpha$
$=$
24.4°, does not remain optimal across a broader range of
$A_\alpha$
in PLL control applications. The PLL control performance of Δ
$\phi$
∈ [0,2
$\pi$
] is shown in figure 27 (red line). The best performance is achieved at Δ
$\phi$
$=$
0.3
$\pi$
(highlighted in yellow), with the resulting limit-cycle amplitude reduced to
$A_\alpha$
$=$
9.13°, corresponding to 37.4 % of the uncontrolled amplitude, which is a reduction of 62.6 %. Significant control effects are observed within the range Δ
$\phi$
∈ (0.083
$\pi$
,1.083
$\pi$
]. Within this interval, control performance is non-monotonic, exhibiting two local minima at Δ
$\phi$
$=$
0.3
$\pi$
and 0.83
$\pi$
, respectively. Comparison with the predictions of the EBM (dashed line in figure 27) shows excellent agreement: the converged limit-cycle states across all tested phase values align closely with the predicted stall flutter boundaries. This confirms the predictive validity of the EBM.
The PLL control performance for A
m
$=$
0.05c: time-domain FSC simulation results (red solid line) and corresponding predictions from the EBM (dashed line,
$E^{\ast}$
$=$
0). The optimal control case is highlighted in yellow.
Figure 27. Long description
A heat map displays the performance of PLL control for Am = 0.05c. The x-axis represents the phase difference (Δφ/π) ranging from 0 to 2, while the y-axis represents the actuation amplitude (Aα in degrees) ranging from 0 to 35. The color scale on the right indicates the energy-based map (E*) values, with blue representing lower values and red representing higher values. The red solid line shows the time-domain FSC simulation results, and the dashed black line shows the corresponding predictions from the EBM. The optimal control case is highlighted in yellow. The heat map is divided into two regions: Topology I and Topology II. The plot shows that the PLL control (red line) generally follows the EBM predictions (dashed line) but with some deviations, especially in the lower actuation amplitude range. The optimal control case is located at the intersection of the red line and the yellow highlight.
Figure 28 further examines the variation of two key pitching response features, oscillation frequency
$f_\alpha$
and mean angle
$\alpha$
0, under PLL control. Both quantities are plotted as percentage deviations from the uncontrolled case. Across the entire range of Δ
$\phi$
, changes remain small, with maximum deviations of 3.3 % for
$f_\alpha$
and 2.2 % for
$\alpha$
0. These results validate a key assumption underpinning the energy map approach: the dynamic features
$f_\alpha$
and
$\alpha$
0 of the pitching motion remain approximately invariant under surface morphing, allowing the energy-based framework to reliably predict control outcomes. The energy map suggests a narrow path in the
$\Delta\phi$
–
$A_\alpha$
space where flutter can be fully suppressed (i.e.
$A_\alpha$
→ 0), but FSC simulations do not achieve complete oscillation elimination. This discrepancy is attributed to the slight variation in
$\alpha$
0, as shown in figure 28, which introduces a minor predictive error in the EBM due to its sensitivity to the mean pitch angle. These deviations introduce a limited degree of predictive uncertainty in the EBM, particularly near the transition between flow topologies. Nevertheless, they do not alter the overall prediction trends nor the identification of the optimal control phase, and the EBM remains a useful tool for systematically assessing control effectiveness across the parameter space.
Pitching response features under PLL control: variation of (a) oscillation frequency
$f_\alpha$
and (b) mean pitch angle
$\alpha$
0 with morphing phase difference Δ
$\phi$
. All values are expressed as percentage changes relative to the uncontrolled baseline.
Figure 28. Long description
The line graph presents two data series: one for oscillation frequency (f) and another for mean pitch angle (alpha). The x-axis represents the morphing phase difference (Delta phi/pi), ranging from 0 to 2. The y-axis on the left shows the percentage change in oscillation frequency relative to the baseline, while the y-axis on the right indicates the percentage change in mean pitch angle. The black line with square markers represents the oscillation frequency, peaking at approximately 2.2 percentage at around 0.25 on the x-axis. The blue line with square markers shows the mean pitch angle, reaching a maximum of 3.3 percentage at around 0.15 on the x-axis. Both lines fluctuate and eventually converge towards the baseline as the morphing phase difference increases. All values are approximated.
From the PLL control results, the optimal control phase is consistently close to Δ
$\phi$
$=$
0.3
$\pi$
, where the aeroelastic response is reduced to near the transition boundary between topologies I and II. How does the current control strategy perform in the topology II regime? Taking the small-amplitude case of
$A_\alpha$
$=$
5° as an example, the flow snapshot in figure 7(a) shows that the characteristic sLEV is absent in topology II. Furthermore, without the presence of the sLEV, the principal pathway, where concave morphing enhances the sLEV to initiate a chain effect that suppresses the TEV, cannot be activated.
(a) Comparison of the total moment and TEV-induced moment between the uncontrolled case and active control case with optimal phases Δ
$\phi$
$=$
0.33
$\pi$
under prescribed pitching motion with A
$\alpha$
$=$
5°. (b) Contributions of TEV to
$E_{f}^{\ast}$
, with blue bars representing the uncontrolled cases and yellow bars representing the active control cases with optimal phases that correspond to
$E_\textit{min}^{\ast}$
for each A
$\alpha$
.
Figure 29. Long description
The image features two distinct graphs. The first graph, labeled (a), is a line graph comparing the total moment and trailing-edge vortex (TEV)-induced moment between uncontrolled and actively controlled cases under prescribed pitching motion with an amplitude of 5 degrees. The x-axis represents time normalized by the period T, and the y-axis represents the moment coefficient CM. The uncontrolled case is depicted with a blue dashed line for total moment and a blue solid line for TEV-induced moment, while the controlled case is shown with a red dashed line for total moment and a red solid line for TEV-induced moment. The graph highlights specific values of TEV-induced moment for both cases, with Ef* values of 0.018 for the uncontrolled case and 0.011 for the controlled case. The second graph, labeled (b), is a bar graph showing the contributions of TEV to Ef* for different angles of attack (Aα) ranging from 5 to 25 degrees. The blue bars represent the uncontrolled cases, and the yellow bars represent the actively controlled cases with optimal phases corresponding to the minimum energy transfer (Emin*). The graph indicates the percentage reduction in TEV-induced energy transfer for each angle, with notable reductions of 61 percent at 5 degrees, 4 percent at 10 degrees, 17 percent at 15 degrees, 30 percent at 20 degrees, and 55 percent at 25 degrees. The graphs collectively illustrate the effectiveness of active surface morphing in mitigating stall flutter by reducing the energy transfer from the unsteady flow to the structure.
Figure 29(a) further compares the total moment and the TEV-induced moment between the uncontrolled case and the active control case with optimal phase Δ
$\phi$
$=$
0.33
$\pi$
at
$A_\alpha$
$=$
5°. While the TEV-induced contribution to
$E_{f}^{\ast}$
decreases from 0.018 to 0.011 (a reduction of about 40 %), as shown in figure 29(b), this reduction is significantly smaller than that observed in the topology I region. For instance, at
$A_\alpha$
$=$
10°, which is the nearest case above
$A_\alpha$
$=$
5°, the TEV-induced
$E_{f}^{\ast}$
drops to only 4 % of its uncontrolled value under the same control strategy. These results demonstrate that, without the mediating role of the sLEV, the current single-mode morphing strategy cannot effectively suppress the TEV-induced moment, resulting in a diminished control performance in the topology II region. Alternative mechanisms may be required to achieve further suppression into this regime. Potential directions include passive morphing surface concepts or modifications of the surface morphing characteristics, such as different standing-wave or travelling-wave modes. A detailed investigation of these approaches lies beyond the scope of the present study and is left for future work.
Test region (outlined in purple) selected according to the FSC results.
Figure 30. Long description
A heat map displays the distribution of values with a color gradient ranging from blue to red. The x-axis represents the phase difference in terms of pi, while the y-axis represents the angle of attack in degrees. The color scale on the right indicates the magnitude of a specific energy metric, with blue representing lower values and red representing higher values. A test region is outlined in purple, highlighting a specific area of interest within the heat map. The data points within this region are compared against uncontrolled and phase-locked loop controlled scenarios, showing variations in the energy metric.
3.5. Generality of the sLEV–TEV control pathway
3.5.1. Effect of surface morphing amplitude
The previous sections have established the sLEV–TEV pathway as the dominant control mechanism and demonstrated its effectiveness through the EBM and FSC validation. A remaining question is how sensitive the control performance is to the morphing amplitude A m , and whether the flow-feature-based indicator introduced earlier remains valid when A m is varied.
Figure 30 shows the test region selected for the following study with the EBM constructed for A
m
$=$
0.05c. To examine the influence of the morphing amplitude, additional EBMs are generated for A
m
ranging from 0.01c to 0.055c, as displayed in figure 31. It is observed that as A
m
increases
$E^{\ast}$
systematically decreases across the test region, and the favourable pathway in the energy map emerges and widens. The stall flutter boundary (
$E^{\ast}$
$=$
0) predicted from the energy map indicates that A
m
$=$
0.055c is expected to deliver improved control performance compared with A
m
$=$
0.05c. Importantly, the predicted optimal phases and the critical indicator based on pre-control vortex dynamics remain nearly unchanged, suggesting that their applicability is not strongly affected by the morphing amplitude.
The EBM in the test region for A
m
ranging from 0.01c to 0.055c. The dashed line denotes the stall flutter boundary where
$E^{\ast}$
$=$
0.
Figure 31. Long description
A heat map displays the EBM in the test region for Am ranging from 0.01c to 0.055c. The x-axis represents the phase difference (Delta phi/pi) ranging from 0.3 to 0.8, and the y-axis represents the angle of attack (A alpha in degrees) ranging from 2 to 20. The color scale at the top indicates the EBM values, with blue representing lower values and red representing higher values. The dashed line denotes the stall flutter boundary where E* equals 0. The heat map shows variations in EBM values across different phase differences and angles of attack, with notable regions of high and low values.
To validate these predictions, FSC simulations are performed for A
m
$=$
0.055c. As shown in figure 32, the PLL-controlled response converges to a reduced limit cycle with minimum amplitude A
$_\alpha$
,
min
$=$
8.68°, corresponding to a 64.4 % reduction relative to the uncontrolled baseline. Compared with A
m
$=$
0.05c, this represents an additional 1.8 % improvement. Figure 33 further compares the aeroelastic response histories for the two cases. Increasing the morphing amplitude by 5‰c shortens the stabilisation time by 44.4 %, consistent with the improved control performance predicted by the EBM. These results confirm that the EBM reliably captures the amplitude dependence of the control effect. While increasing the morphing amplitude beyond A
m
$=$
0.05c further accelerates stabilisation, the associated performance gains become marginal. From a practical perspective, A
m
$=$
0.05c therefore represents a favourable choice that balances control effectiveness and implementation feasibility.
The PLL control performance for A
m
$=$
0.055c: energy map contours in the
$A_\alpha$
–Δ
$\phi$
parameter space, with the dashed line representing the predicted
$E^{\ast}$
$=$
0 boundary and the blue solid line representing FSC simulation results. For comparison, the FSC simulation results for A
m
$=$
0.05c (red solid line) are also included.
Figure 32. Long description
A heat map displays energy map contours in the parameter space for PLL control performance with Am = 0.055c. The dashed line represents the predicted E* = 0 boundary, while the blue solid line shows FSC simulation results. For comparison, the red solid line indicates FSC simulation results for Am = 0.05c. The x-axis is labeled with Δφ/π, and the y-axis is labeled with Aα in degrees. The color scale on the right ranges from -0.08 to 0.08, indicating varying energy levels. The heat map reveals distinct regions of high and low energy, with notable trends and clusters that highlight the relationship between the parameters and energy levels.
Time-domain PLL control results at the optimal phases: red line for A
m
$=$
0.05c and blue line for A
m
$=$
0.055c.
Figure 33. Long description
The line graph presents time-domain PLL control results, comparing two different conditions. The x-axis represents time (t) ranging from 100 to 800, while the y-axis represents the angle of attack amplitude (A alpha) in degrees, ranging from negative 20 to 50 degrees. The red line corresponds to Am = 0.05c with a phase difference (Delta phi) of 0.3 pi, and the blue line corresponds to Am = 0.055c with a phase difference (Delta phi) of 0.29 pi. Key points are highlighted, such as the stall flutter at time 150 with A alpha of 24.4 degrees, control times at 359.7 and 647.0, and specific times at 509.7 and 797.0 with corresponding A alpha values of 8.68 and 9.13 degrees respectively. The graph illustrates the dynamic behavior of the system under these conditions.
The underlying mechanisms are examined in figure 34, which shows the vortex-resolved contributions to
$E_{f}^{\ast}$
for effective morphing amplitude A
m
ranging from 0.03c to 0.055c at
$A_\alpha$
$=$
20°. With increasing A
m
, the disturbance induced by morphing is amplified, leading to a systematic reduction of the TEV-induced moment across all cases. The sLEV–TEV pathway remains the dominant mechanism throughout, with concave morphing consistently strengthening the sLEV and thereby suppressing the TEV.
Contributions of five key vortex structures to
$E_{f}^{\ast}$
at
$A_\alpha$
$=$
20° for surface morphing amplitudes A
m
$=$
0.03c, 0.04c, 0.05c and 0.055c, using their respective optimal phases.
Figure 34. Long description
The bar graph compares the contributions of five key vortex structures to Ef* at A = 20 for different surface morphing amplitudes using their respective optimal phases. The x-axis labels the vortex structures as LEV, sLEV, sV, SLV, and TEV. The y-axis represents the values of Ef*. The graph includes five sets of bars for each vortex structure, each set corresponding to different morphing amplitudes: Uncontrolled, Am = 0.03c, Am = 0.04c, Am = 0.05c, and Am = 0.055c. The bars are color-coded: gray for Uncontrolled, red for Am = 0.03c, blue for Am = 0.04c, green for Am = 0.05c, and yellow for Am = 0.055c. The LEV structure shows the highest contributions across all amplitudes, with values decreasing as the morphing amplitude increases. The sLEV structure shows negative contributions, with values becoming more negative as the morphing amplitude increases. The sV structure shows minimal contributions across all amplitudes. The SLV structure shows negative contributions, with values becoming less negative as the morphing amplitude increases. The TEV structure shows positive contributions, with values decreasing as the morphing amplitude increases. All values are approximated.
Time-domain PLL control results for NACA0015 airfoil.
Figure 35. Long description
The line graph displays the time-domain phase-locked loop control results for a NACA0015 airfoil. The x-axis represents time (t) ranging from 100 to 350 units, while the y-axis represents the amplitude of angle of attack (A alpha) in degrees, ranging from negative 20 to 50 degrees. The blue line indicates the amplitude of angle of attack over time. Initially, the amplitude fluctuates significantly, reaching a peak of approximately 40 degrees. Around the 150-unit mark on the x-axis, a region labeled ‘Stall flutter’ shows an amplitude of 22.3 degrees. After this point, the amplitude gradually decreases, stabilizing around 8.19 degrees as phase-locked loop control is activated, indicated by the shaded blue area labeled ‘PLL control on’.
These above observations indicate that the control mechanism identified for topology I is robust with respect to morphing amplitude, and that the flow-feature-based indicator continues to provide a reliable prediction of the optimal phase.
3.5.2. Effect of airfoil geometry: NACA0015
To verify the generalisability of the control effect and the proposed control mechanism for different airfoil geometries, the additional simulations are implemented for a NACA0015 airfoil at Re
$=$
1000. The other parameters remain unchanged. The optimal control phase predicted by the EBM for the NACA0015 airfoil is Δ
$\phi$
$=$
0.3
$\pi$
. The aeroelastic time response under PLL control with optimal control phase is shown in figure 35. Under the control effect of surface morphing with optimal phase, the pitch amplitude of stall flutter decreases from
$A_{\alpha}=22.3^{\circ}$
to
$A_{\alpha}=8.19^{\circ}$
, and the ratio of stall flutter mitigation is 63.3 %. This level of mitigation is comparable to that obtained for the NACA0012 airfoil.
To verify the scalability of the proposed control mechanism under different airfoil shapes, the
$E_{f}^{\ast}$
reduction mechanism for the prescribed pitch motion case with a pitch amplitude of 20° is analysed using the VILPM. The aerodynamic moment hysteresis curves for the uncontrolled case and controlled case with a control phase of Δ
$\phi$
$=$
0.3
$\pi$
are shown in figure 36(a). Similar to the results of the NACA0012 airfoil cases, the
$E_{f}^{\ast}$
reduction under surface morphing control is induced by two changes of the aerodynamic moment hysteresis curve: (i) a reduction of the moment peak, and (ii) an increase in the moment trough. Furthermore, the contributions of the five key vortex structures to the energy transfer
$E_{f}^{\ast}$
are shown in figure 36(b). Similar to the results of the NACA0012 airfoil cases (figure 15), the
$E_{f}^{\ast}$
reduction is mainly induced by the reduced contributions from LEV, sLEV and TEV. Chief among these is the reduction of the TEV-induced contribution induced by the sLEV–TEV chain effect, which decreases by about 35 %. As shown in figure 36(c,d), the concave morphing enhances the sLEV, which triggers a chain effect that suppresses the TEV, resulting in a reduction of the TEV-induced contribution. This demonstrates that the sLEV–TEV chain effect is also the dominant control mechanism for the NACA0015 airfoil similar to the NACA0012 airfoil.
The NACA0015 airfoil: (a) C
M
hysteresis curves with pitching amplitude
$A_{\alpha}$
$=$
20°, and (b) the corresponding contributions of five key vortex structures to
$E_{f}^{\ast}$
. Flow snapshots comparing the uncontrolled case (c) and active control case (d).
Figure 36. Long description
The NACA0015 airfoil: (a) CM hysteresis curves with pitching amplitude A = 20, and (b) the corresponding contributions of five key vortex structures to Ef*. Flow snapshots comparing the uncontrolled case (c) and active control case (d). The image contains four subfigures. Subfigure (a) shows a line graph depicting CM hysteresis curves with pitching amplitude A = 20 degrees. The x-axis represents the angle of attack in degrees, ranging from -10 to 40 degrees. The y-axis represents the moment coefficient CM, ranging from -0.2 to 0.5. Two curves are shown: a dotted blue line for the uncontrolled case and a solid red line for the active control case with a phase shift of 0.3 pi. Subfigure (b) shows a bar graph depicting the contributions of five key vortex structures to Ef*. The x-axis lists the vortex structures: LEV, sLEV, sV, SLV, and TEV. The y-axis represents the normalized contribution Ef*, ranging from -0.1 to 0.2. Two sets of bars are shown: blue bars for the uncontrolled case and red bars for the active control case with a phase shift of 0.3 pi. Subfigure (c) shows a flow snapshot of the uncontrolled case at time/T = 0.55. The image highlights the positions and interactions of the sLEV, LEV, TEV, and sV structures. Subfigure (d) shows a flow snapshot of the active control case at time/T = 0.55. The image highlights the positions and interactions of the sLEV, LEV, and TEV structures. The flow snapshots use color to indicate velocity magnitude, with a scale ranging from -4.0 to 8.0 meters per second.
3.5.3. Effect of Reynolds number: Re
$=$
5000
To verify the generalisability of the control effect and the proposed control mechanism for higher Reynolds number, the additional simulations are implemented for a NACA0012 airfoil at Re
$=$
5000. The other parameters remain unchanged. The optimal control phase predicted by the EBM is Δ
$\phi$
$=$
0.67
$\pi$
. The aeroelastic time response under PLL control with an optimal control phase is shown in figure 37. Under the control effect of surface morphing with optimal phase, the pitch amplitude of stall flutter decreases from
$A_{\alpha}=29.5^{\circ}$
to
$A_{\alpha}=7.86^{\circ}$
, and the ratio of stall flutter mitigation is 73.4 %. This level of mitigation is comparable to that obtained for the NACA0012 airfoil at Re
$=$
1000.
Time-domain PLL control results for Re
$=$
5000 (NACA0012 airfoil).
Figure 37. Long description
The line graph displays the time-domain PLL control results for a Reynolds number of 5000 using a NACA0012 airfoil. The x-axis represents time (t) ranging from 100 to 300 units, while the y-axis represents the angle of attack amplitude (A alpha) in degrees, ranging from negative 20 to 50 degrees. The graph shows a blue line indicating the amplitude of the angle of attack over time. Initially, the amplitude fluctuates significantly around 40 degrees, indicating stall flutter with an amplitude of 29.5 degrees. At around 150 units of time, the PLL control is activated, leading to a gradual decrease in the amplitude of oscillations. By the end of the graph, the amplitude stabilizes around 7.86 degrees. The graph highlights the effectiveness of PLL control in reducing the amplitude of oscillations over time.
To verify the scalability of the proposed control mechanism at higher Reynolds number, the
$E_{f}^{\ast}$
reduction mechanism for the prescribed pitch motion case with a pitch amplitude of 20° at Re
$=$
5000 is analysed using the VILPM. The aerodynamic moment hysteresis curves for the uncontrolled case and controlled case with a control phase of Δ
$\phi$
$=$
0.67
$\pi$
are shown in figure 38(a). Similar to the results of the NACA0012 airfoil cases at Re
$=$
1000, the
$E_{f}^{\ast}$
reduction under surface morphing control is induced by two changes of the aerodynamic moment hysteresis curve: (i) a reduction of the moment peak, and (ii) an increase in the moment trough. Furthermore, the contributions of the five key vortex structures to the energy transfer
$E_{f}^{\ast}$
are shown in figure 38(b). Similar to the results of the NACA0012 airfoil cases at Re
$=$
1000 (figure 15), the
$E_{f}^{\ast}$
reduction is mainly induced by the reduced contributions from LEV, sLEV and TEV. Chief among these is the reduction of the TEV-induced contribution induced by the sLEV–TEV chain effect, which decreases by about 48 %. As shown in figure 38(c,d), the concave morphing enhances the sLEV, which triggers a chain effect that suppresses the TEV, resulting in a reduction of the TEV-induced contribution. This demonstrates that the sLEV–TEV chain effect is also the dominant control mechanism for the Re
$=$
5000 case similar to the Re
$=$
1000 case. It should be noted that Re
$=$
5000 remains far below realistic aerospace conditions (Re > 106), and the present results therefore represent a step towards generalisability.
The NACA0012 airfoil case at
$\textit{Re}=5000$
: (a) C
M
hysteresis curves with pitching amplitude
$A_{\alpha}=20^{\circ}$
, and (b) the corresponding contributions of five key vortex structures to
$E_{f}^{\ast}$
. Flow snapshots comparing the uncontrolled case (c) and active control case (d).
Figure 38. Long description
The image contains four subfigures. Subfigure (a) shows a line graph depicting CM hysteresis curves with pitching amplitude A = 20 degrees, comparing uncontrolled and controlled cases. The x-axis represents the angle of attack in degrees, and the y-axis represents the moment coefficient CM. The uncontrolled case is shown with a dotted blue line, while the controlled case is shown with a solid red line. Subfigure (b) is a bar graph showing the contributions of five key vortex structures to Ef* at a Reynolds number of 5000. The x-axis lists the vortex structures (LEV, sLEV, sV, SLV, TEV), and the y-axis represents the normalized energy contribution Ef*. The uncontrolled case is shown in blue, and the controlled case is shown in red. Subfigures (c) and (d) are flow snapshots comparing the uncontrolled case and the active control case at Time/T = 0.55. The snapshots highlight the positions and interactions of key vortex structures, including sLEV, LEV, and TEV. The uncontrolled case shows a more chaotic flow pattern, while the active control case shows a more organized flow structure. All values are approximated.
4. Conclusions
This study investigates stall flutter on a NACA0012 airfoil at a Reynolds number of 1000 and demonstrates its suppression through phase-optimised surface morphing. By integrating flow-topology analysis, vortex-based load diagnostics and real-time control, the following conclusions are drawn.
-
(i) Two distinct dynamic stall topologies are identified: (I) boundary-layer-eruption-induced LEV shedding with sLEV formation at large pitching amplitudes, and (II) bluff-body-type shedding without a sLEV at small amplitudes. A new circulation scaling of topology I is introduced, providing a consistent framework for characterising LEV shedding across pitching amplitudes.
-
(ii) A counter-intuitive control pathway is uncovered. Strengthening the sLEV, rather than suppressing dynamic stall vortices, triggers a chain effect that suppresses the TEV, increases the aerodynamic moment trough and disrupts the feedback loop sustaining stall flutter. This sLEV–TEV pathway is shown to dominate control effectiveness in topology I regimes and to hold consistently across pitching amplitudes where a sLEV forms.
-
(iii) By combining vortex-induced load partitioning with an energy-based framework, a flow-feature-based indicator is established. Synchronising actuation with the peak sLEV-induced moment yields near-optimal performance, providing the first demonstration of load partitioning techniques guiding active flow control design.
-
(iv) A PLL control framework is implemented to account for time-varying system responses. It maintains optimal phase relations and achieves complete stall flutter suppression in topology I, in close agreement with EBM predictions. In topology II, since the sLEV is absent, the dominant pathway is inactive, reducing control effect and thus delineating the operational boundary of the present strategy.
-
(v) A key scaling result is demonstrated: active surface morphing with an amplitude an order of magnitude smaller than the characteristic LEV scale can nonetheless achieve complete stall flutter suppression in the topology I regime. It provides a new physical insight into how small-scale perturbations can regulate large-scale vortex dynamics in unsteady aerodynamics.
-
(vi) Varying the morphing amplitude confirmed the robustness of the mechanism. Increasing amplitude strengthened sLEV formation, improved TEV suppression and accelerated LCO stabilisation; however, the benefits diminished beyond a certain limit. At A m
$=$
0.05c we achieved a favourable balance between control effectiveness and practical feasibility. The control effectiveness and mechanism are also verified for a different airfoil geometry (NACA0015) and at a higher Reynolds number (Re
$=$
5000), which further demonstrates the generality of the proposed sLEV–TEV control mechanism. Nevertheless, the Reynolds numbers considered (Re
$=$
1000–5000) remain well below practical aerospace regimes. The results therefore represent a step towards, rather than a proof of, full generalisability, extension to higher Reynolds numbers is deferred to future work.
Acknowledgements
The authors thank the anonymous reviewers whose comments and suggestions greatly improved the quality and clarity of this manuscript.
Funding
This work was supported in part by the China National Natural Science Foundation (Grant No. 12302226), the China Postdoctoral Science Foundation under Grant Number BX20230453, the Fundamental Research Funds for Central Universities (Grant No. YWF-23-SDHK-L-002) and the Academic Excellence Foundation of BUAA for PhD Students.
Data availability statement
Data supporting the findings of this study are available from the corresponding author on a reasonable request.
Declaration of interests
The authors report that they have no conflict of interest.
Appendix A. Numerical verification
Our numerical model is validated through a pitching airfoil at Re
$=$
1000 based on Menon & Mittal (Reference Menon and Mittal2019). The airfoil is NACA0015, the elastic axis location of the airfoil is X
ea
$=$
0.5 and the equilibrium angular position is
$\alpha$
e
$=$
15°. The aerodynamic model is validated for the dynamic stall case of an airfoil subjected to forced oscillation with a pitching amplitude of A
$\alpha$
$=$
25°. Concurrently, three sets of O-topology mesh, as detailed in table 1, are generated to verify mesh independence. The computational results for these three meshes, along with a comparison to literature data, are presented in figure 39. The time step used here is 0.002. It can be observed that the medium grid achieves computational convergence and agrees well with the results from the literature. Subsequently, the time steps dt
$=$
0.001, 0.002, 0.005 are selected for time-step-independence verification, as shown in figure 40. Taking accuracy and computational cost into account, time step dt
$=$
0.002 is used in this paper. Furthermore, the FSI method is validated for the stall flutter cases (Menon & Mittal Reference Menon and Mittal2019). Figure 41 shows the pitching amplitude A
$\alpha$
as a function of U
* for the present study, along with a comparison to the data from the literature. The results show a high degree of agreement between the two.
Mesh parameters.
Computational results of different meshes compared with the literature data.
Figure 39. Long description
The line graph compares computational results of different meshes with literature data. The x-axis represents time (t) ranging from 20 to 36, and the y-axis represents the variable C_M ranging from -0.6 to 0.6. The graph includes three data lines: a black line for coarse grid, a red line for medium grid, and a blue line for fine grid. Additionally, pink diamond-shaped data points represent results from Menon & Mittal 2019. The data lines and points show periodic peaks and troughs, indicating cyclical behavior in the results. All values are approximated.
Computational results of different time steps.
Figure 40. Long description
The line graph illustrates the evolution of flow associated with stall flutter for different time steps. The x-axis represents time (t) ranging from 20 to 36, while the y-axis represents the aerodynamic coefficient (C_M) ranging from -0.6 to 0.6. Three data lines are plotted: one in gray for dt equals 0.001, one in red for dt equals 0.002, and one in blue for dt equals 0.005. Each line shows periodic oscillations with peaks and troughs, indicating the dynamic stall process. The gray and red lines closely follow each other, while the blue line shows slightly different behavior, particularly in the amplitude and phase of the oscillations. All values are approximated.
The pitching amplitude A
$\alpha$
as a function of U
* for the present study, compared with the literature data.
Figure 41. Long description
The line graph presents a comparison of pitching amplitude A as a function of U* between the present study and literature data. The x-axis represents U* values ranging from 2 to 12, while the y-axis represents pitching amplitude A values ranging from 0 to 120. The graph includes two data series: one represented by blue squares and a dashed line, labeled ‘Menon & Mittal (2019)’, and the other by red diamonds and a dashed line, labeled ‘Present study’. Both data series show a similar trend, with pitching amplitude A increasing as U* increases. The data points from both studies closely align, indicating consistency between the present study and the literature data. All values are approximated.
Appendix B. The role of E* in EBM
We perform a derivation from the structural dynamics equation perspective to understand how the key energy parameter
$E^{\ast}$
determines the amplitude response of the FSI system. Multiplying (2.3) by
$\dot{\alpha }$
gives:
\begin{equation} I_{\textit{ea}}\ddot{\alpha }\dot{\alpha }+K(\alpha -\alpha _{e})\dot{\alpha }=M_{\textit{ea}}\dot{\alpha }-D\dot{\alpha }^{2}. \end{equation}
Integrating this equation over an oscillation cycle with time period T, where the beginning of the cycle is at time t n , results in an equation for the change of amplitude over the oscillation cycle:
\begin{align} I_{\textit{ea}}\left[{\dot{\alpha }}_{*}^{2}\big(t_{n}+T\big)-{\dot{\alpha }}_{*}^{2}\big(t_{n}\big)\right]+K\left[{\alpha }_{*}^{2}\big(t_{n}+T\big)-{\alpha }_{*}^{2}\big(t_{n}\big)\right]&=\rho {U}_{\infty }^{2}c^{2}\left({E}_{f}^{*}-{E}_{d}^{*}\right)\nonumber\\&=\rho {U}_{\infty }^{2}c^{2}E^{*}. \end{align}
Here
$\alpha _{*}=\alpha -\alpha _{e}$
and
$E^{*}$
is the net energy transfer over one cycle. Assuming a purely sinusoidal form for
$\alpha _{*}$
, i.e.
$\alpha _{*}=A_{\alpha }\sin (\omega t)$
, and
$\omega =2\pi /T$
, results in the following equation for the change in amplitude over an oscillation cycle:
\begin{equation} {A}_{\alpha }^{2}\left(t_{n}+T\right)-{A}_{\alpha }^{2}\left(t_{n}\right)=\frac{\rho {U}_{\infty }^{2}c^{2}E^{*}}{I_{\textit{ea}}\omega ^{2}\cos ^{2}\left(\omega t_{n}\right)+K\sin ^{2}\left(\omega t_{n}\right)}. \end{equation}
To more clearly see the change in amplitude
$A_\alpha$
, it is further derived as
\begin{equation} A_{\alpha }(t_{n}+T)-A_{\alpha }(t_{n})=\xi E^{*}, \end{equation}
where
\begin{align} \xi =\frac{\rho {U}_{\infty }^{2}c^{2}}{\left[I_{\textit{ea}}\omega ^{2}\cos ^{2}\left(\omega t_{n}\right)+K\sin ^{2}\left(\omega t_{n}\right)\right]\left[A_{\alpha }\left(t_{n}+T\right)+A_{\alpha }\left(t_{n}\right)\right]}. \end{align}
This is the Poincare map of this dynamical system. Equilibrium is attained when
$E^{\ast}$
$=$
0. The equilibrium may be stable or unstable, and a linearised stability condition for this system can be derived using (B4):
\begin{equation} \frac{d}{{\rm d}A_{\alpha }}\left(A_{\alpha }+\xi E^{*}\right)\lt 1\Rightarrow \frac{{\rm d}E^{*}}{{\rm d}A_{\alpha }}\lt 0. \end{equation}
Thus, the above indicates that a system allowed to flutter freely will equilibrate only when condition
$E^{\ast}$
$=$
0 and (B6) are satisfied simultaneously. A map of
$E^{\ast}$
can be used to understand and predict the oscillation amplitude change of the FSI system. Here
$E^{\ast}$
< 0 indicates the decay of oscillation amplitude, while
$E^{\ast}$
> 0 denotes the growth of oscillation amplitude.
Appendix C. Derivation of VILPM
We presented a complete derivation of the VILPM here. First, the aerodynamic moment M ea on the airfoil with surface B, immersed in a flow domain V f bounded by surface Σ, consisting of a pressure component and friction component, can be written as
\begin{equation} M_{\textit{ea}}=\boldsymbol{M}\boldsymbol{\cdot }\boldsymbol{e}_{k}=\underset{\textit{pressure}}{\underbrace{\int _{B}p[(\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}})\times \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{k}]{\rm d}S}}+\underset{\textit{friction}}{\underbrace{\int _{B}\nu \rho [(\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}})\times (\boldsymbol{n}\times \boldsymbol{\omega })\boldsymbol{\cdot }\boldsymbol{e}_{k}]{\rm d}S}}, \end{equation}
where p is pressure, ω is the vorticity vector,
X
is the position vector of a point on the airfoil surface,
n
is the normal unit vector of the airfoil surface and
e
k
is the Cartesian basis vector in the k direction. Due to the direction of M
ea
is defined as positive in the clockwise direction, e
k
is set as (0, 0, –1) in this study. An auxiliary potential
$\psi_{k}$
is defined as
\begin{align} {\nabla} ^{2}\psi _{k}=0\quad\text{ with }\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}=\begin{cases} (\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}})\times \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{k} \quad\text{ on }B,\\ 0 \quad\text{ on }\sum .\end{cases} \end{align}
This auxiliary potential is a function of the instantaneous position and shape of the body as well as the outer domain boundary.
Second, the incompressible Navier–Stokes equation is
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}&=0, \end{align}
\begin{align} \frac{\partial \boldsymbol{u}}{\partial t}+\left(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right)\boldsymbol{u}&=-\frac{1}{\rho }\boldsymbol{\nabla }\!p+\nu {\nabla} ^{2}\boldsymbol{u}, \end{align}
where u denotes the fluid velocity. The Navier–Stokes equation is projected onto the gradient of the auxiliary potential and integrated over the volume of the fluid domain V f:
\begin{align} \underset{\left[a\right]}{\underbrace{\int _{V_{f}}\boldsymbol{\nabla }\!p\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}V}}&=-\rho \underset{\left[b\right]}{\underbrace{\int _{V_{f}}\frac{\partial \boldsymbol{u}}{\partial t}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}V}}-\rho \underset{\left[c\right]}{\underbrace{\int _{V_{f}}\left(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right)\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}V}}\nonumber\\&\quad+\nu \rho \underset{\left[d\right]}{\underbrace{\int _{V_{f}}{\nabla} ^{2}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}V}}. \end{align}
For term [a],
\begin{align} \int _{V_{f}}\boldsymbol{\nabla }\!p\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}V&=\int _{V_{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\left(p\boldsymbol{\nabla }\psi _{k}\right){\rm d}V=\int _{B+\sum }p\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}S\nonumber\\&=\int _{B}p\left[\left(\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}}\right)\times \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{k}\right]{\rm d}S ,\end{align}
which is equal to the pressure component of M ea in (C1).
For term [b],
\begin{align} \int _{V_{f}}\frac{\partial \boldsymbol{u}}{\partial t}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}V & =\int _{V_{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\frac{\partial \boldsymbol{u}}{\partial t}\psi _{k}\right){\rm d}V=\int _{B+\sum }\boldsymbol{n}\boldsymbol{\cdot }\left(\frac{\partial \boldsymbol{u}}{\partial t}\psi _{k}\right){\rm d}S\nonumber\\&=\int _{B+\sum }\boldsymbol{n}\boldsymbol{\cdot }\left[\left(\frac{{\rm d}\boldsymbol{u}}{{\rm d}t}-\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\right)\psi _{k}\right]{\rm d}S\nonumber\\& =\int _{B+\sum }\boldsymbol{n}\boldsymbol{\cdot }\left(\frac{{\rm d}\boldsymbol{u}}{{\rm d}t}\psi _{k}\right){\rm d}S-\int _{V_{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\left[\left(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\right)\psi _{k}\right]{\rm d}V . \end{align}
For term [d], according to the incompressibility constraint (C3), it can be derived as
\begin{align} \int _{V_{f}}{\nabla} ^{2}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}V & =\int _{V_{f}}[-\boldsymbol{\nabla }\times (\boldsymbol{\nabla }\times \boldsymbol{u})]\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}V=-\int _{V_{f}}(\boldsymbol{\nabla }\times \boldsymbol{\omega })\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}V\nonumber\\&=\int _{V_{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{\nabla }\psi _{k}\times \boldsymbol{\omega }){\rm d}V =\int _{B+\sum }\boldsymbol{n}\boldsymbol{\cdot }(\boldsymbol{\nabla }\psi _{k}\times \boldsymbol{\omega }){\rm d}S\nonumber\\&=\int _{B+\sum }(\boldsymbol{\omega }\times \boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm.d}S .\end{align}
Therefore, (C5) can be written as
\begin{align} \underset{\textit{pressure}}{\underbrace{\int _{B}p\left[\left(\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}}\right)\times \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{k}\right]{\rm d}S}}&=-\rho \int _{B+\sum }\boldsymbol{n}\boldsymbol{\cdot }\left(\frac{{\rm d}\boldsymbol{u}}{{\rm d}t}\psi _{k}\right){\rm d}S+\rho \int _{V_{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\left[\left(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\right)\psi _{k}\right]{\rm d}V\nonumber\\&\quad -\rho \int _{V_{f}}\!\left(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\right)\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}V +\nu \rho \!\int _{B+\sum }\!\left(\boldsymbol{\omega }\times \boldsymbol{n}\right)\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}S .\end{align}
As
$\boldsymbol{\nabla }\boldsymbol{\cdot }[(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u})\psi _{k}]=[\boldsymbol{\nabla }\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u})]\psi _{k}+(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}$
, then
\begin{align} \underset{\textit{pressure}}{\underbrace{\int _{B}p\left[\left(\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}}\right)\times \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{k}\right]{\rm d}S}}&=-\rho \int _{B+\sum }\boldsymbol{n}\boldsymbol{\cdot }\left(\frac{{\rm d}\boldsymbol{u}}{{\rm d}t}\psi _{k}\right){\rm d}S\nonumber\\&\quad +\rho \!\int _{V_{f}}\!\left[\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\right)\right]\!\psi _{k}{\rm d}V +\nu \rho \!\!\int _{B+\sum }\!\!\left(\boldsymbol{\omega }\times \boldsymbol{n}\right)\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}S . \end{align}
The second term on the right-hand side can be expressed as
\begin{equation} M_{Q}=-\rho \int _{V_{f}}2Q\psi _{k}{\rm d}V, \end{equation}
where
\begin{equation} Q=-\frac{1}{2}\left(\left|\left|\boldsymbol{S}\right|\right|^{2}-\left|\left|{\Omega} \right|\right|^{2}\right)=-\frac{1}{2}\left[Tr\left(\left(\boldsymbol{\nabla }\boldsymbol{u}\right)^{2}\right)-\left(Tr\left(\boldsymbol{\nabla }\boldsymbol{u}\right)\right)^{2}\right]=-\frac{1}{2}\boldsymbol{\nabla }\boldsymbol{\cdot }\left(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\right). \end{equation}
Therefore, the moment in the k direction is partitioned as
\begin{align} M_{\textit{ea}}&=M_{k}+M_{Q}+M_{vp}+M_{v}+M_{\sum }\nonumber\\& =\underset{M_{k}}{\underbrace{-\rho \int _{B}\boldsymbol{n}\boldsymbol{\cdot }\left(\frac{{\rm d}\boldsymbol{u}}{{\rm d}t}\psi _{k}\right){\rm d}S}} \underset{M_{Q}}{\underbrace{-\rho \int _{V_{f}}2Q\psi _{k}{\rm d}V}} +\underset{M_{vp}}{\underbrace{\nu \rho \int _{B}\left(\boldsymbol{\omega }\times \boldsymbol{n}\right)\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}S}}\nonumber\\&\quad +\underset{M_{v}}{\underbrace{\nu \rho \int _{B}\left(\boldsymbol{n}\times \boldsymbol{\omega }\right)\boldsymbol{\cdot }\left[\boldsymbol{e}_{k}\times \left(\boldsymbol{X}-\boldsymbol{X}_{\textit{ea}}\right)\right]{\rm d}S}}\nonumber\\&\quad \underset{M_{\sum }}{\underbrace{-\rho \int _{\sum }\boldsymbol{n}\boldsymbol{\cdot }\left(\frac{{\rm d}\boldsymbol{u}}{{\rm d}t}\psi _{k}\right){\rm d}S+\nu \rho \int _{\sum }\left(\boldsymbol{\omega }\times \boldsymbol{n}\right)\boldsymbol{\cdot }\boldsymbol{\nabla }\psi _{k}{\rm d}S}} ,\end{align}
where M
k
is the kinematic moment, M
Q
is the vortex-induced moment, M
vp
is the viscous-pressure moment, M
v
is the friction moment and
$M_{\Sigma}$
represents the contribution of the velocity field in the outer domain, which is always found to be negligible.
t∗=tU∞/c
Aα
=
Aα
=
Aα
=
Aα
=
α
E∗
=
Aα
=
ϕ
=
π
Aα
ϕ
=
E∗
=
CM
Aα
Δϕ
=
π
Aα
Δϕ
=
π
Aα
=
Δϕ
=
π
Aα
=
Aα
Δϕ=0.3π
Emin∗
Aα
Aα=20∘
Δϕ
=
π
Δϕ
=
π
ψk
Aα=20∘
Δϕ
=
π
ψk
=
Δϕ
=
π
time/T=0.47
Δϕ
=
π
Aα
Emin∗
Aα
Δϕ
Am=0.05
Emin∗
Aα
Δϕ
π
t∗=tU∞/c=150
Aα
ϕ
=
E∗
=
ϕ
=
π
=
E∗
=
fα
α
ϕ
ϕ
=
π
α
=
Ef∗
Emin∗
α
E∗
=
=
Aα
ϕ
E∗
=
=
=
=
Ef∗
Aα
=
=
Aα
=
Ef∗
=
Re=5000
Aα=20∘
Ef∗
α
