To identify the spatiotemporal coherent structure of compressor tip leakage flow, spectral proper orthogonal decomposition (SPOD) is performed on the near-tip flow field and the blade surface pressure of a low-speed compressor rotor. The data used for the SPOD analysis are obtained by delayed-detached eddy simulation, which is validated against the experimental data. The investigated rotor near-tip flow field is governed by two tip leakage vortices (TLV), and the near-tip compressor passage can be divided into four zones: the formation of main TLV (Zone I), the main TLV breakdown (Zone II), the formation of tip blockage cell (Zone III), and the formation of secondary TLV (Zone IV). Modal analysis from SPOD shows that a major part of total disturbance energy comes from the main TLV oscillating mode in Zone I and the main TLV vortex shedding mode in Zone III, both of which are low-frequency and low-rank; on the contrary, modal components in Zones II and IV are broadband and non-low-rank. Unsteady blade forces are mainly generated by the impingement of the main TLV on the blade pressure surface in Zone III, rather than the detachment of the secondary TLV from the blade suction surface in Zone IV. These identified coherent structures provide valuable knowledge for the aerodynamic/aeroelastic effects, turbulence modeling, and reduced-order modeling of compressor tip leakage flow.

## I. INTRODUCTION

Compressor tip leakage flow has long been known as the dominant mechanism in determining compressor work input, efficiency and stability.^{1} Recent numerical and experimental studies have highlighted the significance of the TLV in the stall inception^{2} and the nonsynchronous vibration^{3} of compressors. Due to the multi-scale nature of the TLV, flow visualization experiments^{3–5} and scale-resolving simulations^{6–8} have been conducted to investigate its flow physics. Flow field data in terabyte-scale are generated, but have these data been exploited effectively?

Modal analysis is a data mining tool that extracts the general spatial and temporal flow structures (i.e., modes) from the snapshots of unsteady flow field data. The identified modes not only help advance the understanding of the flow physics but also serve as the basis for reduced-order modeling of the unsteady flow field.^{9,10}

To date, the most popular mode decomposition methods for compressor flows are proper orthogonal decomposition (POD) proposed by Lumley^{11} and dynamic mode decomposition (DMD) proposed by Schmidt.^{12} Examples include but not limited to POD of compressor wake,^{13} POD and DMD of compressor surge,^{14} DMD of corner separation,^{15} and DMD of spike-type rotating stall.^{16} POD extracts coherent structures in turbulent flows by identifying the optimal set of orthogonal modes ranked by their modal energy. However, the physical meaning of each mode may be difficult to interpret, because the mode shape involves flow field information at different frequencies. DMD extracts the frequency-resolved flow structure with each mode having a characteristic frequency of oscillation and a growth/decay rate. However, DMD lacks a general mode ranking criteria, making it challenging to determine the dominant modes. From the perspective of a turbomachinery aerodynamicist, an ideal mode decomposition method should combine the advantages of both POD and DMD: the identified modes can be ranked according to their energy level with each mode having a characteristic frequency.

Spectral proper orthogonal decomposition (SPOD) is an ideal mode decomposition method for this purpose, whose modes are frequency-resolved and energy-ranked. SPOD was originally proposed by Lumley^{17} and recently revived by Towne *et al.*^{18} It has been successfully applied to extract the spatialtemporal modes of jet and wind turbine flows.^{19–21} The motivation of this paper is to unveil the modal component of compressor tip leakage flow by using SPOD, which will be valuable for understandings of the aerodynamic/aeroelastic effects, turbulence modeling, and reduced-order modeling of compressor tip leakage flow. In the following, the methodologies including the flow solver, the SPOD method, and the case studied will be detailed first. Afterward, the SPOD analysis of the near-tip flow field and the blade surface pressure will be presented and discussed.

## II. METHODOLOGY

### A. Flow solver

The data analyzed in this study are obtained by using the in-house hybrid RANS/LES flow solver HADES. It uses a node-centered finite volume scheme with an edge-based data structure. The numerical fluxes are solved by the Jameson–Schmidt–Turkel scheme^{22} with a Roe's matrix,^{23} leading to the second-order accuracy in space except for the vicinity of a shock. A hybrid low-dissipation scheme^{24} is adopted in the DDES mode of the solver, which reduces the fourth-order dissipation term in the LES region. The second-order implicit time integration scheme is applied for temporal discretization. For simulation of turbulent flows, the standard Spalart–Allmaras turbulence model^{25} and its variants^{26} are available for the RANS branch, and the standard DDES method^{27} and its upgrades^{28} are available for the hybrid RANS/LES branch. Validation and verification of HADES can be found in previous works.^{26,28}

In this work, a variation of the SA-based DDES method^{28} is used. Compared to the original formulation,^{27} an alternative grid spacing is adopted as described in Eq. (1),

where Δ_{max} is the longest edge length of the cell and is used in the original DDES formulation;^{27} Δ_{vol} is the cube root of the cell volume and is widely used in LES;^{29} the shielding function *f _{d}* and the threshold $fd0=0.99$ are used to switch the grid spacing definition between the RANS branch ($fd<fd0$) and the LES branch ($fd\u2265fd0$);

^{30}

^{,}

*F*is an empirical term based on the vortex tilting measure (VTM) to unlock the physics of Kelvin–Helmholtz instability,

_{KH}^{31}as detailed in Eq. (2),

where the constants are $FKHmax=1.0,\u2009FKHmin=0.1,\u2009a1=0.15$, and $a2=0.3$.

The present DDES method achieves LES resolution in far-wall regions but uses RANS to model the near-wall boundary layer flows. Hence, DDES is more numerically efficient but less accurate (mainly in boundary layers) compared to the wall-resolved LES method. Considering compressor tip leakage flows are featured by flow separations due to the sudden change of geometry at the blade tip (i.e., similar to a backward-facing step flow), the advantage of wall-resolved LES in boundary layers is not essential. Therefore, DDES instead of wall-resolved LES is adopted in this work.

### B. Spectral proper orthogonal decomposition

SPOD is a frequency-domain variant of POD^{11} that is designed for statistically stationary flows. It can capture the energy-ranked and mutually orthogonal modes at each frequency. Mathematically, the SPOD modes are eigenvectors of the cross-spectral density (CSD) matrix at each frequency, and the eigenvalues correspond to the energy of each mode at a distinct frequency. Detailed theory and derivation of the SPOD method can be found in Refs. 18 and 21. The Python script of SPOD written for this work is open-accessed [Package available from: https://github.com/HexFluid/spod\_python (accessed June 2021)]. Key procedures of SPOD are illustrated in Fig. 1. Each procedure is detailed as follows.

- A matrix containing the spatiotemporal flow data is constructed first. Let the vector $qk\u2208RNq$ represent the
*k*th time snapshot after subtracting the time-averaged data, and its length $Nq=NgNv$ where*N*and_{g}*N*are the number of grid points and investigated variables, respectively. By assembling the snapshot flow field data in chronological order, the spatiotemporal data matrix is obtained,_{v}(3)$Q=[q1,q2,\u2026,qNt]\u2208RNq\xd7Nt,$where

*N*is the total number of snapshots._{t} - The data matrix
is then segmented into*Q**N*($Nb\u226aNt$) blocks by applying the Welch periodogram method,_{b}^{32}with each block consisting of*N*snapshots. Each block is assumed to be a statistically independent realization under the ergodicity hypothesis. The discrete Fourier transform (DFT) is performed on each block to convert the problem to the frequency domain. To avoid spectral leakage, each block is windowed by the Hamming window function_{f}^{33}and is overlapped with neighboring blocks. The resultant matrix for the*j*block is_{th}(4)$Q\u0302(j)=[q\u03021(j),q\u03022(j),\u2026,q\u0302Nf(j)]\u2208CNq\xd7Nf,$where

*N*is equivalent to the number of resolved frequencies._{f} - Next, the frequency-domain matrices of blocks are reshaped according to the frequency. For the
*k*_{th}frequency, the corresponding matrix is(5)$Q\u0302k=[q\u0302k(1),q\u0302k(2),\u2026,q\u0302k(Nb)]\u2208CNq\xd7Nb,$The weighted CSD matrix for the*k*_{th}frequency is, therefore, obtained as(6)$Sk=1NbW12Q\u0302k*Q\u0302kW*12\u2208CNb\xd7Nb$whereis the weight matrix for scaling of different flow variables. The definition of*W*determines the physical meaning of the SPOD mode energy. This work follows the derivation of total disturbance energy,*W*^{34}(7)$E=12\u222bV(RT\xaf\rho \xaf\rho \u20322+R\rho \xaf(\gamma \u22121)T\xafT\u20322+\rho \xafui\u20322)\u2009dV,=12\u222bV(1\gamma p\xafp\u20322+(\gamma \u22121)p\xaf\gamma R2\Delta s\u20322+\rho \xafui\u20322)\u2009dV,$where $\Delta s$ is the entropy rise,

*V*is the control volume of the grid cell,*R*is the gas constant,*γ*is the specific heat ratio, $(\xb7)\u2032$ denotes the unsteady fluctuation of a quantity, and $(\xb7)\xaf$ denotes the time-average of a quantity. In these equations, the $\rho \xafui\u20322/2$ term corresponds to the turbulence kinetic energy (TKE), and the other terms can be interpreted as potential energy due to fluctuations of density, temperature, pressure and entropy. For flows where the temperature fluctuation is of most interest (e.g., heat transfer in turbines), the quantities of interest $q=[\rho ,ux,ur,ut,T]T$ and the diagonal weight matrix $Wq=diag[RT\xaf\rho \xaf,\rho \xaf,\rho \xaf,\rho \xaf,R\rho \xafT\xaf(\gamma \u22121)]$ are recommended; for flows where the pressure fluctuation is of most interest (e.g., compressor aeroelasticity), $q=[p,ux,ur,ut,\Delta s]T$ and $Wq=diag[1\gamma p\xaf,\rho \xaf,\rho \xaf,\rho \xaf,(\gamma \u22121)p\xaf\gamma R2]$ are preferred. - Finally, the eigenvalue decomposition is performed on the weighted CSD matrix $Sk$ for each frequency,(8)$Sk=\Psi k\Lambda k\Psi k*,\Phi k=Q\u0302k\Psi k\Lambda k\u22121/2,\Phi k=[\varphi k(1),\varphi k(2),\u2026,\varphi k(Nb)]\u2208CNq\xd7Nb,\Lambda k=diag(\lambda k(1),\lambda k(2),\u2026,\lambda k(Nb))\u2208CNb\xd7Nb,$
The matrix $\Phi k$ calculated from the eigenvectors represents the SPOD mode shapes at the

*k*frequency. These modes are ranked according to their corresponding eigenvalues $\Lambda k$, which can be physically interpreted as the total disturbance energy of the mode in the investigated flow domain._{th}

### C. Case description

The axial compressor investigated in this paper is the BUAA Stage B rotor from Du *et al.,*^{4} whose specifications are summarized in Table I. The rotor flow has been investigated experimentally using five-hole probe^{35,36} and stereoscopic particle image velocimetry (SPIV)^{4} measurements under the stage environment. The measurement locations of SPIV are illustrated in Fig. 2(a), including 10 slices between 10% and 100% of blade tip chord near the SS and 11 slices between 5% and 105% of blade tip chord near the PS.

Parameter . | Symbol . | Value . |
---|---|---|

Hub to tip radius ratio | $rh/rt$ | 0.6 |

Tip gap to chord ratio | $\tau /ct$ | 1.8% |

Tip aspect ratio | $h/ct$ | 1.2 |

Flow coefficient | $ux1/Um$ | 0.58 |

Loading coefficient | $\Delta pt/\rho 1Ut2$ | 0.52 |

Tip Mach number | $M1t$ | 0.17 |

Reynolds number | $Re\tau $ | $1.4\xd7104$ |

Number of blade | Z | 17 |

Parameter . | Symbol . | Value . |
---|---|---|

Hub to tip radius ratio | $rh/rt$ | 0.6 |

Tip gap to chord ratio | $\tau /ct$ | 1.8% |

Tip aspect ratio | $h/ct$ | 1.2 |

Flow coefficient | $ux1/Um$ | 0.58 |

Loading coefficient | $\Delta pt/\rho 1Ut2$ | 0.52 |

Tip Mach number | $M1t$ | 0.17 |

Reynolds number | $Re\tau $ | $1.4\xd7104$ |

Number of blade | Z | 17 |

To simplify the numerical model, only one blade passage of the rotor is considered. The flow domain and the boundary conditions are illustrated in Fig. 2(b). At the inlet, the measured radial profile of total pressure, total temperature, and velocity directions are prescribed. Since the measured freestream TKE is two orders of magnitude smaller than those generated by the Kelvin–Helmholtz instability in the main TLV core,^{28} no additional turbulent perturbations are given at the inlet boundary. A choked nozzle^{37} is attached to the rotor exit, whose throat area is adjusted to change the compressor operating condition. This work focuses on the peak efficiency (PE) and the near-stall (NS) conditions at the design speed. Euler equations are solved in the nozzle domain to reduce the calculation cost. A mixing plane is used at the interface between the rotor domain and the nozzle domain.

The rotor domain is meshed by a hexahedral grid with 3.7 × 10^{6} points, which is generated by AutoGrid.^{38} The number of mesh points on the spanwise direction, pitchwise direction, and blade surfaces is 117, 97, and 305, respectively. In order to resolve the tip leakage flow accurately, over 50% of the mesh points are located above the 80% span, including 29 layers inside the tip gap. The average *y*^{+} value for the first-layer mesh points is 1.1, and the linear-law branch of the standard wall function^{39} is used to calculate the shear stress at those points. A zoom-in view of the rotor mesh is presented in Fig. 2(c). The nozzle domain is meshed by a hexahedral grid with 0.2 × 10^{6} points, which is generated by an in-house Python script. The non-dimensional time step is set as $\Delta t\xb7BPF=1/300$. This time step corresponds to a local CFL number of $\Delta t\xb7Ut/\Delta \theta =1/3$, which is finer than the recommended value of unity^{40} and hence achieves good time accuracy. Within each time step, 30 Newton iterations are performed so that the mean residual drops over two orders of magnitude. The converged steady RANS simulation results were used as the initial condition for the DDES simulations. After around two rotor revolutions, the DDES solutions reached a statistically steady state. The simulations were then performed for another five rotor revolutions in which all snapshots were collected for time-averaging and SPOD analysis.

The grid density has a significant effect on the accuracy of scale-resolving simulations. To quantify the grid-induced error, four sets of grids are investigated. These grids have 0.5, 1.2, 3.7, and 24.9 × 10^{6} grid points, respectively, and the grid refinement is performed uniformly in all dimensions. The predicted rotor static pressure rise coefficient $Cp=\Delta p/12\rho Um2$ is plotted against the nominal grid spacing $\Delta =1/Ng1/3$ in Fig. 3(a). The error bars of the DDES results are based on the root mean square of *C _{p}*, which represents the unsteadiness of the flow. Because both the discretization error of a second-order accurate solver and the turbulence modeling error of a Smagorinsky-type subgrid-scale (SGS) model scale with $\Delta 2$, a parabola curve $Cp=k\Delta 2+Cp,ideal$ can be fitted through the DDES results, and its intercept $Cp,ideal$ at Δ = 0 denotes the ideal prediction free of grid-induced errors. The results show that $Cp,ideal$ falls within the measured

*C*range, and the grid used in this work with 3.7 × 10

_{p}^{6}grid points underpredicts

*C*by 0.03. The relative error induced by grid can be defined by $|Cp/Cp,ideal\u22121|$, which is plotted against the total number of grid

_{p}*N*in Fig. 3(b). The grid-induced error of

_{g}*C*for the grid used in this work is about 4.8%, which is sufficiently small for capturing the TLV structures as will be shown later in Sec. III B. Using a finer grid, a high-order numerical scheme and a more accurate SGS model can reduce the grid-induced error further, but this is out of the scope of the paper.

_{p}The accuracy of time-averaging improves with the total number of snapshots *N _{t}* and the sampling frequency

*f*, but using a high value of

_{s}*N*and

_{t}*f*is computationally expensive. In searching for an optimal set of

_{s}*N*and

_{t}*f*considering the trade-off between accuracy and cost, a parametric study has been conducted as presented in Fig. 4. In this plot, the quantities of interest include the rotor static pressure rise

_{s}*C*, the normalized axial velocity

_{p}*u*and the normalized TKE, where the latter two quantities are obtained from six numerical probes evenly distributed between 0% and 100% chord at 90% span and 50% pitch; the relative error is defined by $|q/qref\u22121|$, where

_{x}*q*is the most accurate time-averaged results in this work that is obtained by using the baseline parameters of

_{ref}*N*= 5 and $fs/BPF=300$. The results indicate that using

_{t}*N*= 2 and $fs/BPF=150$ is most cost-effective, which leads to a relative error in time-averaging below 5% for all the quantities of interest. For better accuracy, the baseline parameters of

_{t}*N*= 5 and $fs/BPF=300$ are used for time-averaging throughout the paper.

_{t}## III. FLOW FIELD VISUALIZATION

### A. Transient flow field

The transient flow field of the investigated rotor is illustrated in Fig. 5, where the vortical structures are identified by the iso-surfaces of Q-criterion. At the PE condition, a large-scale TLV is formed immediately after the leading edge due to the shear effect between the tip leakage flow and the incoming main flow. Afterward, the spiral-type vortex breakdown occurs on the main TLV near the mid-pitch, which generates smaller eddies that impinge on the PS of the adjacent blade. Part of the main TLV flow passes through the tip gap of the adjacent blade and forms a secondary TLV, which detaches from the mid-chord SS and propagates toward the rotor exit. The impingement and the detachment of the TLVs exert unsteady forces on the blade. At the NS condition, a similar double-leakage structure can be observed. The difference is that the main TLV breakdown and therefore the formation of the secondary TLV move upstream. Part of the eddies originated from the main TLV breakdown spill over the leading edge of the adjacent blade, which can potentially initiate a spike-type rotating stall cell.^{41} The evolution of the double-leakage structure found in this work is consistent with the previous research.^{7}

To check the quality of the DDES solution, the resolution quality index $IQ\nu $^{42} is plotted on the iso-surfaces of Q-criterion in Fig. 5. The higher the $IQ\nu $, the finer the resolution. The results show that the $IQ\nu $ value near the main and the secondary TLV region is larger than 0.80, indicating a good resolution of DDES. To verify the DDES resolution further, the power spectrum density (PSD) of the fluctuation component of axial velocity is examined in Fig. 6. Data presented in this plot come from the numerical probes at 90% span and 50% pitch. It shows that DDES can resolve the inertial subrange of the turbulence spectrum down to the frequency of 10 BPF, which again demonstrates the good quality of the DDES result. Another observation is that the velocity fluctuation at the NS condition is stronger than the PE condition at 20% and 40% chord length, which is due to the more upstream location of the main TLV breakdown.

### B. Time-averaged flow field

The time-averaged chordwise vorticity fields from the experiment and DDES are compared in Fig. 7(a). It is shown that DDES captures qualitatively the same trajectory of the main TLV core as the experiment does at both the PE and NS conditions. After the breakdown of the main TLV, the vorticity dissipates rapidly. A quantitative comparison on the trajectory of the main TLV core is presented in Fig. 7(b). In this plot, the location of the TLV core is defined by the location of the maximum vorticity magnitude; the shaded area represents the pitchwise grid spacing, which indicates the uncertainty of the DDES results. It is shown that DDES can capture the main TLV trajectory with sufficient accuracy.

The time-averaged chordwise velocity fields from the experiment and DDES are compared in Fig. 8(a). Downstream of the main TLV breakdown, a tip blockage cell is formed at the casing-PS corner for both the PE and NS conditions. The tip blockage cell is composed of low-momentum fluid with small streamwise velocity near the casing wall. It reduces the effective passage area of the compressor, hence given its name “blockage cell.” It is shown that DDES predicts qualitatively the same tip blockage cell at both conditions. To quantify the blockage size, the volume enveloped by the iso-surface of a threshold velocity value is calculated, and the total volume of the compressor passage is used for normalization. The normalized blockage size vs the threshold velocity value is presented in Fig. 8(b). It shows that DDES can predict the tip blockage size with sufficient accuracy.

The entropy field is important for understanding the compressor loss of efficiency. Following the derivation of Denton,^{43} the entropy field can be non-dimensionalized by the compressor work input, yielding a form similar to the isentropic efficiency as denoted by Eq. (9),

where ** x** represents a location in the compressor passage, and the approximation arises due to the assumption of small temperature changes (i.e., $Tt2,isen\u2248Tt2$, which is reasonable especially for low-speed compressors). The time-averaged non-dimensionalized entropy fields at both the PE and NS conditions are compared in Fig. 9. For both conditions, entropy is generated mainly due to the formation of the main TLV (i.e., mixing between the main flow and the tip leakage flow); when the operating point moves from PE to NS, the blade loading at the tip increases, resulting in a stronger mixing process at the main TLV and hence a lower isentropic efficiency. After the main TLV breakdown, the produced entropy propagates downstream and diffuses to lower spans.

### C. Compressor tip leakage flow topology

To conclude the observations from the flow fields, the flow topology of the investigated compressor rotor is illustrated in Fig. 10. The near-tip flow field can be roughly divided into four zones: formation of main TLV (Zone I); main TLV breakdown (Zone II); formation of tip blockage cell (Zone III); formation of secondary TLV (Zone IV). In Sec. IV, the modal component of each zone will be discussed in detail.

## IV. FLOW FIELD MODE DECOMPOSITION

In this section, the SPOD method is performed to extract and analyze the modal components of compressor tip leakage flow. Two case studies are performed. The first case investigates the near-tip flow field quantities $q=[p,ux,ur,ut,\Delta s]T$ at the mesh points between 80% span and the blade tip, which is relevant to compressor aerodynamics and turbulence modeling. The second case investigates the static pressure *p* on the blade surfaces at 90% span, which is of interest to compressor aeroelasticity. For both cases, the sampling frequency of the snapshots is $fs/BPF=50$, which cuts off the high-frequency components of the subgrid scales. Details of the snapshot data and the SPOD parameters for both cases are listed in Table II.

Database . | SPOD . | |||||||
---|---|---|---|---|---|---|---|---|

Case . | Variable . | N
. _{g} | N
. _{t} | $fs/BPF$ . | N
. _{f} | N
. _{o} | N
. _{b} | Norm . |

Near-tip flow field | $p,ux,ur,ut,\Delta s$ | 481493 | 4250 | 50 | 256 | 128 | 32 | Total disturbance energy |

Blade surface pressure | p | 304 | 4250 | 50 | 256 | 128 | 32 | Control volume |

Database . | SPOD . | |||||||
---|---|---|---|---|---|---|---|---|

Case . | Variable . | N
. _{g} | N
. _{t} | $fs/BPF$ . | N
. _{f} | N
. _{o} | N
. _{b} | Norm . |

Near-tip flow field | $p,ux,ur,ut,\Delta s$ | 481493 | 4250 | 50 | 256 | 128 | 32 | Total disturbance energy |

Blade surface pressure | p | 304 | 4250 | 50 | 256 | 128 | 32 | Control volume |

### A. Near-tip flow field

#### 1. Mode energy spectrum of near-tip flow field

For the multi-variable SPOD analysis performed in this case, the SPOD mode energy represents the total disturbance energy which can be subdivided to pressure potential energy, entropy potential energy, and TKE, as shown in Eq. (7). The fraction of each component of the total disturbance energy is presented in Fig. 11. In this plot, the fractions of energy components are calculated at each snapshot; the time variations of the fractions are presented by probability density functions (PDF), which are obtained by using the kernel density estimation^{44} with Gaussian kernels. The results show that the fraction of TKE = $Eux+Eur+Eut$ is roughly three orders of magnitude stronger than the entropy potential $E\Delta s$ and the pressure potential *E _{p}* at both the PE and NS conditions. Therefore, the total disturbance energy is almost the same as the TKE in this case.

The SPOD modal energy spectra of the near-tip flow field for the PE and NS conditions are shown in Figs. 12(a) and 12(b), respectively. In general, the near-tip flow exhibits a broadband characteristic without any spikes in the spectrum; the majority of the flow energy is contained in the large-scale low-frequency components of the flow, and the energy is transferred successively to the smaller scales in the inertial subrange. A low-rank behavior is observed between 0.4 and 3.0 BPF at the PE condition and between 0.2 and 3.0 BPF at the NS condition, where Mode 1 contains a significantly large portion of the energy and thus plays a dominant role. This low-rank behavior indicates that a physically dominant mechanism is present in the flow field, which will be seen later in the mode shapes. To illustrate the low-rank behaviors further, the premultiplied energy spectra are presented in Figs. 12(c) and 12(d), where the energy contained in each mode equals the area beneath the curve of the mode. For the total energy spectrum, the NS condition contains a larger amount of energy than the PE condition especially at low frequencies below 1.0 BPF. This is because the NS condition has a stronger adverse pressure gradient than the PE condition, which leads to more severe flow separation and hence a higher magnitude of flow field fluctuations. For the Mode 1 spectrum, an energy hump is formed near 1.6 BPF at the PE condition, and two energy humps are formed near 0.4 and 1.0 BPF at the NS condition. The energy contained in these humps of Mode 1 is significantly larger than the rest of the modes, and thus the low-rank behavior.

To quantify the contribution of each mode at each frequency to the total energy, the accumulated mode energy fraction for the *n _{th}* frequency and

*m*mode is calculated as follows:

_{th}The accumulated mode energy fractions of both the PE and NS conditions are plotted as isolines in Fig. 13, where the horizontal axis is the non-dimensionalized frequency and the vertical axis is the number of mode. When using a large number of modes ($m\u226520$), the NS condition requires fewer frequencies than the PE condition to restore the same fraction of total energy. It indicates that the low-frequency components of the NS condition are more energetic than the PE condition, which is consistent with the observation in Fig. 12. When using a large number of frequencies ($fn\u226510$ BPF), the NS conditions requires more modes than the PE condition to restore the same fraction of total energy, indicating the modal component is more complex (i.e., non-low-rank) in the NS condition than the PE condition. To restore 90% of the total energy, it requires at least $fn>4$ BPF and *m *>* *21, which is equivalent to using $mn/NbNf>10.5%$ of the total modes.

#### 2. Mode shape of near-tip flow field

To unveil the physical relevance of the above SPOD modes, the Mode 1 shapes at 1.0, 1.6, and 5.1 BPF for the PE condition and at 0.4, 1.0, and 5.1 BPF for the NS condition are plotted in Fig. 14. In this figure, the wavepackets are represented by the iso-surfaces of mode shapes at 50%, 67%, and 100% of the minimum/maximum values of the color range. The frequencies of 1.0 and 1.6 BPF of the PE condition and 0.4 and 1.0 BPF of the NS condition correspond to the low-rank behavior (i.e., energy humps); the frequency of 5.1 BPF at both the PE and NS conditions represents the high-frequency non-low-rank region. For brevity, high-order modes are not presented since the energy contained in these modes is significantly smaller than the leading mode; only the mode shapes of *u _{x}* and

*p*are analyzed since they are of most interest for compressor aerodynamics and aeroelasticity.

For the PE condition, the mode shapes of *u _{x}* and

*p*at both 1.0 and 1.6 BPF show two coherent wavepackets: the dominant wavepackets is located in Zone III (Fig. 10, same reference for the following zones), which corresponds to the vortex shedding downstream the main TLV breakdown; a second wavepackets can be observed on the trajectory of the secondary TLV in Zone IV, but it is much weaker than the primary wavepackets. At 5.1 BPF, two small-scale wavepackets occur in Zones II and IV, which corresponds to the main TLV breakdown and the secondary TLV shedding, respectively. These results indicate that the main TLV shedding in the tip blockage cell zone is most energetic in the flow domain, which creates the energy hump of Mode 1 near 1.6 BPF. For the NS condition, a similar main TLV shedding mode in Zone III can be observed at 0.4 and 1.0 BPF; an oscillating mode of the main TLV around its time-averaged trajectory is also observed in Zone I at 0.4 BPF. These modes lead to the two energy humps of Mode 1 near 0.4 and 1.0 BPF. At 5.1 BPF, small-scale wavepackets occur in Zones II and IV, which is similar to the PE condition.

Given the physical phenomena behind the energy hump is the main TLV shedding, the next question is: what determines the frequency range of the energy hump? The shedded vortices create wave-like pulsations of flow quantities in blade passages: the wave speed is the propagation speed of the vortices, which equals to the local relative velocity *W*; the wavelength is approximated by the distance *l* between two adjacent vortex cores, which can be measured by a pair of wavepacket shown in Fig. 10; thus, the wave frequency $f=W/l$. This relation can be non-dimensionalized by using the BPF, the blade tip speed *U _{t}*, and the blade tip pitch

*s*(note BPF=$Ut/s$), as denoted in Eq. (11),

It shows that the frequency range of the energy hump depends on two terms: $W/Ut$ which depends purely on the operating condition of the compressor, and *l*/*s* which depends on both the operating condition and the design of the compressor. By substituting the values in Figs. 8 and 14 into Eq. (11), this relation is checked to be correct in principle.

The $W/Ut$ term decreases when throttling the compressor from PE to NS (from analysis of the inlet velocity triangle, or more intuitively from the observation in Fig. 8), leading to the shift of the energy hump toward the low-frequency direction shown in Fig. 12. When operating at the same operating condition, the wavelength decreases when increasing the frequency, which is captured in Fig. 14.

The minimum *l*/*s* value is determined by the Kolmogorov scale of the vortices, which scales with $Re\u22123/4$. For large-scale high-speed industrial machines with high $Re$ numbers, more high-frequency components are expected in the energy hump. The maximum *l*/*s* value depends on the trajectory of the main TLV, which changes with the operating condition and the design of the blade. Based on the 2D blade geometry of Fig. 15, the maximum possible *l*/*s* value is $\sigma 2+2\sigma \u2009sin\u2009\xi +1$, where $\sigma =c/s$ is the solidity and *ξ* is the stagger angle. For compressor designs with high values of solidity and stagger angle, more low-frequency component is expected in the energy hump.

### B. Blade surface pressure

#### 1. Mode Energy Spectrum of Blade Surface Pressure

The premultiplied SPOD modal energy spectra of the blade surface pressure are shown in Figs. 16(a) and 16(b). Similar to the spectrum of the near-tip flow field, a low-rank behavior is observed between 0.4 and 3.0 BPF at the PE condition and between 0.2 and 3.0 BPF at the NS condition, and the peak frequencies are 1.6 BPF for the PE condition and 0.4 BPF and 1.0 BPF for the NS condition. The difference is that the low-rank behavior is more evident in the spectrum of blade surface pressure compared to the near-tip flow field.

The more evident low-rank behavior is further confirmed in Fig. 17, where the accumulated mode energy fractions of both the PE and NS conditions are plotted. To restore 90% of the total energy, it requires at least $fn>3$ BPF and *m *>* *4, which is equivalent to using only $mn/NbNf>1.5%$ of the total modes. Comparing the NS condition with the PE condition, the low-frequency components are more energetic and the high-frequency components are more non-low-rank, which is consistent with the analysis of the near-tip flow field.

#### 2. Mode shape of blade surface pressure

The Mode 1 shapes in the frequency range of the low-rank energy hump are presented in Fig. 18. To help correlate the mode shape with the mode energy, the premultiplied energy spectrum is plotted on the side for reference. Since the mode shape on the SS is comparatively insignificant, only the PS is shown. The results show that coherent wavepackets exist on the PS at both conditions, with the affected region more upstream in the NS than the PE condition. Combined with the mode shape shown in Fig. 14, these pressure waves correspond to the vortex shedding of the main TLV in Zone III. When throttling the compressor toward stall, the main TLV trajectory moves toward the leading edge [Fig. 7(b)] and hence the affected region by the pressure waves moves upstream.

#### 3. Reconstruction of blade forces

To leverage the low-rank feature of the blade surface pressure data, this section explores the feasibility of reconstructing the pressure data using low-order SPOD modes at low frequencies. The frequency-domain reconstruction method from Nekkanti and Schmidt^{45} is performed. This method is equivalent to the inverse of SPOD, which can reconstruct the solution by using selected modes and frequencies. In the following, the reconstruction results using the leading *m* modes (*m *=* *1, 2, 5) and *n* frequencies (*n *=* *102, *f _{n}* = 10 BPF) will be presented, where the choice of

*m*and

*n*are based on the observation from Fig. 17. These results will be denoted by SPOD(

*m*×

*n*).

The time-averaged static pressure rise coefficient *C _{p}* on the 90% span of the blade surfaces is presented in Fig. 19(a). The overall loading is similar for both the PE and NS conditions, but the NS condition has a higher loading between 0 and 20% chord and hence a stronger main TLV. The root mean square (RMS)

*C*profiles are presented in Figs. 19(b) and 19(c). For both operating conditions, the unsteadiness is stronger on the PS than the SS, which is consistent with the previous observations. When throttling the compressor from the PE to NS condition, the unsteadiness intensifies, and the affected region moves toward upstream. Comparing DDES with SPOD reconstruction results, it is seen that the $Cp,RMS$ profiles of both conditions can be reconstructed with good accuracy using the leading

_{p}*m*=

*5 modes and*

*n*=

*102 frequencies (i.e., $fn<10$ BPF).*

The primary modes of vibration for a 2D airfoil can be divided into a plunging, edgewise and pitching motion (Fig. 20), and the corresponding dimensionless aerodynamic forces are *C _{L}*,

*C*, and

_{D}*C*. The moment is calculated about the aerodynamic center (a.c.). The premultiplied PSD of

_{M}*C*,

_{L}*C*, and

_{D}*C*is presented in Fig. 21, where the frequency in engine order ($EO=Z\xb7f/BPF$) is also shown on the top axis for reference. At the PE condition, an energy hump near 1.0 BPF is observed in the spectra, but its intensity of the premultiplied PSD is comparable to the noise near 6.0 BPF and thus small. At the NS condition, the energy hump has a frequency range between 0.3 and 1.0 BPF, within which the intensity of the premultiplied PSD is evidently larger than the high-frequency noise. For both operating conditions, using the leading

_{M}*m*=

*5 modes and*

*n*=

*102 frequencies can reconstruct the premultiplied PSD of*

*C*,

_{L}*C*, and

_{D}*C*accurately. The success of SPOD-based low-rank reconstruction of blade forces is the first step toward a reduced-order model.

_{M}## V. FURTHER DISCUSSIONS

### A. RANS turbulence modeling of tip leakage flow

From the SPOD mode energy spectra and mode shapes of the near-tip flow field, it is suggested that the turbulent fluctuations in Zones I and III are conceptually easy to model because of the low-rank behavior at low frequencies. These fluctuations form a hump in the energy spectrum, which corresponds to turbulence production. Previous RANS simulations have inferred that the SA model underpredicts the eddy viscosity in the tip blockage cell and hence leads to premature stall.^{28,46} This work confirms that modifying the production term of the SA model in the tip blockage cell is conceptually correct.

On the contrary, modeling the turbulent fluctuations in Zones II and IV is conceptually difficult due to the non-low-rank modal component at high frequencies. Zone II involves the laminar-to-turbulent transition process, which is related to the transition term in RANS turbulence models. The amplification mechanism during the turbulence transition process is not discussed in this work, which needs to be investigated for future research. Zone IV shows the interaction between the secondary TLV and the SS boundary layer flows, the latter of which is related to the destruction term in RANS turbulence models. Such an observation indicates the difficulty of modeling the destruction term in internal flows.

### B. Nonsynchronous vibration and tip leakage flow

From previous measurements, nonsynchronous vibration usually occurs on the first or second torsional mode whose frequency is between 5 and 20 EO.^{3,47} The numerical results in this work (i.e., single-passage DDES simulations) suggest that strong unsteady blade forces between 5 and 20 EO can occur at the NS condition, which are mainly induced by the vortex shedding downstream the main TLV breakdown. When considering a full-annulus flow domain in the simulation, the vortex shedding in one blade passage can interact with that in another, which potentially forms circumferentially traveling cells that intensify the low-frequency unsteady blade forces further. When considering the blade structural dynamics in the simulation, the low frequencies of aerodynamic blade forces can lock-in with the natural frequency of the blade, eventually triggers nonsynchronous vibration. To consolidate the above analysis, further simulations with a full-annulus flow domain and blade motion are needed for future research.

## VI. CONCLUSIONS

In this paper, the unsteady modal component of compressor tip leakage flow has been analyzed for the first time by using the SPOD method. The investigated case is a low-speed compressor rotor, and its flow field data were obtained by DDES simulations that were validated against the SPIV measurements. Several conclusions can be drawn as follows.

The tip leakage flow in the investigated low-speed compressor rotor is governed by the double-leakage flow structure. The main TLV forms immediately after detaching from the leading edge (Zone I), breaks down into smaller eddies near the mid-pitch (Zone II), and then impinges on the PS of the adjacent blade and forms a tip blockage cell (Zone III); part of the main TLV flow passes through the tip gap of the adjacent blade, leading to the formation of a secondary TLV (Zone IV). When changing the operating condition from PE to NS, the trajectory of both TLVs moves toward upstream, and the areas of these zones change accordingly.

The SPOD modal analysis on the near-tip flow field revealed that a major part of the total disturbance energy originates from the fluctuations in Zones I and III, which are low-frequency (i.e., $0.2<f/BPF<3.0$ in this case), low-rank, and thus conceptually easy to model by RANS turbulence models. These fluctuations correspond to the oscillating mode of the main TLV in Zone I and the vortex shedding mode of the main TLV in Zone III. In particular, the frequency range of the shedding mode can be estimated based on the blade geometry and the operating condition. On the contrary, the fluctuations in Zones II and IV are broadband, non-low-rank, and thus conceptually difficult to model by RANS turbulence models. These fluctuations correspond to the laminar-to-turbulent transition in Zone II and the interaction between the secondary TLV and the SS boundary layer in Zone IV. When operating from the PE to NS condition, the total disturbance energy increases, and the energy hump corresponding to the main TLV modes moves toward the low-frequency; however, the low-rank behavior of the main TLV modes becomes less evident.

The SPOD modal analysis on the blade surface pressure showed that the unsteady blade forces due to the impingement of the main TLV on the PS in Zone III are more significant than that created by the detachment of the secondary TLV from the SS in Zone IV. When operating at the NS condition, the unsteady components of the blade forces become strong at low frequencies (i.e., $5<EO<20$ in this case), which can potentially lead to nonsynchronous vibration. SPOD-based low-rank reconstruction can restore most of the blade force signals by using only the leading 5 modes at frequencies below 10 BPF, which highlights the opportunity of a SPOD-based reduced-order model for predicting the unsteady blade forces.

## ACKNOWLEDGMENTS

We thank Professor Nick Cumpsty of Imperial College for helpful comments and feedback, Professor Baojie Liu of Beihang University for sharing the Stage B rotor data, Professors Oliver Schmidt and Akhil Nekkanti of UC San Diego for the input of SPOD reconstruction, and NUMECA for providing the student license of its mesh generator AutoGrid. Xiao He greatly acknowledges the Imperial College President Ph.D. Scholarship for funding this research.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## NOMENCLATURE

- BPF
blade passing frequency

- CSD
cross-spectral density

- DDES
delayed-detached eddy simulation

- EO
engine order

- LES
large eddy simulation

- NS
near-stall

- PE
peak efficiency

- PS
pressure surface

- PSD
power spectrum density

- RANS
Reynolds-averaged Navier–Stokes

- SPIV
stereoscopic particle image velocimetry

- SPOD
spectral proper orthogonal decomposition

- SS
suction surface

- TLV
tip leakage vortex

*c*blade chord length (

*m*)*C*_{p}pressure rise coefficient

*f*_{s}sampling frequency (Hz)

- M
Mach number

*N*_{b}number of blocks

*N*_{f}number of frequencies

*N*_{g}number of grid points

*N*_{o}number of overlapped snapshots

*N*_{t}number of snapshots

- Re
Reynolds number

*U*_{m}blade speed at mid-span (

*m*/*s*)*U*_{t}blade speed at tip-span (

*m*/*s*)- Δ
grid spacing (

*m*)*λ*SPOD mode energy

*τ*tip gap size (

*m*)- $\varphi $
SPOD mode shape

- Ω
_{rot} blade angular velocity ($s\u22121$)