1. Introduction
Aircraft engine blades, and their incorporation into a turbofan, turboprop or any other engine configuration, are critical components that have significant effects on the performance and safety of aviation. Concurrently, efforts are made to develop engines that meet the growing demand for enhanced performance while ensuring structural integrity and safety assurance, represented by a series of modern engine development plans proposed by the US military, including the Integrated High-Performance Turbine Engine Technology (IHPTET), Versatile Affordable Advanced Turbine Engines (VAATE) and Adaptive Versatile Engine Technology (ADVENT); see Ballal & Zelina (Reference Ballal and Zelina2004) and Hong & Collopy (Reference Hong and Collopy2005). Thus, blade design prioritises characteristics such as lightweight construction, thin profiles and the capacity to withstand higher pressures, significantly increasing the potential for aeroelastic issues that lead to blade-off accidents. As exemplified by incidents such as the engine failure of British Midland Flight 92 in 1989, resulting in 47 fatalities, blade failures have occurred with concerning frequency in recent years, including the blade-related accidents that have been reported in cases involving Southwest Airlines Flight 1380 in 2018, Delta Air Lines Flight 1425 in 2019 and China Eastern Airlines Flight MU712 in 2023, indicating the necessity for comprehensive investigations and scholarly attention to the aeroelasticity of engine blades. The main phenomenon of aeroelastic instability is flutter, characterised by self-excited vibrations that pose an obstacle to the progress of advanced aircraft engines.
The complex flow around the blade produces a wide range of difficult-to-identify aerodynamic characteristics. Saunders & Marshall (Reference Saunders and Marshall2015) delved into vortex reconnection during vortex cutting by a blade, revealing complex phases and mechanics. Leclercq & De Langre (Reference Leclercq and De Langre2018) investigated the reconfiguration of elastic blades in oscillatory flow and identified four kinematic regimes. Kingan & Parry (Reference Kingan and Parry2020) addressed acoustic theory in contrarotating propellers, emphasising the effects of blade sweep on wake interaction noise. Lengani et al. (Reference Lengani, Simoni, Pralits, Durovic, De Vincentiis, Henningson and Hanifi2022) focused on laminar–turbulent transition around low-pressure turbine blades using proper orthogonal decomposition (POD). Together, these investigations underscore the complexity and further emphasise the challenges associated with identifying properties in coupled aeroelastic problems.
Various investigations have highlighted the influence of shock waves on aeroelastic stability, with findings ranging from destabilising effects behind the shock to stabilising effects of the shock itself (Carstens & Belz Reference Carstens and Belz2001). Vahdati et al. (Reference Vahdati, Sayma, Marshall and Imregun2001) came to the conclusion that the separation region behind the shock had a destabilising effect while the shock itself had a stabilising effect. Srivastava & Keith (Reference Srivastava and Keith2012) discovered that the location of the shock wave and the interblade phase angle (IBPA) of blade motion determined the shock wave's effect on blade stability. However, Aotsuka & Murooka (Reference Aotsuka and Murooka2014) demonstrated the stabilising effect of the separation area brought on by the interaction between the shock and boundary layer. In addition, the tip clearance flow, as demonstrated by researchers, significantly affects aeroelastic stability. Fu et al. (Reference Fu, Wang, Jiang and Wei2015) showed the variation of aerodynamic damping is not monotonic, but a trend of first decrease and then increase with the rising of tip gap size, whereas Sun et al. (Reference Sun, Hou, Zhang and Petrie-Repar2019) found that the global aerodynamic damping for the least-stable IBPA increases with the tip gap height. Zheng et al. (Reference Zheng, Gao, Yang and Xu2020) discovered the intensity of the tip vortex and shock wave are the key factors affecting the aeroelastic stability of the rotor when the tip gap size increases. These factors, alongside considerations such as blade mode and reduced velocity (Zhang, Wang & Xu Reference Zhang, Wang and Xu2013) and twist-induced pressure leading to destabilisation (Dong, Zhang & Lu Reference Dong, Zhang and Lu2023) contribute to the complexity of aeroelasticity analysis in this field. In summary, despite ongoing research efforts, the intricate interplay of these factors in aeroelastic behaviour and flutter mechanisms presents a challenging and evolving area of study.
As determined by the National Transportation Safety Board (NTSB), the 2018 Southwest Airlines case illustrated premature fatigue cracks in CFM56 engines, resulting in blade failure (NTSB 2019). NASA has identified blade vibration-induced fatigue issues that contribute to 10–50 % of engine-related accidents (Ramsey Reference Ramsey2006). Therefore, it is imperative to address a critical structural issue: engines, despite having stability margins, are susceptible to failure caused by the development of fatigue cracks. Such damage within the structure plays a fundamental factor in blade flutter, with the flow field serving as a source of energy absorption for vibrations. Thus, the influence of structural damage on aeroelasticity represents a key consideration deserving further investigation. Kim, Vahdati & Imregun (Reference Kim, Vahdati and Imregun2001) analysed aeroelastic stability for damaged blades with bird-strike impairments, noting stability reductions for moderate-to-severe damage. Muir & Friedmann (Reference Muir and Friedmann2016) investigated the aeroelastic response in fan stages with strike impact deformation. Mai & Ryu (Reference Mai and Ryu2020) examined the effect of damage location on gas turbine performance, observing a thermal stress surge near damage sites. The majority of research on damaged blade aeroelasticity focuses on post-accident (e.g. bird-strike) performance, where the blade's geometry is altered. However, limited attention has been devoted to the aeroelastic instability process induced by internal structural fatigue damage.
The presence of such damage can be defined as changes in the dynamic system and will adversely affect the current or future performance of the system. Thus, research typically begins by focusing on two characteristics of material properties: structural stiffness and strength. Hassiotis (Reference Hassiotis2000) simulated the damage with a known decrease in the stiffness of elements in the frame and developed an algorithm to detect and locate the magnitude and the location of such damage. Xie, Xu & Dai (Reference Xie, Xu and Dai2019) explored the effect of structural damage on the aeroelastic stability of supersonic plates, where structural damage is a local bending stiffness loss with various levels, extents and positions. Fayyadh & Raza (Reference Fayyadh and Raza2022) investigated how different damage levels and scenarios affect dynamic-based and static-based stiffness indices representing stiffness deterioration for reinforced-concrete structures. Regarding strength, Zhang et al. (Reference Zhang, Guo, Xu, Li and Yang2021) summarised the strength loss correction formulae based on stiffness loss under different damage conditions of the hull girder. Gao et al. (Reference Gao, Zhu, Yuan, Wu and Xu2022) presented a pair of residual strength and residual stiffness models based on a unified fatigue damage formula, and predicted fatigue life under variable amplitude cyclic loading for composite laminates. Due to the primary focus on the aeroelasticity of damaged blades, the current research in this paper does not consider changes in strength. Instead, it concentrates on localised damage in specific areas (such as cracks or fatigue spots) that can reduce the stiffness of those particular regions. This approach ensures that the blade maintains its functional shape and performance characteristics.
Achieving these goals depends on precise and efficient numerical simulation to study flutter. Current blade flutter analysis technologies fall into three categories: classic uncoupled, fully integrated and partially integrated (Marshall & Imregun Reference Marshall and Imregun1996). The classic method deals with this issue in an uncoupled manner, known as the energy method, via calculating the sum of work done by unsteady aerodynamic forces when the blade vibrates in a certain mode (Carta Reference Carta1967). Researchers have developed energy methods by solving the Euler/Navier–Stokes (N–S) equations and achieved success (see Clark & Hall Reference Clark and Hall2000; Sanders, Hassan & Rabe Reference Sanders, Hassan and Rabe2004; Vahdati, Simpson & Imregun Reference Vahdati, Simpson and Imregun2011). An increase in the coupling influence among the blades within the airflow necessitates a more comprehensive analysis of the uncoupled approach (see Hall & Silkowski Reference Hall and Silkowski1997; Namba & Nishino Reference Namba and Nishino2006; Kubo & Namba Reference Kubo and Namba2011). The fully integrated method streamlines the aeroelastic problem-solving process by synchronising the integration of fluid and structural equations. Advances in rotor blade research have led to more precise analyses (see Yamamoto & August Reference Yamamoto and August1992; Forsching Reference Forsching1994; Vahdati, Sayma & Simpson Reference Vahdati, Sayma and Simpson2005), yet a simulation requires extensive computational effort and faces convergence issues. The partially integrated method solves the fluid and structural equations through iterative data exchange, hereby optimising computational efficiency. Efforts (see Gnesin, Kolodyazhnaya & Rzadkowski Reference Gnesin, Kolodyazhnaya and Rzadkowski2004; Rzadkowski, Gnesin & Kolodyazhnaya Reference Rzadkowski, Gnesin and Kolodyazhnaya2010; Im & Zha Reference Im and Zha2014) have advanced in understanding the aeroelastic behaviour of turbines under various conditions. The partially integrated method strikes a balance between computational cost and accuracy, making it the preferred high-fidelity computational fluid dynamics (CFD)/computational structural dynamics (CSD) coupling method in this paper.
In summary, there are still several difficulties in analysing the aeroelastic behaviour and flutter mechanism/boundary. Commonly used public examples, such as Rotor37 and Rotor67, along with low-speed experimental benches designed by universities, typically do not exhibit flutter problems. The lack of flutter data may produce conflicting or even differing findings in studies on the main causes of flutter. Some studies only analysed aerodynamic damping decreases. Furthermore, real accidents highlight the need to investigate flutter under fatigue damage. Consequently, this paper first employs the partially integrated CFD/CSD coupling method, starting with static aerodynamic elastic deformation to construct a more accurate dynamic aeroelastic analysis. Subsequently, the damaged blade is modelled with local equivalent stiffness loss, and an investigation is conducted into the effect of damage parameters on aeroelastic stability. Finally, the SPOD method is employed to compare flow field characteristics under flutter and stable conditions, thereby helping to better understand the flutter mechanism of transonic damaged blades. The paper is organised as follows. In § 2, the methodology and process for subsystem modelling and coupling are introduced. Section 3 provides the numerical results and discussions, encompassing validation of models, and static and dynamic aeroelastic analysis conducted with a statically deformed blade aerodynamic configuration. In addition, the influence of damage parameters on the aeroelastic behaviour of rotating blades is explored and the flutter mechanisms under destabilising conditions are investigated via the SPOD method. Finally, § 4 provides the main conclusions of the present study.
2. Numerical methodologies
2.1. Geometry model
NASA Rotor67, a low-aspect-ratio transonic axial-flow fan rotor, is used to investigate the flutter mechanism of damaged blades in this work. Since the healthy Rotor67 never fluttered, it can serve as an ideal baseline for a steady aeroelastic system, facilitating the isolation and analysis of the effect of introduced damage without the confounding effects of previous flutter incidents. According to Strazisar et al. (Reference Strazisar, Wood, Hathaway and Suder1989), a comprehensive dataset is presented, encompassing flow conditions, blade coordinates, and experimental data. The fundamental design parameters are listed in table 1.
Table 1. Basic design parameters of NASA Rotor67.
The flow field boundaries of the blades can be considered both periodic and symmetrical, owing to the rotational symmetry exhibited by the blade turbine. The flow characteristics observed on both sides of a single passage, exhibiting identical periodic relationships, can be regarded as forming the periodic boundary. Hence, the flow field of the whole ring cascade is simplified to a single passage, as depicted in figure 1. In addition, the lengths of the upstream and downstream meshes are extended and coarsened to eliminate the pressure wave reflection at the inlet and outlet boundaries (Dong et al. Reference Dong, Zhang, Zhang, Zhang and Lu2020). It is important to acknowledge that the blade model employed in this study does not address the influence of the hub/disc construction or its method of attachment. The model solely focuses on analysing the response of the blade structure itself. In addition, it is assumed that the grid on the hub remains stationary in the present investigation. Unless otherwise specified, the structural parameters of the blade in this study are taken as density 
$\rho =4440\ {\rm kg}\ {\rm m}^{-3}$, Young's modulus 
$D=142$ GPa and Poisson's ratio 
$\nu =0.3$.
Figure 1. Single-passage mesh overview with an extend coarse inlet and outlet.
2.2. Subsystem model
2.2.1. Aerodynamic model
The three-dimensional Reynolds-averaged Navier–Stokes (RANS) equations in the rotating coordinate system attached to the rotating blade are employed in both steady and unsteady simulations. The rotational angular velocity vector is 
${\boldsymbol {\omega }} = {[\omega,0,0]^{\rm T}}$. The finite-volume method is used to discretise the foregoing governing equation,
\begin{equation} \frac{\partial }{{\partial t}}\int_\varOmega {\boldsymbol{W}} \,{\rm d}\varOmega + \oint_S {{\boldsymbol{F_c}}}\, {\rm d}S + \oint_S {{\boldsymbol{F_v}}}\,{\rm d}S = \int_\varOmega {\boldsymbol{Q}}\,{\rm d}\varOmega, \end{equation}
where 
$\boldsymbol {\varOmega }$ represents the fixed control volume with boundary 
$\boldsymbol {S}$, and the conservative variables 
$\boldsymbol {W}$, the convection flux 
$\boldsymbol {F_c}$, the viscous flux 
$\boldsymbol {F_v}$ and the source term 
$\boldsymbol {Q}$ can be written in the following component forms:
\begin{equation} \left.\begin{gathered}
{{\boldsymbol{W}}} = \left[ {\begin{array}{@{}c@{}} \rho \\ {\rho
{u_i}}\\ {\rho E} \end{array}} \right],\quad {{{\boldsymbol{F}}}_c} =
\left[ {\begin{array}{@{}c@{}} {\rho ({u_i} -
{\omega_j})}\\ {\rho {u_i}({u_i} - {\omega_j}) + \rho
{\delta_{ij}}}\\ {\rho E({u_j} - {\omega_j}) + \rho {u_j}}
\end{array}} \right],\quad {{{\boldsymbol{F}}}_v} = \left[
{\begin{array}{@{}c@{}} 0\\ {{\sigma_{ij}}}\\
{{u_k}{\sigma_{ik}} + \kappa \dfrac{{\partial T}}{{\partial
{x_i}}}} \end{array}} \right],\\ {{\boldsymbol{Q}}} = \left[
{\begin{array}{@{}c@{}} 0\\ 0\\ {\omega \rho {u_3}}\\ { -
\omega \rho {u_2}}\\ 0 \end{array}} \right],
\end{gathered}\right\}
\end{equation}
where 
$\rho$, 
$p$, 
$E$, 
$K$ and 
$T$ represent density, pressure, relative total energy per unit mass, thermal conductivity and temperature, respectively. In addition, 
$r_i$, 
$u_i$ and 
$\omega _i$ refer to the displacement, velocity and velocity relative to the rotating frame of reference in three directions. We use 
${\sigma _{ij}}$ to represent the viscous stress tensor expressed with dynamic viscosity coefficient 
$\mu$, which can be approximated by Sutherland's law as a perfect gas:
\begin{equation} {\sigma_{ij}} = \mu \left( {\frac{{\partial {u_i}}}{{\partial {r_j}}} + \frac{{\partial {u_j}}}{{\partial {r_i}}}} \right) - \frac{2}{3}{\sigma_{ij}}\mu \frac{{\partial {u_k}}}{{\partial {r_k}}}. \end{equation} The discretised equations are temporally integrated using the second-order backward Euler technique, and spatial discretisation is achieved through the utilisation of the second-order upwind approach. All simulations are solved using double precision and the shear stress transport 
$k$–
$\omega$ turbulence model. The stagnation parameter and flow angle are given at the inlet, whereas the radial pressure balancing condition is specified at the outflow. It is worth noting that a dependable determination of the stall point can be achieved by employing the outflow mass boundary condition.
2.2.2. Mechanical model
The structural dynamic response can be achieved by a finite-element analysis (FEA) across spatial variables and implementing temporal solutions with respect to the time variable of the global aeroelastic equations of motion:
\begin{equation}
{\boldsymbol{M}}\ddot{\boldsymbol{u}} + {\boldsymbol{C}}\dot{\boldsymbol{u}} +
{\boldsymbol{K}}\boldsymbol{u} = {{\boldsymbol{F}}_e},
\end{equation}
where 
$\boldsymbol {M}$, 
$\boldsymbol {C}$ and 
$\boldsymbol {K}$ represent the system mass, damping and stiffness matrix, respectively, and 
$\boldsymbol {F}$ is the external load vector of system nodes, which are obtained by integrating their respective unit matrices and vectors. We use 
$\ddot {\boldsymbol{u}}$, 
$\dot {\boldsymbol{u}}$ and 
$\boldsymbol{u}$ to represent the acceleration, velocity and displacement of the structure, and 
$\boldsymbol {C}$, which includes both inherent structural damping and Coriolis force effects, has a limited effect and thus is mostly neglected for a greater safety design margin. Responses are computed at these discrete time points, initially utilising the central difference method, and subsequently refined and structured through the application of the Newmark Beta method.
As illustrated in figure 2, the damage under consideration is equivalent to a stiffness degradation by introducing a modified Young's modulus, denoted as 
$\bar {D}$, in comparison with the healthy one 
$D$. Consequently, a stiffness degradation factor 
$Sr$ is pronounced and defined with 
$Sr = \bar {D}/D$, serving as a quantitative measure of the damage level. It is important to note that a smaller 
$Sr$ value corresponds to a lower stiffness and a more severe level of damage. Moreover, it is essential to identify the potential damage location, which is denoted as 
$h_0$, and the damage extent is defined by 
$h$. To provide a non-dimensional assessment of the damage, two parameters are introduced: 
$h_d$, expressed as the starting position of the damage along the blade height 
$H$, denoted as 
$h_d = h_0/H$; and 
$l_d$, which signifies the damage length relative to the blade height as 
$l_d = h/H$. These three parameters, namely 
$Sr$, 
$h_d$ and 
$l_d$, collectively serve to define and quantify the damage level, position and extent within the blade.
Figure 2. Schematic diagram of structural damage: (a) different fatigue failure positions of blades; (b) definition of damage parameters.
What determines whether the time-domain iteration continues or not is the convergence of displacement in time history. The load on blades encompasses both the inertial force and aerodynamic force. The structural analysis commences by examining the displacement response of the structure following the application of minor disturbances. Notably, the damaged model solely accounts for the deterioration of local material properties, without altering the structural appearance. It is important to emphasise that the primary objective here is to investigate the effect and fundamental aeroelastic mechanisms of the blade and the effect of damage parameters, rather than delving into the intricacies of the damage processes.
2.3. Coupling model
The CFD and CSD are integrated into the time domain via their own independent classic solvers, as described in figure 3. The interface data between the fluid and solid domains are exchanged at each integrated time step. The flow solution is utilised to generate the unsteady aerodynamic load vector 
${{\boldsymbol {F}}_e}$, which is subsequently utilised as a boundary condition for the structural model in order to update the blade position. The aerodynamic mesh is afterwards adjusted to align with the structural motion by employing a grid interpolating technique. The iterative process culminates by computing the updated solution for unsteady flow at the new time step. This makes it possible to determine the unsteady pressures, which thus can be used as boundary conditions for the following time step.
Figure 3. Schematic diagram of time-domain coupling solution.
The entire aeroelastic stability is determined through a two-part process, which is namely the analysis of static and dynamic aeroelasticity. The initial stage involves identifying the optimal aerodynamic configuration for the static aeroelastic analysis, referred to as the configuration during operational conditions while disregarding the thermal effect. Prior to initiating the dynamic aeroelastic computation cycle, a slight perturbation deformation is executed on the static equilibrium shape. The investigation and evaluation of the blade flutter characteristics are conducted by analysing the displacement time history.
2.4. Spectral proper orthogonal decomposition method
The spectral proper orthogonal decomposition (SPOD) method (see Towne, Schmidt & Colonius Reference Towne, Schmidt and Colonius2018; Schmidt & Colonius Reference Schmidt and Colonius2020; Ji et al. Reference Ji, Xie, Zhang and Maqsood2023) achieves unsteady flow analysis by extracting the main components of flow field snapshots from time series data. The transition process from time to frequency domain is shown in figure 4, which consists of overlapping time-domain segmentation, discrete Fourier transform (DFT) and frequency domain reordering. Taking the unsteady pressure data of the flow field as input, the mode extraction starts from the construction of the snapshots matrix, which can be expressed in the form of column vector 
${\boldsymbol {\bar Q}} = [{{{\boldsymbol {q}}_{\boldsymbol {1}}}{{,}}{{\boldsymbol {q}}_{\boldsymbol {2}}} {{,}}\ldots {{,}}{{\boldsymbol {q}}_{\boldsymbol {J}}}}]$. With Welch's method, the snapshot matrix is divided into multiple overlapping block matrices, and the 
$n$th block snapshot matrix can be defined as
\begin{equation} {{\boldsymbol{\bar{Q}}}^{(n)}} = \left[ {{\boldsymbol{q}}_1^{(n)},{\boldsymbol{q}}_2^{(n)}, \ldots ,{\boldsymbol{q}}_{{N_f}}^{(n)}} \right],\quad n=1,2,\ldots,{N_b}, \end{equation}
where 
$N_b$ is the number of blocks and 
$N_f$ is the number of snapshots in each block. Taking 
$N_o$ as the number of snapshots overlapped between each block, the 
$k$th snapshot of the 
$n$th block is
\begin{equation} {\boldsymbol{q}}_k^{(n)} = {{\boldsymbol{q}}_{k + (n - 1)\left( {{N_f} - {N_o}} \right)}},\quad k=1,2,\ldots,{N_f}. \end{equation}
Figure 4. Illustration of SPOD method modal extraction (Ji et al. Reference Ji, Xie, Zhang and Maqsood2023).
In the follow-up solution, 
$N_f=2^n$, 
$n=4,5,\ldots 8$ and 
$N_o=N_f/2$, which change with snapshot data length with the need to guarantee 
$N_b>20$. For each block, DFT is used to obtain the new snapshots matrix, which is recorded as
\begin{equation} {{\boldsymbol{\hat{Q}}}^{(n)}} = \left[ {{\boldsymbol{\hat{q}}}_1^{(n)},{\boldsymbol{\hat{q}}}_2^{(n)}, \ldots ,{\boldsymbol{\hat{q}}}_{{N_f}}^{(n)}} \right],\quad n=1,2,\ldots,{N_b}. \end{equation}The present snapshots of the blocks are
$$\begin{gather} \hat{\boldsymbol{q}}_{k}^{(n)}=\frac{1}{\sqrt{N_{f}}} \sum_{j=1}^{N_{f}} w_{j} \boldsymbol{q}_{j}^{(n)} \exp({-\mathrm{i} 2 {\rm \pi}(k-1)})\left[\frac{(j-1)}{N_{f}}\right] ,\quad k=1, \ldots, N_{f}, n=1, \ldots,N_{b}, \end{gather}$$
$$\begin{gather}w(j) = {a_0} - (1 - {a_0})\cos \left( {\frac{{2{\rm \pi} j}}{{{N_f} - 1}}} \right),\quad {0 \leqslant j \leqslant {N_f} - 1}. \end{gather}$$ The weight 
$w_j$ is the value of the window function, which is used to reduce the spectrum leakage caused by the aperiodicity of each block. Here, 
${\boldsymbol {\hat {q}}}_k^{(n)}$ is the Fourier transform at the frequency 
$f_k$ in the 
$n$th block, the corresponding frequency 
$f_k$ is
\begin{equation} f_{k}= \begin{cases}\dfrac{k-1}{N_{f} \Delta t} & k \leqslant N_{f} / 2 ,\\ \dfrac{k-1-N_{f}}{N_{f} \Delta t} & k>N_{f} / 2.\end{cases} \end{equation} The correlation matrix at 
$f_k$ is
$$\begin{gather} {{\boldsymbol{S}}_{{f_k}}} = \frac{{\Delta t}}{{s{N_b}}}\sum_{n = 1}^{{N_b}} {{\boldsymbol{\hat{q}}}_k^{(n)}} {\left( {{\boldsymbol{\hat{q}}}_k^{(n)}} \right)^{{H}}}, \end{gather}$$
$$\begin{gather}s = \sum\limits_{j = 1}^{{N_f}} {w_j^2}. \end{gather}$$ The Fourier components of each block are arranged in ascending order of frequency 
$f_k$, and a new snapshot matrix is introduced:
\begin{equation} {{\boldsymbol{\hat{Q}}}_{{f_k}}} = \sqrt{\chi} \left[ {{\boldsymbol{\hat{q}}}_k^{(1)},{\boldsymbol{\hat{q}}}_k^{(2)}, \ldots ,{\boldsymbol{\hat{q}}}_k^{\left( {{N_b}} \right)}} \right], \end{equation}
where 
$\chi = \Delta t/(s{N_b})$ with 
$\Delta t$ as the sampling time step, and thus the correlation matrix at 
$f_k$ is
\begin{equation} {{\boldsymbol{S}}_{{f_k}}} = {{\boldsymbol{\hat{Q}}}_{{f_k}}}{\left( {{{{\boldsymbol{\hat{Q}}}}_{{f_k}}}} \right)^{{{H}}}}. \end{equation} The spectral proper orthogonal decomposition modes (SPOMs) are obtained by solving the eigenvalue problem for each frequency-dependent matrix 
${\boldsymbol {S}}_{{f_k}}$.
3. Results and discussion
In order to examine the aerodynamic stability of NASA Rotor67, it is necessary to simulate the working line at the designated operational velocity of 16 043 rpm. Based on the current velocity, the Mach number at the blade tip is calculated to be 1.38. At the design condition, the flow rate is recorded as 
$33.15\ \textrm {kg}\ \textrm {s}^{-1}$, whereas the pressure ratio is measured to be 1.63. Furthermore, the flow rate at the choke point on the work line is 
$34.96\ \textrm {kg}\ \textrm {s}^{-1}$. A time-domain analysis is performed to investigate the effect of different operating circumstances on the working thread by altering the exit pressure. The temporal increment employed for the subsequent computation is experientially taken as 
$\Delta t = 1.25 \times 10^{-5}\ \textrm {s}$, which denotes 1/200 of the reciprocal of the first-order natural frequency.
3.1. Validation of coupling model
Initially, a comprehensive examination of mesh dependency is undertaken, concurrently with an assessment of solver accuracy. This entails the utilisation of multiple meshes, which are subsequently juxtaposed against standard experimental conditions (Strazisar et al. Reference Strazisar, Wood, Hathaway and Suder1989). The primary focus is on two pivotal parameters: the total pressure ratio 
${P_{rc}}$ and adiabatic efficiency 
${\eta _c}$, as illustrated in figure 5:
\begin{equation} {P_{rc}} = \frac{{{P_{t2}}}}{{{P_{t1}}}},\quad {\eta_c} = \frac{{1 - {{({{{P_{t2}}}/{{P_{t1}}}})}^{{{{(\gamma - 1)}}/{\gamma }}}}}}{{1 - {{{T_{t2}}}/ {{T_{t1}}}}}}. \end{equation}The prevailing approach utilised in the examination of aerodynamic performance involves the treatment of structures as rigid entities, facilitating the exploration of flow field characteristics around them through steady simulations. Consequently, the initial processing of the conventional static tip clearance of 1.016 mm has been executed. The outcomes demonstrate that the mesh comprising 665 414 cells attains a favourable balance between computational precision and cost. Nonetheless, a disparity is shown while contrasting the obtained results with the experimental data. The findings reveal that the rigid blades portray a decreased total pressure level near the stall point (figure 5a), a diminished choke mass flow rate and a forward shift of the near-peak efficiency point (figure 5b), which does not account for aeroelastic effects present in the actual experiments.
Figure 5. Mesh-dependence test based on a rigid blade: (a) total pressure ratio; (b) adiabatic efficiency.
Figures 6 and 7 delineate the flow field characteristics observed on the typical surface, i.e. 10 % span from the shroud. These figures facilitate a comparative analysis between the computed results and empirical observations. An evident rise of Mach number distribution is observed at the leading edge (LE; figure 6b) for working status near the peak efficiency point. Furthermore, complications are observed in the maximum Mach number between blades for the conditions near stall (figure 7b). These circumstances directly contribute to discrepancies in the placement of shock waves and the magnitude of aerodynamic pressure. As a result, the deformation of the elastic blade is ascertained through the application of static aeroelasticity, which hinges on the integration of CFD and CSD. In distinction from conventional methodologies, which often provide a one-step solution for force-induced deformations, the current approach employs the developed coupling model to advance the analysis in the time domain. Three operating conditions as specified in table 2 are chosen to represent typical performance for the static aeroelastic solution.
Figure 6. Contour plots of Mach number near the peak efficiency point (10 % span from the shroud): (a) experiment (Strazisar et al. Reference Strazisar, Wood, Hathaway and Suder1989); (b) uncoupled blade; (c) coupled blade at Case B.
Figure 7. Contour plots of Mach number near the stall point (10 % span from the shroud): (a) Experiment (Strazisar et al. Reference Strazisar, Wood, Hathaway and Suder1989); (b) uncoupled blade; (c) coupled blade at Case B.
Table 2. Solution cases of static aeroelastic analysis.
The deformation solution of the coupled blade is presented in figure 8. It is important to emphasise that the discrepancy between Cases B and C is negligible, suggesting a convergence in the findings. It is apparent that, apart from undergoing radial stretching, the local blade section also undergoes bending and untwisting. This phenomenon becomes increasingly pronounced as the radius increases and leads to a non-uniform distribution of tip clearance. Case B is thoughtfully selected as the instructive exemplar. Despite the relatively modest magnitude (with a maximum deformation of 
$\varDelta (x,y,z) = (0.3302, 2.4984, -2.3240),\ \textrm {mm}$), the incorporation of this coupled solution adeptly resolves the concerns associated with the flow field around rigid blades. Figure 6(c) elucidates the expanded range and heightened intensity of a larger transonic flow region at the LE, whereas figure 7(c) portrays an augmented maximum Mach number.
Figure 8. Solution of static aeroelastic configurations of Cases A, B and C: (a) 
$z$–
$x$; (b) 
$y$–
$z$.
Furthermore, it is imperative to emphasise that numerous studies operate under the assumption that blades exhibit reduced and uniform tip clearances. Hence, the comparison with the simulation under the presumptions of 0.6-mm tip clearance (Zhang Reference Zhang2017) is depicted in figure 9. It becomes evident that the shape obtained through coupled simulation bears a closer resemblance to the experimental data, particularly in the proximity of peak efficiency. The present model has been validated and, thus, will subsequently serve as the foundation for the simulations in the following section.
Figure 9. Aerodynamic performance comparison: (a) total pressure ratio; (b) adiabatic efficiency.
3.2. Effect of ‘hot’ aerodynamic configuration
3.2.1. Static aeroelastic analysis
The variation in compressor performance can be attributed to the contrasting aerodynamic behaviours displayed by uncoupled and coupled blades. A fluid–structure coupling method has been applied to study the static aeroelastic problems of the ‘hot’ aerodynamic definition, where ‘hot’ means the updated blade shape undergoing design speed and temperature (Ramsey Reference Ramsey2006), whereas the rigid blades are referred to as a ‘cold’ shape. It is crucial to recognise that even minor deformations possess the potential to influence the overall structural performance, as they can affect the flow field. Consequently, we delve deeper into analysing the deviations arising from the ‘hot’ aerodynamic configuration. Figure 10 provides an illustration of the pressure distribution along chord-wise, whereas figure 11 delineates the corresponding interblade streamline diagram.
Figure 10. Surface pressure distribution of ‘cold’ and ‘hot’ shapes for: (a) Case A; (b) Case B; and (c) Case C.
Figure 11. Contours of the pressure and surface streamline of ‘cold’ and ‘hot’ shapes for: (a) Case A; (b) Case B; and (c) Case C.
Regarding Case B, an observable phenomenon is the heightened prominence of pressure distribution towards the LE in the ‘hot’ shape (figure 10b). Due to the slightly smaller clearance at the LE, fewer areas of tip vortex leakage are generated in the ‘hot’ shape from LE, allowing for a sustained flow velocity at the trailing edge (TE). In contrast, the tip vortex leakage in the ‘cold’ shape traverses the entire tip flow field until it reaches the TE of adjacent blades (figure 11b). Another distinction emerges in the distribution of shock waves along the blade surface. It is notable that the streamline on the ‘hot’ blade's surface illustrates a forward shift of the separation zone, signifying a significant reduction in the pressure distribution along the suction surface. These differences gradually decrease as the position approaches the root.
The same phenomenon exists in Case C (figures 10c and 11c). However, the larger difference appears at the 30 % span compared with the 10 % span from the shroud, which can be attributed to variations in the span-wise extent of the shock wave separation zone on the blade surface under different cases. In Case A, the pressure differential near the LE showcases a relatively subdued variation but a greater difference at the TE as shown in figure 10(a). This peculiarity is a consequence of the tip leaking vortex originating in closer proximity to the centre of the chord direction, as depicted in figure 11(a). The performance evaluation of a system based on a ‘cold’ shape can lead to inaccuracies due to the influence of ‘hot’ configurations. Since the operational point typically resides in the vicinity of peak efficiency, Case B, representing the fixed ‘hot’ shape, was selected for subsequent exploration.
3.2.2. Dynamic aeroelastic analysis
Temporal aeroelastic characteristics of aeroelasticity are scrutinised to discern disparities and influential factors between the ‘cold’ and ‘hot’ conditions. Quantitative evaluation of aerodynamic stability under various operational scenarios entails the analysis of displacement convergence rates, serving as a crucial indicator of damping magnitude within an oscillating system. The exponential function 
$y = a + b\textrm {e}^{ - \varsigma t}$ is employed to model the maximum value of the displacement response curve, with 
$\varsigma$ signifying the damping coefficient. The dynamic aeroelastic solution commences by applying minor perturbations to the ‘hot’ shape derived from the static aerodynamic results. Subsequently, the analysis progresses into the time domain to capture time-varying phenomena.
The maximum deformation occurs at the LE of the blade tip. Figure 12 illustrates the time-domain motion for the near-stall operating state and provides results of aeroelastic stability verification within the operational range. Significantly, compared with the ‘hot’ shape, the ‘cold’ shape exhibits smaller vibration amplitudes and faster decay rates (figure 12a), primarily attributed to its lower pressure levels. However, the stall point resulting from vibration is notably advanced across the entire operational range when compared with the behaviour observed in the ‘hot’ shape. Moreover, although a sudden drop in aeroelastic stability near stall is observed in both cases, the ‘cold’ shape experiences an untimely increase before reaching this condition (figure 12b).
Figure 12. Dynamic aeroelastic response: (a) time history for ‘cold’ and ‘hot’ shape near stall; (b) damping coefficient of the working area.
The distribution of displacement and equivalent stress at the extreme position of structural motion are depicted in figure 13, considering both the span-wise and chord-wise directions. Observing the deformation at the LE, it is notable that it experiences a rapid change with an increase in blade height. However, the change is less pronounced at lower positions but exhibits a significant increase after surpassing 40 % of the blade height at TE (figure 13a). When examining the equivalent strain energy and stress at the blade tip, where deformation is most significant, a stress concentration effect extending from the root to 70 % in the span-wise direction is observed (figure 13b,d). In terms of chord distribution, the stress level is notably higher at 30 % and 70 % than that at 10 % span from the shroud, with the peak value occurring between 1/3 and 2/3 from the LE (figure 13c).
Figure 13. Structural load at the maximum point: (a) deformation at the LE, 1/3 from the LE, 2/3 from the LE and the TE; (b) equivalent elastic strain at the LE, 1/3 from the LE, 2/3 from the the LE and the TE; (c) equivalent stress at 10 %, 30 % and 70 % span from the shroud; (d) equivalent stress at 1/3 and 2/3 from the LE of the pressure surface and the suction surface.
Based on the simulation, it can be concluded that there is no danger of flutter before the stall threshold within the design speed for NASA Rotor67. Furthermore, the presence of vibration induces stress concentration inside the blade structure, extending from the root to the middle section, as well as the central chord section. These phenomena elucidate the potential locations of fatigue-induced cracks.
3.3. Dynamic aeroelastic analysis of damaged blade
3.3.1. Effect of damage parameters
Representative damage parameters are selected for investigation based on the available results in figure 2. In the chord-wise distribution, the studied damage is limited to a distance of 
$1/3$ to 
$2/3$ of the LE span-wise damage starts at 
$30\,\%$, 
$60\,\%$ and 
$80\,\%$ span from the shroud and ends at the root, i.e. 
$h_d=0.3$, 
$0.6$ and 
$0.8$. Take 
$h_d=0.3$ as an example, damage extent 
$l_d$ is determined by dividing the remaining 70
$\,\%$ blade height into three parts, i.e. 
$l_d=0.3$, 0.5 and 0.7. Similarly, 
$h_d=0.6$ corresponds to 
$l_d=0.2$ and 
$0.4$, whereas 
$h_d=0.8$ corresponds to 
$l_d=0.2$. Finally, the damage level increases sequentially as 
$Sr=0.1$–
$1$ with 
$\Delta Sr=0.1$. Effects of damage level, extent and position on the aeroelastic behaviours and stability boundaries are discussed in detail. The near-stall working state is employed for simulation with a total pressure ratio of 
$P_{rc}=1.688$ and a non-dimensional flow rate of 0.94. All figures are plots of a typical point at the LE of the blade tip.
Figures 14 and 15 present the schematic diagram of span-wise damage and the damping coefficient 
$\zeta$ of above-mentioned cases for variable 
$Sr$, 
$h_d$ and 
$l_d$. In this study, 
$h_d = 0.3$, representing the minimum damage location studied, corresponds to a position closer to the blade tip. For 
$l_d=0.3$, the results indicate that the blade remains stable across all tested scenarios. In contrast to initial predictions, damaged blades with higher damage levels (lower 
$Sr$) exhibit a larger 
$\zeta$ within 
$Sr=0.2\sim 1$ upon examination of the effect of 
$Sr$. For 
$l_d=0.5$, 
$\zeta$ displays similar behaviour within the range of 
$Sr=0.4$–1 but experiences a sharp reduction from 
$Sr=0.4$ to 
$0.1$. The turning point, marking the transition from increasing to decreasing damping, is identified as the critical 
$Sr$ here. Upon further expansion to 
$l_d=0.7$, the critical 
$Sr$ increases to 
$Sr=0.6$. In the above cases with 
$h_d=0.3$, the damaged blade with 
$Sr=0.2$, 
$l_d=0.5$ and 
$Sr=0.1-0.3$, 
$l_d=0.7$ will experience aeroelastic instability. Under conditions of 
$h_d = 0.6$ and 
$0.8$, 
$\zeta$ of damaged blades are generally lower than those of healthy blades (
$Sr=1$), yet the aforementioned anomaly persists. In addition, as 
$l_d$ increases, the critical 
$Sr$ also gradually increases under the cases with 
$h_d=0.3$ and 
$0.6$. Before the critical 
$Sr$ threshold, the optimal 
$\zeta$ is observed under the condition of 
$h_d=0.3$, 
$l_d=0.5$. It is evident that before reaching the critical 
$Sr$ threshold, damage regions that exhibit superior stability tend to experience a more significant decrease in stability after surpassing the critical 
$Sr$. Once the 
$Sr$ threshold is surpassed, 
$\zeta$ undergoes a sharp decline, ultimately leading to the emergence of flutter.
Figure 14. Damping coefficient 
$\zeta$ for varying 
$S=0.1$–
$1$, 
$l_d=0.3$, 0.5 and 
$0.7$ with 
$h_d=0.3$ fixed.
Figure 15. Damping coefficient 
$\zeta$ for varying 
$S=0.1$–
$1$, 
$l_d=0.2$ and 
$0.4$ with 
$h_d=0.6$ fixed and 
$l_d=0.2$ with 
$h_d=0.8$.
To elucidate the effects of varying damage positions, a comparative analysis is conducted between the cases where 
$h_d = 0.6$, 
$l_d = 0.2$ and 
$h_d = 0.8$, 
$l_d = 0.2$. Specifically, the value of 
$h_d = 0.6$ signifies a location closer to the span-wise midpoint, whereas 
$h_d = 0.8$ indicates a position nearer to the root. The results show that damages occurring for 
$h_d = 0.6$ exhibit a more gradual change rate, lower damping levels and a smaller critical 
$Sr$ before reaching the critical 
$Sr$. Thus, the combined effects of 
$Sr$, 
$h_d$ and 
$l_d$ reveal that damage occurring nearer the blade root poses a greater risk of instability, as evidenced by comparing the longitudinal distribution of the damaged blocks in the span-wise direction. Hence, the observed difference in the rate of change of the damaged blade's vibrational damping 
$\zeta$ with 
$Sr$ can be primarily attributed to the combined effect of the location 
$l_d$ and extent 
$h_d$ of the damage incurred.
To sum up, the increase of damage level, as observed within a certain range can fortunately enhance aeroelastic stability. However, once it surpasses the critical 
$Sr$, the aeroelastic stability presents a rapid decrease or even results in instability. Moreover, the study demonstrates that the location and extent of damage mainly affect the aeroelastic stability by influencing the critical 
$Sr$ value and the change rate of 
$\zeta$ with respect to 
$Sr$. Damage closer to the blade tip enhances stability, whereas damage nearer the blade root increases the risk of instability. This interesting finding could be utilised for weight optimisation of the rotor blade by adjusting an appropriate local material.
3.3.2. Exploration of blade flutter mechanism
In this section, the investigation into the distinctions between aeroelastic stability and instability conditions of damaged blades with 
$Sr=0.2$, 
$l_d=0.7$ and 
$h_d=0.3$ will be conducted to explore the potential flutter mechanism for the blade. The response time history is depicted in figure 16 for varying total pressure ratios 
$P_{rc}=1.668$, 1.675, 1.682 and 1.688. It is evident that with an increase in 
$P_{rc}$, aeroelastic stability decreases, with a relatively unchanged vibration frequency. Figure 17 illustrates the aerodynamic forces in the 
$x$ (span-wise), 
$y$ (circumferential-wise) and 
$z$ (chord-wise) directions (
$F_{x, y, z}$). Their amplitudes consistently vary with displacement magnitude, but much more notable differences are observed in the phase portraits of these forces. In the cases investigated, the phase of 
$F_x$ exhibits a generally consistent pattern, with notable variations observed in its magnitude. When considering the variable 
$F_y$, a little phase advancement is seen as the value of 
$P_{rc}$ increases. Significantly, the phase progression of 
$F_z$ is particularly prominent. From the perspective of energy transfer, the phase discrepancies between unstable pressure and blade vibration result in variations in the work of the aerodynamic force on the blade, which have a negative effect on the aeroelastic stability at 
$P_{rc}=1.688$.
Figure 16. Dynamic aeroelastic for varying 
$P_{rc}=1.668$–
$1.688$ with 
$Sr=0.2$, 
$h_d=0.3$ and 
$l_d=0.7$ fixed: (a) displacement in time history; (b) damping coefficient.
Figure 17. Aerodynamic force for varying 
$P_{rc}=1.668$–
$1.688$ with 
$Sr=0.2$, 
$h_d=0.3$ and 
$l_d=0.7$ fixed.
To understand the mechanism behind the observed phase shift, a detailed examination on the flow field is necessary. The most significant phase change occurs in the 
$z$-direction, prompting us to focus on the flow field information at the blade's chord. Figures 18 and 19 present contour plots depicting the pressure and Mach number distributions in the flow field between the vibrating blades, respectively. At time 
$t=0.0200\ \textrm {s}$, the displacements for the three conditions exhibit minimal differences. However, the flow field reveals substantial variations in the high-pressure region at the TE, which extends directly to the pressure surface of the neighbouring blade as displacement progresses, particularly notable at 
$t=0.0205\ \textrm {s}$, under 
$P_{rc}=1.688$. This phenomenon leads to a heightened pressure difference between the blade's upper and lower surfaces, subsequently resulting in an increased displacement at 
$t=0.0210 \textrm {s}$. In particular, distinctions in the flow field around the blades’ LE become evident at this juncture, characterised by the appearance of a smaller low-pressure zone for the 
$P_{rc}=1.688$ case, which is absent in the other two cases.
Figure 18. Pressure distribution at 
$t=0.0200$, 
$0.0205$, 
$0.0210$, 
$0.0215$ and 
$0.0220$ s of 10 % span from the shroud.
Figure 19. Mach number distribution at 
$t=0.0200$, 
$0.0205$, 
$0.0210$, 
$0.0215$ and 
$0.0220$ s and vorticity distribution at 
$t=0.0225$ and 
$0.0230$ s of 10 % span from the shroud.
From 
$t=0.0210$ to 
$0.0215\ \textrm {s}$, the low-pressure conditions near the LE exhibit a modest expansion, whereas there is a notable reduction in pressure conditions near the TE. By 
$t=0.0220\ \textrm {s}$, the three displacements had converged. However, significant changes in the TE low-pressure conditions for 
$P_{rc}=1.668$ were detected, with a high-pressure regime in the blade's middle impacting nearby blades. These phenomena are reflected in the Mach number and vorticity distribution, as illustrated in figure 19, corresponding to a subsonic regime at the TE due to separation and a transonic regime at the LE In the vorticity plots, it is evident that from 
$t=0.0025$ to 
$0.0023\ \textrm {s}$, the higher-intensity shock wave at the LE under the 
$P_{rc}=1.688$ condition significantly affects convergence. Simultaneously, towards the TE, there is a larger region of vortex detachment within the flow field between the blades, intensifying its impact on adjacent blades.
Upon observing the phenomena, it is noteworthy that the TE subsonic region (high-pressure region) stemming from boundary layer separation exhibits a dual functionality. Under stable conditions, this region serves to dampen vibrations. Conversely, during fluttering, as the high-pressure region expands onto adjacent blade surfaces, it paradoxically acts as a stimulus for vibrations, ultimately leading to instability. Thus, the distributions of the LE and TE lead to differences in the flow field divisions within the blade's midsection. Utilising vorticity plots to characterise this phenomenon can effectively elucidate the disparities between stability and instability. The aerodynamic instability condition at the LE is characterised by a more pronounced shock wave that converges towards the LE. Moreover, there are vortices with a wide range and high intensity located between the blades near the TE, contributing to instability and uneven energy distribution within the flow field. With the foregoing explanations, it is possible to conduct a preliminary analysis of the physical causes behind aeroelastic instability from two perspectives: aerodynamic work and flow field characteristics. The analysis serves to unveil the instability mechanisms and offers guidance for subsequent research.
3.3.3. SPOD modal identification of flutter mechanism
In a previous study (see Ji et al. Reference Ji, Xie, Zhang and Maqsood2023), some of the present authors clarified the ability of the SPOD method to extract global modes and their spectral advantages, which will not be further elaborated on here. Thus, as shown in § 3.3.2, the divergent and convergent cases at 
$P_{rc}=1.688$ and 
$P_{rc}=1.675$ are chosen to investigate the special characteristics via the SPOD method. A total of 500 snapshots of pressure distribution were collected in the time domain of the three-dimensional flow field of the damaged blade with 
$Sr=0.2$, 
$l_d=0.7$ and 
$h_d=0.3$. Each block contains 64 snapshot data, with adjacent block matrices overlapping by 50 %.
The SPOD method decomposes the data into a set of orthogonal modes, each associated with a specific frequency and spatial structure. The resulting SPOD spectra (eigenvalues 
$\lambda$) are plotted vs (non-dimensional) frequency in figure 20. For each case, three frequencies are selected for characteristic SPOMs. Here 
$f_{1}$–
$f_{3}$ are 0.1328BPF (blade passing frequency) and its multiples. Both cases are dominated by the first-order mode, whereas higher frequencies and higher-order modes serve as complementary factors, enhancing and refining the effect. In terms of frequency, a small decrement in the energy levels associated with 
$f_{1}$–
$f_{3}$ is observed for the flutter case 
$P_{rc}=1.688$. However, 
$f_{3}$ from stable case 
$P_{rc}=1.675$ no longer serves as the primary distinguishing characteristic, because of its energy reduction of multiple orders of magnitude. To facilitate a comparative analysis, a representation of 
$f_{3}$ for 
$P_{rc}=1.675$ is also included. Figures 21 and 23 show the first two SPOMs at 
$f_{1}$–
$f_{3}$ for 
$P_{rc}=1.688$ and 
$P_{rc}=1.675$ of 5 % span from the shroud. The span-wise stream surface is determined by the periodic boundary between the blades, shown in figures 22 and 24. Due to the representative changes near the blade tip, the chord position is given by 5 % span from the shroud.
Figure 20. SPOM energy spectrum: (a) 
$P_{rc}=1.688$; (b) 
$P_{rc}=1.675$.
Figure 21. The first- and second-order SPOMs of 5 % span from the shroud at 
$P_{rc}=1.688$: (a) first order at 
$f_{1}$; (b) first order at 
$f_{2}$; (c) first order at 
$f_{3}$; (d) second order at 
$f_{1}$; (e) second order at 
$f_{2}$; ( f) second order at 
$f_{3}$.
Figure 22. The first- and second-order SPOMs of span-wise stream surface at 
$P_{rc}=1.688$: (a) first order at 
$f_{1}$; (b) first order at 
$f_{2}$; (c) first order at 
$f_{3}$; (d) second order at 
$f_{1}$; (e) second order at 
$f_{2}$; ( f) second order at 
$f_{3}$.
In the span-wise direction of 
$f_{1}$-
$1\textrm {st}$ for two cases, it is observed that the leakage vortex occurs at the tip gap near the LE of the blade (figures 22a and 24a), which interferes with the shock wave and generates disturbance chord-wise (figures 21a and 23a). Thus, the 
$f_{1}$-
$2\textrm {nd}$ reveals the presence of small-scale disturbance clusters downstream of the shockwave chord-wise (figures 21d and 23d). However, these disturbance clusters dissipate rapidly within the flow passage as they propagate downstream. Moreover, 
$f_{1}$-
$2\textrm {nd}$ for the flutter case 
$P_{rc}=1.688$ indicates that the extent of the tip leakage flow results in a more extensive influence on the blade tip region. For 
$f_{2}$-
$1\textrm {st}$ (figures 21b and 23b), the SPOMs show that the interaction between the tip clearance leakage vortices and the mainstream flow results in the formation of large-scale disturbance. Such disturbance propagates downstream along the chord-wise direction from the LE and concurrently spreads from the tip to the root of the blade (figures 22b and 24b). These occurrences in both flutter and stable cases suggest a potential correlation with the unsteady periodic fluctuations in the blade flow field.
Figure 23. The first- and second-order SPOMs of 5 % span from the shroud at 
$P_{rc}=1.675$: (a) first order at 
$f_{1}$; (b) first order at 
$f_{2}$; (c) first order at 
$f_{3}$; (d) second order at 
$f_{1}$; (e) second order at 
$f_{2}$; ( f) second order at 
$f_{3}$.
Figure 24. The first- and second-order SPOMs of span-wise stream surface at 
$P_{rc}=1.675$: (a) first order at 
$f_{1}$; (b) first order at 
$f_{2}$; (c) first order at 
$f_{3}$; (d) second order at 
$f_{1}$; (e) second order at 
$f_{2}$; ( f) second order at 
$f_{3}$.
Different modes appear in 
$f_{2}$-
$2\textrm {nd}$ and 
$f_{3}$. For the stable case, 
$f_{2}$-
$2\textrm {nd}$ (figure 23e) has a phase difference in disturbance clusters compared with 
$f_{1}$-
$1\textrm {st}$. However, 
$f_{2}$-
$2\textrm {nd}$ becomes even more intricate for flutter (figures 21e and 22e). This complexity is evident in the increased quantity of disturbance clusters. Despite being smaller scale compared with 
$f_{2}$-
$1\textrm {st}$, these disturbances maintain their intensity along both chord-wise and span-wise flow channels and engage in interactions among the blades. As expected, the 
$f_{3}$ for 
$P_{rc}=1.675$ (figures 23c, f and 24c, f) exhibits a flow field structure similar to that of 
$f_{1}$ due to its lower energy. The phenomenon distribution of 
$f_{3}$ (figure 21c, f) resembles that of 
$f_{2}$ for 
$P_{rc}=1.688$, yet the influence scope of the 
$f_{3}$ disturbance exhibits a backwards shift. These disturbance groups exist through the entire interblade channel in the chord direction, but no longer maintain their strength in the span-wise direction.
Overall, the SPOMs observed between the flutter and stable cases can be primarily summarised into two crucial factors. First, the high-frequency perturbation with a larger scale exhibits a distinct frequency characteristic that is intrinsic to the SPOD method. Second, the high-order small-scale perturbation clusters, which propagate across the flow channels, maintain their strength over a broader area, thus exerting a more widespread influence between blades.
4. Conclusions
This study was mainly conducted to predict the flutter boundary of transonic rotor blades such as the representative NASA Rotor67 with possible structural damage, which is parameterised with varying damage levels, extents and locations. Furthermore, the physical mechanism of rotor blade flutter has been explored in time and frequency domains. Specifically, a high-fidelity CFD/CSD coupling model in the time domain has been established for both static and dynamic aeroelastic analysis at a predetermined rotational speed (16 043 rpm), and SPOMs in the frequency domain have been identified for the evaluation of factors inducing fluttering instability, respectively.
The main conclusions are as follows.
(i) The ‘hot’ configuration results in more accurate steady-state and aeroelastic results compared with the ‘cold’ shape. The differences are ascribed to variations in the flow field vortices, which are induced by the variable clearance between the blade tips.
(ii) Curves of damping coefficient with damage level indicate that there is a critical threshold, below which aeroelastic stability is improved positively with increased damage level (
$Sr$ decreases accordingly), but inversely, a significant loss of stability is induced. The damage extent changes the critical $Sr$ value and the damage location affects the change rate of the curve in terms of the damping coefficient with respect to $Sr$. In particular, damage located near the blade tip can enhance the aeroelastic stability, whereas the damage near the blade root will more likely induce instability.(iii) The numerical results on flutter analysis revealed that the boundary layer separation played inverse roles in stable and flutter conditions. In stable conditions, the high-pressure subsonic region induced from boundary layer separation provides damping and prevents vibrations from spreading. However, during flutter, the high-pressure region diffuses to adjacent blade surfaces, stimulating vibrations and finally triggering instability.
(iv) The exploration of the flutter mechanism via the SPOMs reveals that the involvement of high-frequency large-scale disturbances are the key factor inducing flutter, and the high-order small-scale disturbances with slow dissipation further stimulate this phenomenon. These disturbances maintain their intensity and interact with each other, persistently causing extensive influence both chord-wise and span-wise.
In summary, by investigating the hot configuration, aeroelastically damaged blade and the flutter mechanism, this study has provided important insights into the aeroelastic behaviour of a transonic rotor blade. These results offer valuable information that can be used to guide the design of blade structures to reduce flutter risk and forecast rotor blade damage.
Acknowledgements
The authors gratefully acknowledge the technical support provided by the Multidisciplinary Flight Dynamics and Control Laboratory of the School of Astronautics, Northwestern Polytechnical University. The authors would also like to thank the Aeroelasticity Research Group of the Department of Mechanical Engineering and Materials Science, School of Engineering, Duke University for invaluable suggestions, insightful discussions and mentorship throughout this study.
Funding
This study was supported by the NSFC International (Regional) Cooperation and Exchanges Program (12211540709) and the National Natural Science Foundation of China (11972294).
Declaration of interests
The authors report no conflict of interest.
