2.1. Brief introduction to the model

This study proposes a time-varying system identification modelling framework to address the nonlinear and time-dependent dynamic challenges of transonic buffet outlined in § 1. Combining linear reduced-order modelling with an adaptive time-varying parameter tracking strategy, piecewise state-space equations are constructed to investigate buffet behaviour for the NACA0012 aerofoil under dynamic angle-of-attack sweep conditions. As illustrated in figure 1, the proposed TV-ARX workflow follows three successive functional stages: time-varying input/output signal preparation, adaptive reduced-order identification, and buffet stability evaluation. Time-varying input/output time series comprising current excitation, delayed inputs and previous aerodynamic responses provide the essential information for system identification. These time histories are fed into the reduced-order model (ROM), where the adaptive recursive algorithm continuously updates model coefficients to capture the time-varying aerodynamic characteristics induced by angle-of-attack variation. This coupled procedure produces time-varying state-space representations, whose system matrix varies with time and reflects the evolving aerodynamic stability characteristics. The least stable eigenvalue is extracted at each time instant to quantitatively track how system stability evolves. Migration of this eigenvalue across the stability boundary directly indicates the onset of buffet instability. The present framework enables reliable buffet boundary predictions for continuously swept operating conditions.

Figure 1.

Schematic diagram of the TV-ARX modelling framework.

Figure 1. Long description

Panel A: The diagram begins with time-varying input/output signals, including current excitation, delayed inputs, and previous responses. Panel B: These signals feed into a ROM combined with an adaptive algorithm. Panel C: The output is a time-varying state equation, showing the least stable eigenvalue and the evolution trajectory, indicating stability or instability. Panel D: The process flows into a TV-ARX model, which uses input/output signals to predict buffet boundaries.

2.2. Numerical method and its validation

The URANS equations are employed with an in-house hybrid-unstructured flow solver (Gao et al. Reference Gao, Zhang and Ye2016) for flow field simulation of the NACA0012 aerofoil. The Spalart–Allmaras (S–A) turbulence model is employed for closure, and the governing equations are solved using a cell-centred finite-volume method. The integral formulation of the two-dimensional URANS equations is given by

(2.1) \begin{eqnarray} \frac {\partial }{\partial t}\iiint _{\varOmega }W\,\text{d}\varOmega + \iint _{\partial \varOmega }F(W)\boldsymbol{\cdot }n\,\text{d}S = \iint _{\partial \varOmega }F_{v}n\,\text{d}S . \end{eqnarray}

Here, $W=[\rho ,\rho u,\rho v,\rho e]^{\rm T}$ represents a conserved variable, where $\rho$ denotes fluid density, $u$ and $v$ are the velocities in the $x$ and $y$ directions, respectively, and $e$ is the total internal energy of fluid per unit mass; $F(W)$ represents the convective flux, while $F_{v}$ denotes the viscous flux. Also, $n$ is the normal unit vector outside the boundary of the control volume $\varOmega$ , and $\partial \varOmega$ refers to the cell boundary corresponding to the control volume. Specifically, the diffusion term is discretised using the AUSM+ scheme, whereas the destruction term is approximated by a central scheme. A dual time stepping method is adopted for the unsteady solution of the governing equations. To accelerate convergence during sub-iterations, a four-stage Runge–Kutta scheme is employed along with local time stepping and residual smoothing.

The computational hybrid-unstructured flow solver has 25 361 nodes and 40 layers of structured viscous grid around the aerofoil. The distance between the first layer and the wall in the perpendicular direction is $5\times 10^{-6}$ chords $(y+\sim 1)$ . The time step for the computational fluid dynamics (CFD) simulation is $2.94\times 10^{-4}$ s. The corresponding non-dimensional time step ( $\text{d}t=\text{d}t_{\textit{physics}}a_{\infty }/c$ ) is set to 0.1, where $a_{\infty }$ is the speed of sound, and $c$ indicates the chord of the aerofoil. Further details about grid and computational step convergence are provided in Appendix A. Figure 2 illustrates the buffet onset boundary calculated using the URANS method, which is compared with both simulation results of Soda & Voss (Reference Soda and Voss2005) and experimental data of Doerffer et al. (Reference Doerffer, Hirsch, Dussauge, Babinsky and Barakos2010). Figure 3 compares the root mean square of the fluctuating pressure on the upper profile surface of the OAT15A supercritical aerofoil with the zonal DES results of Deck (Reference Deck2005) and the experimental results of Jacquin et al. (Reference Jacquin, Molton, Deck, Maury and Soulevant2009). Here, $(P_{o})$ is the freestream dynamic pressure. Numerical results obtained with the S–A turbulence model agree well with experimental data for both the buffet onset boundary and aerodynamic loads.

Figure 2.

Comparison of the buffet onset boundary of the NACA0012 aerofoil, $Re=3\times 10^6$ .

Figure 3.

Root mean square (RMS) comparison of fluctuating pressure on the aerofoil upper surface at $ \textit{Ma} = 0.73, \alpha =3.5^{\circ },\ Re=3\times 10^6$ . Here DES is the detached-eddy simulation.

2.3. The TV-ARX identification strategy

A TV-ARX model, which relaxes the constant-coefficient constraint and allows model parameters to evolve dynamically with time, is adopted for the identification of time-varying nonlinear systems.

First, the linear time-invariant ARX (LTI-ARX) model is introduced. It describes the input–output relationship of a system via a discrete difference equation, assuming constant coefficients that are independent of time. For a single-input single-output (SISO) system free from noise interference, its mathematical description is given as

(2.2) \begin{eqnarray} \hat {y}(k) = \sum _{i=1}^{na} A_i\, y(k-i) + \sum _{i=0}^{nb-1} B_i\, u(k-i), \end{eqnarray}

where $u(k)$ and $y(k)$ denote the system excitation input and corresponding output at time step $k$ , respectively. Here, $A_i$ and $B_i$ stand for constant model coefficients to be identified. The integers $na$ and $nb$ denote the output-model order and input-model order, which govern the dependence of the predicted output $\hat {y}(k)$ upon the preceding $na$ output samples and $nb$ input samples. Prior to parameter estimation, it is necessary to remove the constant level from the initial data to ensure that the average data point is zero.

To capture the time-varying behaviour, the LTI-ARX model is extended into a time-dependent modelling framework. A local linear ARX representation is identified at each discrete time step $k$ to match the instantaneous input–output response, producing a family of TV-ARX descriptions. Within this formulation, the model coefficients are no longer constant but evolve over time, as expressed in

(2.3) \begin{eqnarray} \hat {y}(k) = \sum _{i=1}^{na} a_k(i)\, y(k-i) + \sum _{i=0}^{nb-1} b_k(i)\, u(k-i), \end{eqnarray}

where $a_k(i)$ and $b_k(i)$ represent the time-varying coefficients of the TV-ARX model at time step $k$ .

Specifically, a VFF-RLS algorithm (Fortescue et al. Reference Fortescue, Kershenbaum and Ydstie1981; Leung & So Reference Leung and So2005; Paleologu et al. Reference Paleologu, Benesty and Ciochina2008) is employed to track the time-varying model coefficients. The forgetting factor $\mu$ is adaptively tuned according to the model residual, so as to maintain a dynamic balance between accurate tracking of time-varying coefficient behaviour and algorithm numerical robustness. The adaptive update law for the forgetting factor $\mu (k)$ is given as

(2.4) \begin{eqnarray} \mu (k) = \begin{cases} \mu _{\textit{max}}, & \gamma _e(k) \lt 10^{-6}, \\[4pt] \mu _{\textit{min}} + (\mu _{\textit{max}} - \mu _{\textit{min}}) \exp \left ( -{\gamma _e(k)}/{\gamma } \right ), & \text{otherwise}, \end{cases} \end{eqnarray}

where $\mu _{\textit{min}}=0.9$ and $\mu _{\textit{max}}=0.998$ are the lower and upper bounds of the forgetting factor, respectively, $\gamma _e(k)$ denotes the sliding-window standard deviation of noise-free residuals at time $k$ , using window size $N_w=20$ , and $\gamma$ is the scaling coefficient for the residual standard deviation and serves as a sensitivity-tuning parameter.

The regressor vector for VFF-RLS is defined as $\boldsymbol{\varPhi }(k) = [y(k-1), \ldots , y(k-na), u(k),\ldots , u(k-nb+1)]^{\textit{T}}$ . The model predicted output reads $\hat {y}(k) = \hat {\boldsymbol{\theta }}^{\textit{T}}(k-1) \boldsymbol{\varPhi }(k)$ , where $\hat {\boldsymbol{\theta }}(k-1)$ stands for the time-varying coefficient vector of the TV-ARX model obtained at the preceding time step. The main implementation steps are summarised below, centred on recursive iteration for dynamic coefficient updating and correction.

Step 1. Initialisation phase. The time-varying parameter vector $\hat {\boldsymbol{\theta }}(0)$ is initialised to a zero vector to ensure unbiased start-up, and the covariance matrix $\boldsymbol{P}(0)$ is set to a scaled identity matrix $\boldsymbol{P}=10\,000\boldsymbol{I}$ to endow the algorithm with fast initial tracking performance.

Step 2. Recursive update phase. First, the regressor vector $\boldsymbol{\varPhi }(k)$ and physical system output $y(k)$ are acquired at the current time step. The observation error $E(k)$ is then computed to quantify the discrepancy between the model prediction using the preceding coefficient estimates and the actual system response, as defined in

(2.5) \begin{eqnarray} E(k) = y(k) - \hat {\boldsymbol{\theta }}^{\textit{T}}(k-1)\,\boldsymbol{\varPhi }(k) . \end{eqnarray}

Based on the covariance matrix $\boldsymbol{P}(k-1)$ and the regression vector $\boldsymbol{\varPhi }(k)$ , the gain matrix $\boldsymbol{K}(k)$ of the RLS algorithm is updated with the forgetting factor $\mu (k)$ , as given in (this update ensures that the parameter correction is responsive to real-time data while maintaining robustness against noise)

(2.6) \begin{eqnarray} \boldsymbol{K}(k) = \frac {\boldsymbol{P}(k-1)\,\boldsymbol{\varPhi }(k)}{\mu (k) + \boldsymbol{\varPhi }^{\textit{T}}(k)\,\boldsymbol{P}(k-1)\,\boldsymbol{\varPhi }(k)} . \end{eqnarray}

Subsequently, the parameter estimate is corrected using the updated gain matrix and observation error, enabling iterative updates of the TV-ARX model coefficients to capture the time-varying dynamics of the system, as shown in

(2.7) \begin{eqnarray} \hat {\boldsymbol{\theta }}(k) = \hat {\boldsymbol{\theta }}(k-1) + \boldsymbol{K}(k)\,E(k) . \end{eqnarray}

Finally, the covariance matrix is updated to preserve the statistical properties of parameter estimation, provide a reliable basis for the subsequent recursion, and guarantee the numerical stability and convergence of the RLS algorithm, as formulated in

(2.8) \begin{eqnarray} \boldsymbol{P}(k) = \mu (k)^{-1}\left [\boldsymbol{I} - \boldsymbol{K}(k)\,\boldsymbol{\varPhi }^{\textit{T}}(k)\right ]\boldsymbol{P}(k-1) . \end{eqnarray}

Step 3. Result output phase. Upon completion of each recursive iteration, the updated coefficient vector $\hat {\boldsymbol{\theta }}(k)$ is exported. This vector consists of the instantaneous time-varying coefficients $a_k(i)$ and $b_k(i)$ of the TV-ARX model at time step $k$ , directly supporting the modelling and analysis of time-varying systems.

To facilitate systematic stability analysis, the TV-ARX model is reformulated into a state-space representation. According to (2.3), the state vector $\boldsymbol{x}(k)$ is defined as $[y(k-1), \ldots , y(k-na), u(k-1), \ldots , u(k-nb+1)]^{\textit{T}}$ . Accordingly, the state and output governing equations of the TV-ARX model are established in discrete state-space form as

(2.9) \begin{eqnarray} \left . \begin{aligned} \boldsymbol{x}(k+1) &= \tilde {\boldsymbol{A}}_k\, \boldsymbol{x}(k) + \tilde {\boldsymbol{B}}_k\, u(k), \\[4pt]\hat {y}(k) &= \tilde {\boldsymbol{C}}_k\, \boldsymbol{x}(k) + \tilde {\boldsymbol{D}}_k\, u(k), \end{aligned} \right \} \end{eqnarray}

where $\tilde {\boldsymbol{A}}_k, \tilde {\boldsymbol{B}}_k, \tilde {\boldsymbol{C}}_k, \tilde {\boldsymbol{D}}_k$ represent the time-varying state transition matrix, input matrix, output matrix and feedforward matrix at time step $k$ , respectively, and their elements are determined by the time-varying autoregressive coefficients $a_k(i)$ and exogenous input coefficients $b_k(i)$ of the TV-ARX model. The time-varying matrices are defined as follows:

(2.10) \begin{align} \tilde {\boldsymbol{A}}_k = \small \begin{bmatrix} a(1) & a(2) & \ldots & a(na-1) & a(na) & b(1) &b(2) & \ldots & b(nb-2) & b(nb-1) \\ 1 & 0 & \ldots & 0 & 0 & 0 & 0 & \ldots & 0 & 0 \\ \vdots & 1 & \ldots & 0 & 0 & 0 &0 & \ldots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots &\vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \ldots & 1 & 0 & 0 &0 & \ldots & 0 & 0 \\ 0 & 0 & \ldots & 0 & 0 & 0 &0 & \ldots & 0 & 0 \\ 0 & 0 & \ldots & 0 & 0 & 1 &0 & \ldots & 0 & 0 \\ 0 & 0 & \ldots & 0 & 0 & 0 &1 & \ldots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots &\vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \ldots & 0 & 0 & 0 &0 & \ldots & 1 & 0 \end{bmatrix}, \\[-28pt] \nonumber \end{align}

(2.11) \begin{align} \tilde {\boldsymbol{B}}_k = \small \begin{bmatrix} b(0) & 0 & \ldots & 0 & 0 & 1 & 0 &\ldots & 0 & 0 \end{bmatrix}^{\textit{T}}\!, \\[-28pt] \nonumber \end{align}

(2.12) \begin{align} \tilde {\boldsymbol{C}}_k = \small \begin{bmatrix} a(1) & a(2) & \ldots & a(na-1) & a(na) & b(1) & b(2) & \ldots & b(nb-2) & b(nb-1) \end{bmatrix} \!, \\[-28pt] \nonumber \end{align}

(2.13) \begin{align} \tilde {\boldsymbol{D}}_k = b(0). \\[0pt] \nonumber \end{align}

The discrete-time state-space model is further transformed into its continuous-time representation via the bilinear transformation, as expressed in (2.14) below. This transformation allows the system stability to be assessed directly from the eigenvalues of matrix $\boldsymbol{A}$ . Specifically, the stability characteristics of the system are governed by the sign of the real part of the eigenvalue with the largest real part, namely the least stable eigenvalue (Schmid Reference Schmid2010; Kaczorek Reference Kaczorek2016). The system is stable if the real part of the least stable eigenvalue is negative, and unstable otherwise. We have

(2.14) \begin{align} \left . \begin{aligned} &\dot {\boldsymbol{x}}(t) = \boldsymbol{A}\, \boldsymbol{x}(t) + \boldsymbol{B}\, u(t), \\[4pt] &y(t) = \boldsymbol{C}\, \boldsymbol{x}(t) + \boldsymbol{D}\, u(t). \end{aligned} \right \} \end{align}

3.1. Stability analysis criteria

This subsection develops stability analysis criteria for boundary prediction of time-varying systems. Physically reasonable constraints and supporting quantitative evaluation and scoring metrics are integrated to ensure reliable stability tracking and physically consistent analysis results. Let $\lambda _1(t),\lambda _2(t),\ldots ,\lambda _{na+nb-1}(t)$ denote the identified eigenvalues of the TV-ARX model at time $t$ , where $\lambda _i(t)=\sigma _i(t)+j\,\omega _i(t)$ , with $\sigma _i(t)$ and $\omega _i(t)$ being its real and imaginary parts, respectively. As one key physical validity constraint, only eigenvalues falling within a physically meaningful frequency range are retained. Specifically, the imaginary part, which represents the modal frequency, must satisfy $0 \leq \omega _i(t) \leq f_{\textit{max}}$ , where the lower bound is $f_{\textit{min}}=0$ . Among these physically admissible eigenvalues, the least stable eigenvalue is defined as the one with the maximum real part, as given in

(3.1) \begin{align} \lambda _{\textit{max}}(t) = \arg \max _{\lambda _i(t) \in \varOmega _{v\textit{alid}}} \left \{ \text{Re}(\lambda _i(t)) \right \}, \end{align}

where $\varOmega _{v\textit{alid}}$ denotes the set of eigenvalues that satisfy the frequency constraint. Its real and imaginary components are expressed as

(3.2) \begin{eqnarray} \sigma _{\textit{max}}(t) = \text{Re}(\lambda _{\textit{max}}(t)), \quad \omega _{\textit{max}}(t) = \text{Im}(\lambda _{\textit{max}}(t)) . \end{eqnarray}

Considering the initial transient convergence of the TV-ARX model, the frequency constraint is relaxed during the first 30 % of the total time period $t \leq 0.3T_{\textit{total}}$ to accommodate the initial parameter convergence phase. Once the model parameters stabilise, strict physical screening is restored. Finally, the least stable eigenvalue is used to determine system stability, and the criterion is defined as

(3.3) \begin{align} \text{sys}_{stab}(t) = \begin{cases} \text{Stable}, & \sigma _{\textit{max}}(t) \lt 0, \\ \text{Unstable}, & \sigma _{\textit{max}}(t) \geq 0, \end{cases} \end{align}

where

(3.4) \begin{align} \begin{cases} \text{no frequency constraint}, & t \leq 0.3T_{\textit{total}}, \\ 0 \leq \omega _i(t) \leq f_{\textit{max}}, & t \gt 0.3T_{\textit{total}}. \end{cases} \end{align}

Figure 4 schematically depicts the temporal evolution of the real part of the least stable eigenvalue during an upward angle-of-attack sweep, together with the detailed screening procedure for this eigenvalue at a representative time instant. The horizontal axis $t$ in figure 4(a) corresponds to a monotonic upward angle-of-attack sweep, in which the angle of attack increases continuously over time. At early low-angle stages (the green stable region), the real part of the least stable eigenvalue remains negative, indicating stable aerodynamic behaviour. Within this schematic, $t=0.3T_{\textit{total}}$ serves as the activating switch for the frequency constraint. The stability transition occurs when the real part of the least stable eigenvalue crosses zero from negative to positive. The red-shaded region denotes the unstable buffet regime after this zero-crossing event. Figure 4(b) illustrates eigenvalue processing at a representative time step during the upward angle-of-attack sweep. The procedure first applies frequency-filtered screening to discard spurious eigenvalues, and the least stable eigenvalue at this instant is then selected according to its real part value.

Figure 4.

Schematic of the time-varying trajectory for the real part of the least stable eigenvalue and the stepwise screening procedure.

Using the dynamic evolution curve for the real part of the least stable eigenvalue subjected to time-varying input excitation, this paper proposes a comprehensive quantitative index for assessing stability prediction performance. This index incorporates both sign-crossing consistency and smoothness. As one fundamental quantity, the actual number of zero crossings $N_{cross}$ is calculated by detecting sign changes between adjacent points, as given in (3.5) below. In addition, $N_{\textit{ideal}}$ denotes the predefined ideal number of stability boundary, which is determined from the physical mechanism governing the system response. Owing to zero initialisation of the parameter vector $\hat {\boldsymbol{\theta }}$ , at the initial time step, the system eigenvalues start within the stable region with negative real parts. With recursive updates of the VFF-RLS algorithm, the model parameters converge rapidly, and the eigenvalues begin to evolve following the system dynamics. Values of $N_{\textit{ideal}}$ corresponding to the two directions of excitation variation are presented in figure 5. To quantify the agreement between predicted and ideal behaviour, an exponential matching metric is constructed, as given in (3.6):

(3.5) \begin{align} N_{\textit{cross}} = \sum _{t=na+1}^{\textit{total}} \mathbf{1}\left \{ \operatorname {sign}(\sigma _{\textit{max}}(t)) \neq \operatorname {sign}(\sigma _{\textit{max}}(t-1)) \right \}, \\[-28pt] \nonumber \end{align}

(3.6) \begin{align} M = \exp \!\left (-0.5 \left | N_{\textit{cross}} - N_{\textit{ideal}} \right | \right )\!. \\[0pt] \nonumber \end{align}

Figure 5.

Ideal zero-crossing count of the least stable eigenvalue for two excitation ramping directions.

A smoothness score is defined to quantify the fluctuation level of the predicted eigenvalue sequence. To mitigate interference induced by abrupt variations near stability boundaries, data points within a tolerance $tol=5$ around each sign-crossing event are masked. Evaluations are performed only over the remaining unmasked data. The local slope variance is computed via a sliding-window approach with window size $N_{w}=5$ , and further normalised by the squared amplitude of the global mean slope to yield a normalised jitter metric. A Gaussian nonlinear mapping transforms this normalised jitter into the final smoothness score, as given in

(3.7) \begin{align} S = \exp \left ( -\left ( \frac {1}{C} \frac {({1}/{N_w}) \sum _{k=1}^{N_w} \mathrm{Var}\left ( \Delta \sigma _{\textit{max},k:k+w-1} \right )}{\left ( \mathrm{mean}\left ( \left | \Delta \sigma _{\textit{max},v\textit{alid}} \right | \right ) \right )^2} \right )^2 \right )\!, \end{align}

where the coefficient $C=60$ is empirically determined to achieve appropriate sensitivity and discriminative performance.

The overall stability prediction score is the product of the matching metric $M$  and the smoothness metric $S$ , defined as in (3.8) below. A high score under these criteria indicates that the evolution of eigenvalues is physically consistent and realistic, which ensures the validity and fidelity of the identified model. The specific TV-ARX modelling framework and the developed stability analysis criteria are illustrated in figure 6. Here,

(3.8) \begin{align} Score = M \times S. \end{align}

Figure 6.

Schematic of the modelling framework and stability analysis criteria.

3.2. Model verification

To verify the effectiveness and reliability of the proposed method, a time-varying damped single-degree-of-freedom second-order oscillator is employed as a benchmark mathematical model. As a classical prototype of self-excited and resonant dynamics, the time-varying damped second-order oscillator exhibits typical dynamic behaviours such as Hopf bifurcation and hysteretic response during parameter sweeps. The governing equation of the system is expressed as

(3.9) \begin{align} \begin{aligned} m\ddot {y} + c(t)\,\dot {y} + ky = F(t), \end{aligned} \end{align}

where $m$ , $k$ and $c(t)$ denote the mass, linear stiffness and time-varying damping coefficient, respectively; $\ddot {y}(t)$ is the instantaneous acceleration, $\dot {y}(t)$ is the instantaneous velocity, $y(t)$ is the oscillator displacement, and $F(t)$ represents the external excitation.

The time-varying damping coefficient $c(t)$ is defined as the core control term governing the stability changes of the control system, which exhibits a linear attenuation trend relative to the input $U(t)$ . In aeroelastic systems, $U(t)$ denotes the time-varying incoming flow velocity, whose acceleration or deceleration evolution directly governs the variation of structural damping characteristics. Its specific expression is given as

(3.10) \begin{align} \begin{aligned} c(t) = c_0 - k_c\, U(t), \end{aligned} \end{align}

where $c_0$ is the initial damping coefficient, and $k_c$ is the linear attenuation coefficient. The stability boundary of the system is determined by the core condition $c(t)=0$ . This condition depends solely on the intrinsic system parameters $c_0$ and $k_c$ , together with the input $U(t)$ , and is independent of the excitation term $F(t)$ . The theoretical critical input $U_{\textit{theory}} = c_0/k_c$ , which corresponds to the condition of zero system damping, remains constant and provides an undeviated reference for the fundamental stability transition boundary.

To ensure that the system response exhibits favourable linear characteristics within the stable region and clearly shows divergent behaviour in the unstable region, a piecewise nonlinear excitation term $F(t)$ is designed, with its core expression as

(3.11) \begin{align} \begin{aligned} F(t)= \begin{cases} F_0, & U(t) \leq U_{\textit{theory}}, \\[5pt] F_0\left (\dfrac {U(t)}{U_{\textit{theory}}}\right )^2, & U(t) \gt U_{\textit{theory}}, \end{cases} \end{aligned} \end{align}

where $F_0$ is an initial small-amplitude excitation that acts as the trigger for system response, avoiding the zero-response problem caused by the absence of excitation. The parameters of the oscillator system are listed in table 1.

Table 1.

Parameters of the time-varying damping oscillator system.

Table 1. Long description

A table listing parameters of a time-varying damping oscillator system. The table has two columns: Parameter and Value. It contains seven rows with the following data: Row 1: Parameter, Mass m; Value, 2. Row 2: Parameter, Stiffness k; Value, 16. Row 3: Parameter, Initial damping coefficient c 0; Value, 5. Row 4: Parameter, Damping attenuation coefficient k c; Value, 0.8. Row 5: Parameter, Initial small-amplitude excitation F 0; Value, 0.1. Row 6: Parameter, Total simulation time t; Value, 20.

A pair of linearly increasing and decreasing input signals $U(t)$ is designed to induce stability variations through velocity excitation. The core objective is to cover the theoretical critical value $U_{\textit{theory}}=6.25$ , and fully capture the displacement response of the second-order damped oscillator during system stability transitions. The system is initialised with zero displacement $y(0)=0$ and zero velocity $\dot {y}(0)=0$ . The time step is set to $\text{d}t=0.004$ , corresponding to a velocity increment ${\rm d}U=0.002$ per time step. The model is numerically solved using the fourth-order Runge–Kutta (RK4) method. This scheme delivers high computational accuracy and good numerical stability, and it can effectively suppress numerical oscillations triggered by the time-varying damping term. In the numerical simulation, a linearly varying signal $U(t)$ is supplied to the system as the input, and the corresponding displacement response $y(t)$ is recorded as the output. The response curves are shown in figure 7, from which a discrepancy can be observed between the actual stability boundary (identified by the onset of marked response divergence) and the theoretical reference.

Figure 7.

The RK4 simulation results with increasing and decreasing input signals and displacement response.

Table 2.

Modelling errors of time-varying damping oscillator with increasing and decreasing input signals.

Table 2. Long description

The table presents the modeling errors and stability prediction scores for a time-varying damping oscillator under increasing and decreasing input signals, for two model settings, na = nb = 80 and na = nb = 120. For each setting, three quantities are reported: the modeling error (error), the critical error percentage e_crit,1 (%), and the overall stability prediction Score. The rows correspond to the two input signals: Increasing and Decreasing. For na = nb = 80: Increasing; error, 8.65E-4; e_crit,1 (%), 4.48; Score, 0.108; Decreasing; error, 2.42E-4; e_crit,1 (%), −4.16; Score, 0.955. For na = nb = 120: Increasing; error, 7.66E-4; e_crit,1 (%), 3.38; Score, 0.991; Decreasing; error, 2.41E-4; e_crit,1 (%), −21.26; Score, 0.312.

Figure 8.

Evolution of the real part of the least stable eigenvalue for the oscillator system.

Based on the input–output data, time-varying system identification is performed using the TV-ARX approach. The model orders are set to $na=nb=80$ or $120$ . The physically valid frequency range is given by $[f_{\textit{min}},f_{\textit{max}}]=[0,5]$ , and the scaling coefficient is set to $\gamma =0.1$ . Figure 8 shows the evolution of the real part of the least stable eigenvalue throughout the modelling process. Herein, the critical input ( $U_{\textit{crit}}$ ) is defined as the input velocity at which the real part of the least stable eigenvalue crosses zero. This threshold corresponds to the transition boundary between the stable and unstable states of the dynamic system. Table 2 summarises the relative identification errors, boundary prediction errors and corresponding stability scores. The identification error ( $ \textit{error}$ ) is defined in (3.12) below, while the boundary prediction error ( $e_{\textit{crit},1}$ ) is given by (3.13). Here, $y_{\textit{identified}}$ represents the response derived from the ROM at time step $k$ , and $L$ is the length of the input–output dataset. The results indicate that for TV-ARX system identification with both increasing and decreasing input signals, the algorithm maintains satisfactory identification accuracy, and yields accurate predictions of the stability boundary. Under high stability scores, the relative errors remain below 5 %. Furthermore, this actual stability boundary deviates from the theoretical reference in the same sense as the stability boundary predicted by the model. We have

(3.12) \begin{align} \begin{aligned} \textit{error} = \frac {\sum _{k=na+1}^L \bigl | y(k) - y_{\textit{identified}}(k) \bigr |}{\sum _{k=na+1}^L \bigl | y(k) \bigr |}, \end{aligned} \\[-28pt] \nonumber \end{align}

(3.13) \begin{align} \begin{aligned} e_{\textit{crit},1} = \frac {U_{\textit{crit}}-U_{\textit{theory}}}{U_{\textit{theory}}}\ast 100\,\%. \end{aligned} \\[-24pt] \nonumber \end{align}

4.1. Linear evolution

For signal A, the system is treated as a SISO system, where the input is the pitching angle $(\theta )$ , and the system output is the corresponding moment coefficient $(C_m)$ . By constructing the continuous-time state-space equation of the TV-ARX model as (2.14), the least stable eigenvalue of the state matrix can be computed to characterise the stability behaviour of system evolution. For the buffet phenomenon, the imaginary part of the least stable eigenvalue that characterises the oscillatory frequency, is uniformly converted into the dimensionless reduced frequency ( $k_b$ ). Its effective range is limited to $[0,0.5]$ , as specified in (3.1). To evaluate the tracking capability of the TV-ARX model against time-varying system dynamics, two distinct values are assigned to the scaling coefficient $\gamma$ , namely 0.01 and 0.00001. A smaller value $\gamma = 0.00001$ increases sensitivity of the model to time-varying variations, while a larger value $\gamma = 0.01$ yields smoother tracking performance.

Figure 12 illustrates the variation of system identification error and real part of the least stable eigenvalue at the last time step with respect to the delay order under the two scaling coefficients. It can be observed that in the low-order model with $na=nb=5$ , the identification error corresponding to a small $\gamma$ is significantly smaller than that for a large $\gamma$ . Moreover, as the model order increases, the trend of the identification error for the two scaling coefficients is exactly reversed. In terms of stability prediction, when a larger $\gamma$ is adopted, the real part of the least stable eigenvalue tends to converge to the ARX result obtained from the full input–output signal indicated by the dashed line as the model order increases, leading to higher accuracy in stability prediction. A trade-off exists between identification accuracy and stability prediction performance, where a larger $\gamma$ combined with a higher model order yields accurate stability prediction while maintaining a relatively small identification error.

Figure 12.

Identification errors and real parts of the least stable eigenvalue at the last time step for two scaling coefficients with varying model orders.

Figure 13 illustrates the stability prediction scores at different delayed orders and $\gamma$ values. When $na=nb=80$ and $\gamma =0.01$ , the TV-ARX model achieves the highest score in boundary prediction. Under this parameter setting, the real part of the least stable eigenvalue from the time-varying model gradually converges as excitation time elapses, and abrupt jumps in stability vanish. This observation agrees with the time-domain responses. In contrast, low-order modelling fails to capture the global characteristics of the system, yielding repeated transitions between stable and unstable predicted states. The specific time-varying eigenvalue trajectory is shown in figure 31. Compared to monitoring only the final-step eigenvalue, the scoring requirement is more stringent, considering the evolution of the entire time-varying process to determine the accuracy of the prediction. Figure 14(a) compares the predicted moment coefficient obtained with optimal model parameters against the CFD reference solution. Both models achieve good agreement with CFD results, with the identification errors of the two models being 6.58 % and 4.59 %, respectively. Figure 14(b) and table 3 further compare the least stable eigenvalues at the last time step of the TV-ARX and time-invariant ARX models. It is evident that the stability variation identified by the time-varying model eventually converges to the results of the time-invariant model. This demonstrates that this time-varying model achieves accurate identification while exhibiting excellent performance in capturing the linear evolution of buffet dynamics. Furthermore, tests are performed under noisy-measurement conditions. The TV-ARX model can accurately track the evolution of system stability even under considerable measurement noise, as detailed in Appendix C.

Table 3.

Comparison of the least stable eigenvalue for signal A at $ \textit{Ma}=0.7,\ \alpha =5.5^{\circ }$ , $Re=3\times 10^6$ .

Table 3. Long description

A table comparing stability prediction scores for ARX and TV-ARX models. The table has two rows and three columns. The columns are labeled with sigma, omega (k sub b), and error (percent). The row labels are ARX and TV-ARX. Row 1: sigma, 1.30e-2; omega (k sub b), 0.199; error (percent), 6.58. Row 2: sigma, 1.28e-2; omega (k sub b), 0.197; error (percent), 4.59.

Figure 13.

Stability prediction scores of two scaling coefficients with different model orders.

Figure 14.

Comparison of results from the two ARX models for signal A at $ \textit{Ma}=0.7,\ \alpha =5.5^{\circ }$ , $Re=3\times 10^6$ .

Figure 14. Long description

Panel A: A line graph titled ‘Moment coefficient result’ shows the comparison of moment coefficient results over time. The x-axis represents time (t) and the y-axis represents the difference in moment coefficients (Cm - Cm0). Three lines are plotted: a solid black line for CFD, a dashed red line for ARX, and a dashed blue line for TV-ARX. The lines show fluctuations over time, with the TV-ARX line closely following the CFD line. Panel B: A scatter plot titled ‘Eigenvalue distribution’ displays the distribution of eigenvalues. The x-axis represents the real part (σ) and the y-axis represents the imaginary part (ω/kb). Red squares represent ARX and blue diamonds represent TV-ARX. The plot includes a vertical dashed line at σ = 0 and annotations indicating specific points for ARX and TV-ARX.

4.2. Nonlinear evolution

Signal B characterises the evolution of the buffet system following short-duration pulse excitation, which transitions to a nonlinear limit cycle oscillation under zero-input free evolution conditions. For such evolutionary processes, traditional linear models fail to capture the dynamics of linear and nonlinear mechanisms, nor can they accurately describe stability transitions and the emergence of self-sustained oscillations. Accordingly, analysis of dynamical behaviours under signal B offers valuable insights for time-varying modelling.

Based on the results presented in the previous subsection, the parameters $na=nb=80$ and $\gamma =0.01$ are adopted for the TV-ARX modelling. Figure 15(a) presents the predicted evolution of aerodynamic moment throughout the whole process. The identification error for the piecewise input is merely $0.38\,\%$ , demonstrating that the proposed method can accurately capture the time-varying variations in aerodynamic moment. In contrast, the time-invariant ARX model fails to capture the growing trend of the moment coefficient amplitude over time. Its predictions rapidly deviate from the CFD reference data, maintaining only low-amplitude steady oscillations, and failing to reflect the evolution of buffet from the linear stage to the strongly nonlinear stage. Meanwhile, the least stable eigenvalues for both models at the final time step are illustrated in figure 15(b). The time-invariant ARX model, derived from the full dataset, yields a weakly unstable yet nearly neutrally stable system. This static eigenvalue distribution cannot support the self-sustained growth of buffet oscillations. In contrast, the proposed TV-ARX model dynamically adapts to capture the limit cycle characteristics of buffet flow. Notably, the least stable eigenvalue yields a near-zero real part at the final time step, which precisely captures the time-varying aerodynamic dynamics of fully developed transonic buffet.

Figure 15.

Comparison of results from the two ARX models for signal B at $ \textit{Ma}=0.7,\ \alpha =5.5^{\circ }$ , $Re=3\times 10^6$ .

The time evolution of the least stable eigenvalue is presented in figure 16. The TV-ARX model updates time-varying parameters over a short initial period, and thereafter tracks the dynamic variations of the flow system efficiently. After excitation by signal B, transonic buffet develops and evolves naturally. In the linear growth stage, the real part of the least stable eigenvalue is approximately $0.01$ , corresponding to a distinctly unstable flow state that drives a rapid increase in oscillation amplitude. Meanwhile, the imaginary part rapidly converges towards the inherent buffet frequency, indicating successful capture of the primary oscillatory mode. The system subsequently enters the transition stage, where the real part decreases gradually under the effect of nonlinear saturation until fully developed limit cycle oscillation is formed. Under steady limit cycle conditions, the mean value of this real part reaches $6.0\times 10^{-4}$ under the steady limit cycle state, which reflects the marginally stable character of sustained buffet response.

This validates that the developed TV-ARX model possesses excellent accuracy and dynamic tracking performance in characterising transonic buffet nonlinear evolution. Different from time-invariant ARX models, which are only applicable to fixed flow conditions under small-amplitude perturbations and cannot capture the full-range dynamic transition of buffet, the strategy for iterative updating of time-varying parameters embedded in the TV-ARX model can synchronously track the fast linear growth of oscillation amplitude and the slow nonlinear saturation process towards limit cycle oscillation. This yields high-consistency fitting with CFD-derived buffet amplitude and stability evolution results. Such unique dynamic adaptation capability fundamentally enables the model to cope with continuously varying aerodynamic conditions, which is exactly the core technical requirement for dynamic buffet characteristic analysis under angle-of-attack sweeping.

Figure 16.

Time evolution of the least stable eigenvalue for signal B at $ \textit{Ma}=0.7,\ \alpha =5.5^{\circ }$ , $Re=3\times 10^6$ .

Figure 16. Long description

The image contains two line graphs labeled Panel A and Panel B. Panel A, titled ‘Real part,’ shows the real part of the least stable eigenvalue over time. The x-axis represents time (t) with units in seconds, and the y-axis represents the maximum real part of the eigenvalue (σ_max) with units in percent. The graph is divided into three regions: Linear, Transitional, and Nonlinear, marked by vertical dashed lines. The data series are color-coded: blue for the real part and red for the threshold value. Key time points are annotated at t = 190 and t = 291. Panel B, titled ‘Imaginary part,’ shows the imaginary part of the least stable eigenvalue over time. The x-axis represents time (t) with units in seconds, and the y-axis represents the imaginary part of the eigenvalue (k_b) with units in percent. The graph includes a horizontal dashed line at k_b = 0.199 and an inset scatter plot showing the relationship between σ and k_b for different regions: Linear, Transitional, and Nonlinear.

5.1. Rate-dependent hysteresis

To explore the effect of pitch rate on the dynamic hysteresis of the buffet onset boundary, this subsection characterises the full-range dynamic response during the transition from stable conditions to unstable transonic shock buffet, and employs an excitation scheme with a linearly increasing angle of attack. Taking the NACA0012 aerofoil at freestream Mach number $ \textit{Ma}=0.70$ as an example, signal training simulations are conducted by the URANS method based on the steady flow. The angle of attack is linearly increased from $3.8^{\circ }$ to $5.5^{\circ }$ with three distinct pitch rates, including $\Delta \alpha =0.068^{\circ }/(c/a_{\infty })$ , $\Delta \alpha =0.034^{\circ }/(c/a_{\infty })$ and $\Delta \alpha =0.0085^{\circ }/(c/a_{\infty })$ , which cover the range across the shock buffet onset. The non-dimensional time step of the CFD simulation is set to 0.1. Two typical time-varying input–output response curves are shown in figure 17. Figure 18 reveals a clear difference in nonlinear coupling between the two cases by bicoherence analysis. For the high pitching rate case, the mean bicoherence is only $0.27$ , with sparse high-coherence points concentrated at low frequencies. This indicates a low level of overall nonlinear coupling, with no broad regions of strong frequency–phase interaction observed. In contrast, the low pitching rate case yields a significantly higher mean bicoherence of $0.66$ . The spectrum exhibits large, continuous regions of strong coherence, indicating widespread and intense frequency–phase interactions across this frequency range.

Figure 17.

The URANS simulation results with increasing angle of attack at two pitch rates.

Figure 18.

Bicoherence spectrum of the moment coefficient for two sweeping cases.

Figure 18. Long description

Panel A: A heat map showing the bicoherence spectrum of the moment coefficient for a sweeping case with a mean bicoherence of 0.27. The x-axis is labeled f1 (kb) ranging from 0 to 20, and the y-axis is labeled f2 (kb) ranging from 0 to 20. The color scale on the right ranges from 0.2 to 1.0, with blue indicating lower values and red indicating higher values. The data shows a concentration of higher values along the bottom and left edges, with sparse higher values scattered throughout the rest of the map. Panel B: A heat map showing the bicoherence spectrum of the moment coefficient for another sweeping case with a mean bicoherence of 0.66. The x-axis is labeled f1 (kb) ranging from 0 to 20, and the y-axis is labeled f2 (kb) ranging from 0 to 20. The color scale on the right ranges from 0.2 to 1.0, with blue indicating lower values and red indicating higher values. The data shows a large concentration of higher values forming a triangular region in the upper right portion of the map, with lower values concentrated along the bottom and left edges.

Figure 19.

Evolution of the real part of the least stable eigenvalue with increasing angle of attack at $ \textit{Ma}=0.70$ .

A SISO TV-ARX model is developed, with angle of attack ( $\alpha$ ) as the input, and moment coefficient ( $C_m$ ) as the output. The predicted critical angle of attack is defined as the first intersection between the real part of the eigenvalue and the zero axis, as the system evolves from a sustained stable state to an unstable state. Since the onset of transonic aeorfoil buffet is closely associated with global instability (Crouch et al. 2007, Reference Crouch, Garbaruk and Magidov2009), global stability analysis serves as an essential method for investigating the buffet phenomenon. For comparison with the time-varying cases, global stability analysis results are obtained at discrete angles of attack under various Mach numbers (Ma et al. Reference Ma, Gao, Xiong and Zhang2024). The critical boundary prediction deviation is defined as

(5.1) \begin{align} \begin{aligned} e_{crit,2} = \frac {\alpha _{\textit{crit}}-\alpha _{global}}{\alpha _{global}}\ast 100\,\%. \end{aligned} \end{align}

The corresponding real-part trajectories for three pitch rates are presented in figure 19, with a stability prediction score ( $Score$ ) higher than 0.99 achieved primarily by tuning the model order. As the angle of attack increases, the real part of the eigenvalues representing the least stable mode gradually approaches the imaginary axis from the stable region, and finally crosses this stability boundary into the unstable region. Figure 20 presents the predicted buffet onset angles of attack at various pitch rates. It is found that the buffet onset angle of attack predicted at high pitch rate deviates significantly from the global stability result $4.58^{\circ }$ , and instead approaches the maximum value $5.5^{\circ }$ during the angle-of-attack sweep. This indicates that when the pitch rate exceeds a critical threshold, the flow field exhibits strong dynamic lag, and cannot respond promptly to rapid changes in the aerofoil attitude, leading to pronounced unsteady delay. For lower pitch rates, the predicted buffet onset angles of attack are $4.74^{\circ }$ and $4.61^{\circ }$ , which agree well with the global stability prediction, with relative deviations of only 3.5 % and 0.7 %, respectively. Comparative results confirm that a smaller pitch rate advances the onset of shock buffet. This trend is in excellent agreement with the findings of Crouch et al. (Reference Crouch, Ahrabi and Kamenetskiy2024). Relevant tests under varying pitch rates show that discernible oscillations first emerge at a larger angle of attack for higher pitch rate.

Figure 20.

Predicted buffet onset angles of attack with various pitch rates at $ \textit{Ma}=0.70$ .

The pressure distributions at $\alpha = 4.5^\circ$ (steady state) and at buffet onset conditions are presented for three different pitching rates in figure 21. At the fixed angle of attack ( $\alpha =4.5^\circ$ ), the steady-state pressure distributions exhibit negligible differences across varying sweep rates, with only minor deviations localised to the shock region. In contrast, at buffet onset conditions, the pressure profiles differ significantly with sweep rate. As the sweep rate decreases, the shock wave shifts downstream, and the suction peak becomes more pronounced. This behaviour directly demonstrates that higher pitching rates delay buffet onset to higher angles of attack, leading to a more forward shock location at onset compared with lower pitching rates. Collectively, these results indicate that while steady-state flow characteristics at a given angle are largely insensitive to the rate of angle-of-attack variation, buffet onset conditions are strongly rate-dependent.

Figure 21.

Pressure distributions at $\alpha = 4.5^\circ$ (steady state) and the predicted buffet onset conditions with $ \textit{Ma}=0.70$ .

Furthermore, verification tests at various Mach numbers are conducted by varying the pitching rate. For each tested Mach condition, the initial flow field is initialised at angle of attack approximately $1.0^\circ$ below the static buffet onset of the corresponding case. Figure 22 compares the buffet onset boundaries predicted at different pitching rates with the results from global stability analysis and wind tunnel experiments (Doerffer et al. Reference Doerffer, Hirsch, Dussauge, Babinsky and Barakos2010). It is clearly observed that decreasing the pitching rate shifts the buffet onset boundary to smaller angles of attack, and this rate-dependent onset exhibits nonlinearity against Mach number. The maximum discrepancy of buffet onset angle of attack between the highest and lowest pitching rates reaches $0.7^\circ$ across all tested Mach numbers. Once the pitching rate decreases below a critical threshold, the dynamic hysteresis triggered by unsteady pitching motion is substantially suppressed, especially for $ \textit{Ma}=0.70$ and $0.75$ , where the boundary deviation between two disparate pitching rates contracts to below $0.2^\circ$ .

Figure 22.

Comparison of predicted buffet onset boundaries during ascending angle-of-attack sweeps.

5.2. Path-dependent hysteresis

To investigate the path-dependent hysteresis effect induced by the angle-of-attack sweep direction, this subsection characterises the full-range dynamic response during the transition of shock buffet from unstable to stable states, and employs an excitation scheme featuring a linearly decreasing angle of attack. Taking the incoming flow conditions with $M=0.70$ and $0.80$ as examples, signal training simulations are conducted based on the unstable steady flow field. In the simulations, three pitch rates, $\Delta \alpha =0.068^{\circ }/(c/a_{\infty })$ , $\Delta \alpha =0.034^{\circ }/(c/a_{\infty })$ and $\Delta \alpha =0.0085^{\circ }/(c/a_{\infty })$ , are adopted respectively. A linearly decreasing excitation is applied to drive the angle of attack across the buffet onset boundary, and the corresponding typical dynamic response curves from URANS simulations at $M=0.70$ are illustrated in figure 23.

Figure 23.

The URANS simulation results with decreasing angle of attack at two pitch rates for $M=0.70$ .

Figure 24.

Evolution of the real part of the least stable eigenvalue with decreasing angle of attack at $ \textit{Ma}=0.70$ .

Under the condition of $Score$ greater than 0.99, the trajectories of the least stable eigenvalues of the system state matrix under different pitch rates are presented in figure 24. Distinct differences are observed in the predicted buffet boundaries between the ascending and descending angle-of-attack sweeps, as shown in figure 25 and table 4. Specifically, the red diamonds and black squares in the figure denote the buffet onset and extinction during ascending and descending angle-of-attack sweeps, respectively. Notably, the predicted buffet boundaries show a more distinct difference between the two variation directions at large pitch rate, where the onset angle of attack during the ascending sweep is obviously higher than that during the descending sweep. This trend coincides with the results of Nitzsche (Reference Nitzsche2009), who investigated transonic buffet characteristics of the supercritical aerofoil BAC 3-11/RES/30/21 at $ \textit{Ma} = 0.75$ via ascending and descending angle-of-attack sweeps. The previous study also identified higher buffet onset angle of attack during ascending sweep relative to descending sweep, confirming the inherent path-dependent hysteresis of transonic shock buffet stability boundaries. A higher pitch rate restricts the full development of flow nonlinearities, and amplifies hysteretic discrepancies between sweep directions, whereas reducing the pitch rate enables sufficient flow adaptation to angle-of-attack variations, and effectively weakens the hysteresis effect.

Table 4.

Comparison of buffet onset and extinction under ascending and descending angle-of-attack sweeps.

Table 4. Long description

The table presents data on pitch rates and their effects on buffet onset and extinction under ascending and descending angle-of-attack sweeps at two different Mach numbers (Ma = 0.70 and Ma = 0.80). The table has four rows and five columns. The columns are labeled Pitch rate, Ma = 0.70 Ascending (3.8° -> 5.5°), Ma = 0.70 Descending (5.5° -> 3.8°), Ma = 0.80 Ascending (0.2° -> 1.9°), and Ma = 0.80 Descending (1.9° -> 0.2°). The rows are labeled with different pitch rates: Δα = 0.068°/(c/a∞), Δα = 0.034°/(c/a∞), and Δα = 0.0085°/(c/a∞). The values in the table represent the pitch rates under different conditions. For Ma = 0.70 Ascending, the values are 5.34, 4.74, and 4.61. For Ma = 0.70 Descending, the values are 4.11, 4.11, and 4.49. For Ma = 0.80 Ascending, the values are 1.55, 1.22, and 0.94. For Ma = 0.80 Descending, the values are 0.43, 0.85, and 0.92.

Figure 25.

Predicted buffet boundaries between the ascending and descending angle-of-attack sweeps.

Wavelet transform analysis of the moment coefficient during angle-of-attack sweeps at $ \textit{Ma}=0.70,\ \Delta \alpha =0.0085^{\circ }/(c/a_{\infty })$ is presented in figure 26. During the angle-of-attack ascending sweep, the flow gradually evolves from a stable steady regime to an unstable buffet state. As the angle of attack increases and buffet occurs, the dominant frequency rapidly rises from its baseline to an intermediate value $k_b = 0.173$ . With further increase in the angle of attack, the buffet mode continues to strengthen, and shifts its spectral signature, eventually converging to the dominant frequency $k_b = 0.228$ at the maximum angle of attack $5.5^\circ$ . In contrast, during the descending sweep, the flow starts from an unstable steady solution at high angles of attack. As the angle of attack decreases, the dominant buffet frequency ( $k_b = 0.185$ ) is fully established, and peaks in amplitude at the middle stage of the sweep. As the sweep continues further, the frequency gradually reduces, and the flow returns to a stable steady state. The distinct spectral characteristics observed in the ascending and descending sweeps directly reflect the contrasting flow evolution paths under increasing and decreasing angles of attack.

Figure 26.

Wavelet analysis of the moment coefficient during angle-of-attack sweeps at $ \textit{Ma}=0.70,\ \Delta \alpha =0.0085^{\circ }/(c/a_{\infty })$ .

Figure 27 illustrates the evolution of shock wave position during ascending and descending angle-of-attack sweeps. The overall range of shock motion for the bidirectional sweeps remains nearly consistent across different pitch rates. The envelope of shock oscillation is inherently constrained by the fixed angle-of-attack range, and remains insensitive to the temporal flow scaling induced by variations in pitching rate. By contrast, the number of shock oscillation cycles shows a strong dependence on pitching rate, evidently increasing at lower rates. A slower sweep provides sufficient time for unsteady flow evolution and shock buffet development. Noticeable phase offsets and differences in the enclosed area of shock oscillation cycles are observed between ascending and descending sweeps, which are governed primarily by aerodynamic hysteresis. At the same angle of attack, the shock wave sustains distinctly different flow states during upward and downward sweeps, leading to asymmetric cyclic shock motions and rate-dependent discrepancies in the hysteresis loop area. The pitching rate serves as a temporal regulator of unsteady shock oscillations, whereas flow hysteresis essentially dominates the irreversible and asymmetric shock behaviour during bidirectional angle of attack variation. As the pitching rate decreases, the buffet onset shock position captured in the ascending sweep gradually converges to the reference location resolved from global stability analysis of static unstable steady base flows ( $\alpha _{global}=4.58^\circ$ ). With the reduction of pitching rate, unsteady flow perturbations weaken, and the flow field progressively reverts towards quasi-static characteristics. Once the pitching rate falls below a critical magnitude, the buffet onset state during ascending sweep nearly coincides with the static analytical result. Meanwhile, the positional deviation of shock-induced buffet onset between ascending and descending sweeps decreases markedly, which evidences the pronounced attenuation of dynamic hysteresis behaviour.

Figure 27.

Shock wave characteristics of three pitch rates during angle-of-attack sweeps at $ \textit{Ma}=0.70$ .

Figure 27. Long description

Three line graphs depict the shock wave characteristics of three pitch rates during angle-of-attack sweeps. Panel A: A line graph shows the relationship between the angle of attack (alpha) in degrees and the position of the shock wave (x/c). The red solid line represents the ascending pitch rate, and the blue dashed line represents the descending pitch rate. The graph includes labels for buffet onset and static conditions, with a global angle of attack (alpha_global) of 4.58 degrees. Panel B: Another line graph with similar axes and labels as Panel A, but with a different pitch rate. Panel C: A third line graph with the same axes and labels, representing a different pitch rate. Each graph shows the variation in shock wave position with changes in the angle of attack for both ascending and descending pitch rates.

Figure 28.

DMD Analysis of pressure field characteristics during angle-of-attack sweeps at $ \textit{Ma}=0.70,\ \Delta \alpha =0.0085^{\circ }/(c/a_{\infty })$ .

To investigate the flow mechanism of transonic aerofoil buffet hysteresis, DMD is applied to transient pressure field snapshots from both ascending and descending angle-of-attack sweeps at rate $\Delta \alpha =0.0085^{\circ }/(c/a_{\infty })$ . One snapshot segment is defined for each sweep branch to mitigate unsteady effects induced by slow angle-of-attack ramping, and lower the risk of buffet-mode aliasing. These two segments are centred on key flow states, namely buffet onset for the ascending branch, and buffet disappearance for the descending branch. Each segment comprises snapshots covering approximately two dominant buffet cycles. This window selection mitigates the frequency shift and modal distortion arising from the limited spectral resolution of single-cycle segments, while suppressing the low-frequency drift associated with excessively long time spans. The DMD is carried out separately for the two segments to obtain modal features illustrated in figure 28. Clear differences are observed in the zero frequency (time averaged) modes obtained from the ascending and descending sweeps. Under the ascending sweep, the zero-frequency mode captures the time-averaged baseline flow within the transition regime from attached to separated conditions. In contrast, under the descending sweep, the zero-frequency mode exhibits pronounced shock-related pressure perturbations on the aerofoil upper surface, with both suction and compression signatures observable. These features directly preserve the shock position and strength characteristics obtained during the strong buffet stage at high angles of attack, clearly reflecting the hysteresis memory effect of the flow field. During the ascending sweep, the flow evolves from attached to shock oscillation with gradually strengthening nonlinear effects, without forming a mature fundamental-harmonic coupling. The third-order mode dominates the primary buffet frequency ( $k_b=0.190$ ), and high-order disturbances show no effective harmonic matching, resulting in an incomplete multiple-frequency mechanism. The descending sweep preserves the memory of nonlinear coupling from the high-angles-of-attack stage. The second- and third-order modes constitute the dominant fundamental buffet frequencies ( $k_b=0.195$ and $k_b=0.147$ ), and the fourth- and fifth-order modes exhibit typical harmonic frequency responses ( $k_b=0.330$ and $k_b=0.393$ ). These harmonic components closely couple with primary shock oscillation, yielding far stronger nonlinear interaction compared with the ascending branch. Such modal discrepancy essentially accounts for the pronounced hysteretic behaviour of transonic buffet.