Nomenclature
- θ
Circumferential angle
- k
Number of decomposed modes
- p
Static pressure
- f
Frequency
- rot
Rotor
- t
Time
-
${\omega _k}$
-
Centre frequency
- α
-
Penalty factor
- λ
-
Lagrangian multiplication operator
-
${\eta _{\rm{e}}}$
-
Energy ratio
-
${\eta _{\rm{s}}}$
-
Similarity factor
- BPF
-
Blade pass frequency
- EMD
-
Empirical mode decomposition
- HHT
-
Hilbert–Huang transform
- L
-
Extended Lagrange expression
- IMF
-
Intrinsic mode function
- RI
-
Rotating instability
- RF
-
Rotational frequency
- VMD
-
Variational mode decomposition
-
$H\left( {\omega ,t} \right)$
-
Hilbert amplitude spectrum
-
$h(\omega)$
-
Hilbert marginal spectrum
-
${u_k}(t)$
-
Each decomposed IMF
-
${z_k}(t)$
-
Analytic signal
-
${\phi _k}(t)$
-
Instantaneous phase
1. Introduction
Nowadays, turbomachines play a significant role in industrial installations, such as power plants, compressed air energy storage, automotive vehicles and small aircraft. Generally, turbomachines are grouped into two categories. One mainly includes the steam and turbines, which are used to produce power. The other group comprises machines that are used to increase the flow fluid pressure, such as compressors, pumps and fans. The compressor may suffer from different flow instabilities under a low mass flow rate. These type flow phenomena, such as stall and surge, have been investigated through experiments and numerical simulations by several researchers.
Surge is a global, large-amplitude flow oscillation in the compressor system, which causes a loud noise and mechanical vibration. Most notably, surge is a low-frequency and axisymmetric flow. Unlike surge, the rotating stall is a local and circumferential flow disturbance. This phenomenon exists widely in both axial and centrifugal compressors. When a compressor enters rotating stall, the rotating cell propagates from blade to blade at some frequency, loading and unloading each blade. Many researchers have recognised two types of stall inception: spike and modal wave inception [Reference McDougall, Cumpsty and Hynes1, Reference Day2]. The modal wave is a small amplitude and low frequency disturbance. Greitzer et al. [Reference Greitzer3] used the Moore-Greitzer model to experimentally investigate modal oscillation behaviour. Based on fast wavelet analysis, Liu et al. [Reference Liu, Li, Du, Li and Zhang4] achieved low-frequency reconstruction and detected modal wave inception at 110 revolutions before stall in an axial compressor. Arshad et al. [Reference Arshad, Li, Li and Pan5] reported that modal waves were related to the separation flow near the hub region, whereas they were not related to the origin of the stall disturbance location. Vo et al. [Reference Vo, Tan and Greitzer6] summarised two criteria for judging the occurrence of spike-type stall: the spillage and the backflow of tip leakage flow occurs at the blade leading edge and trailing edge, respectively.
For centrifugal compressors, researchers have done many studies on the causes of and location of stall inception. Spakovszky et al. [Reference Spakovszky and Roduner7] observed that the spike-type stall occurred the in the vaneless space between the impeller and vane diffuser. Everitt et al. [Reference Everitt and Spakovszky8] revealed that vortical structures increased diffuser inlet blockage and triggered short-wavelength stall in an isolated diffuser. Bousquet et al. [Reference Bousquet, Binder, Dufour, Carbonneau, Trebinjac and Roumeas9, Reference Bousquet, Binder, Dufour, Carbonneau, Trebinjac and Roumeas10] found a modal wave stall inception in a centrifugal compressor with a vane diffuser (without a volute). They suggested that the stall wave can propagate upstream and interact with the flow structure at the impeller inlet. Owing to the complex inlet bent torsional pipes and volute structure, the internal flow field become more complicated, thereby affecting the initial occurrence position of stall. Yang et al. [Reference Yang, Wang, Lao, Wang, Yang, Chen and Zhang11–Reference Yang, Zhang, Yang, Lao and Yang13] conducted a lot of experimental and numerical studies on the stall in a centrifugal compressor with a volute. They suggested that volute tongue induced an adverse pressure wave, and the wave propagated backwards into a certain circumferential position at the impeller inlet. The stall first occurred at this position. Zhang et al. [Reference Zhang, Yang, Yang, Zhang, Wang and Chen14] summarised the influence of inlet and outlet asymmetric structure on flow field distortion, and they suggested that the compressor stability was reduced when the high static pressure region induced by the inlet bent torsional pipe was located at the same circumferential positions of the high static pressure region induced by the volute.
Apart from surge and rotating stall, rotating instability (RI) is another focus of the flow instability. Kameier et al. [Reference Kameier and Neise15, Reference Kameier and Neise16] suggested that RI can be decomposed by the spatial Fourier transform. RI only appears when there were a reverse flow in the blade gap. Mailach et al. [Reference Mailach, Lehmann and Vogeler17, Reference Mailach, Sauer and Vogeler18] found that blade tip vortices caused RI. Vo et al. [Reference Vo19] suggested that the RI was due to flow oscillations caused by the tip clearance backflow. Based on acoustic measurements, Sun et al. [Reference Sun, Zou and Zheng20] found that there were four instability phenomena in a centrifugal compressor: rotating stall, RI, mild and deep surge. Li et al. [Reference Li, Wu and Ouyang21, Reference Li, Wu and Ouyang22] suggested that surface feature was related to the RI behaviours. For a transonic axial compressor, Ye et al. [Reference Ye, Zhao, Zhou, Xi and Xu23] investigated the origin and propagation of RI and discussed the influence of circumferential casing grooves on this instability phenomenon. Zou et al. [Reference Zou, He, Zhang, Niu and Zheng24] experimentally investigated the dynamic pressure for two types of centrifugal compressors. At low and middle speeds, the RI was not almost affected by diffuser vanes. In addition, He [Reference He and Zheng25] also adopted a casing treatment to eliminate the rotating instability. For a transonic rotor, Wang et al. [Reference Wang, Lu, Liu and Hu26] revealed the RI phenomenon with a circumferential propagation “multi-passage structure.” The RI appeared at the rotor tip region when the compressor operates in the stable range. Other studies [Reference Holzinger, Wartzek, Schiffer, Leichtfuss and Nestle27–Reference Wang, Yang, Hu and Zhang30] also have shown the relation between tip leakage flow and RI.
Numerous experiments have indicated that flow instabilities can significantly decrease the performance and operating range of compressors. Therefore, it is necessary to detect the pressure signal from compressor experiments under different operating conditions effectively and quickly. Researchers have used different methods for instability detection, such as time series analysis, time and spatial Fourier transform, and time-frequency analysis methods including short time Fourier transform and wavelet transform. Time series analysis can only study the time domain characteristics of signals, and the Fourier transform can only study the frequency characteristics. The temporal and spatial characteristics of instability signals has been conveniently obtained using a time-frequency method. The one-dimensional continuous wavelet transform has been utilised in studies on compressor stall [Reference Lin, Chen and Li31, Reference Lin, Li and Chen32]. Cameron et al. [Reference Cameron and Morris33] employed different techniques to investigate two representative stall inception events, and they also proposed a windowed correlation technique for the signal detection. For a low-speed contra rotating axial flow fan, Toge et al. [Reference Toge and Pradeep34, Reference Dejene Toge and Pradeep35] studied the stall inception mechanism from unsteady flow measurements using the wavelet transform. Recently, Liu et al. [Reference Liu, Li, Du, Zhang and Nie36] employed various techniques to investigate stall warnings, such as root mean square, auto-correlation, cross-correlation, and fast wavelet analysis. Champ et al. [Reference Niccolini Marmont Du Haut Champ, Ferrari, Silvestri and Massardo37] suggested that wavelet analysis and short-time Fourier transform can effectively detect incipient surge conditions. Li et al. [Reference Li, Zheng, Chen, Duan, Jiang and Benini38] also adopted the wavelet transform technique to divide the rotating stall process into RI and stall inception. Moreover, Li et al. [Reference Li, Du, Liu, Zhang and Nie39] and Manas et al. [Reference Manas and Pradeep40, Reference Manas and Pradeep41] also employed the discrete spatial Fourier series or wavelet transform to study the aerodynamic stability of compressor. Wang et al. [Reference Wang, Chen, Dong and Hua42] observed the RI using the wavelet features.
Recently, novel nonlinear analysis methods have developed rapidly, and can give a better choice to decompose nonlinear signals. Huang et al. [Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu43] proposed creatively a revolutionary time-frequency method called the Hilbert-Huang transform (HHT) combining the empirical mode decomposition (EMD) with Hilbert transform. By this method, a non-stationary and nonlinear signal was be decomposed into different intrinsic mode functions (IMFs), which is an appropriate method for a pressure signal measurement of a compressor, various types of instability signals are typically nonlinear and non-stationary. Li et al. [Reference Li, Lei and Fu44] employed the HHT method based on EMD and obtained IMFs in different frequency bands, which can be used to detect aerodynamic instability. Wu et al. [Reference Wu, Liu, Liu and Zhao45] suggested that the EMD technique can discover the compressor surge earlier than the Morlet continuous wavelet transform. Zhao et al. [Reference Zhao, Xi, Wang and Zhao46] employed this method to show the propagation from the normal condition to spike-wave stall inception. Liu et al. [Reference Liu, Ma, He and Gao47] proposed a method, namely EEMD, which combines an improved multi-scale fuzzy entropy with an improved EMD, to investigate unstable operating states of a centrifugal compressor. Guan et al. [Reference Guan, Sun, Xu, Bishop, Sun, Ni, Du and Zhao48] applied the EMD to observe the compressor stall mode. Yue et al. [Reference Yue, Wang, Wei, Zhang and Wang49] proposed a novel technique which combines an EMD with a local mean decomposition technique. They applied it to investigate the stall data which is nonlinear and non-stationary. However, EMD has shortcomings, such as the mode mixing and the insufficiency of the mathematical theory. Recently, based on variational theory, Dragomiretskiy [Reference Dragomiretskiy and Zosso50] proposed a new, non-recursive, adaptive technique which is called variational mode decomposition (VMD), to overcome these drawbacks. Compared with EMD, VMD has a strict mathematical theory. Specifically, it is an appropriate approach for investigating nonlinear and non-stationary signals. Therefore, for compressors operating in an unstable operating range, the VMD method can be considered an effective means to investigate the nonlinear characteristics. In addition, some novel methods were introduced to analyse the flow field in compressors such as dynamic mode decomposition [Reference Hu, Yang, Yi, Zheng, Zou and Zhou51] and global stability approach [Reference Yang, Fu, Hu and Shi52].
Most studies have applied traditional methods to study flow instability. However, it is necessary to apply more advanced signal techniques to investigate the instability owing to the presence of the nonlinear and non-stationary characteristics of the signal. The current study aimed to investigate the evolution and circumferential nonuniformity characteristics of flow instabilities under various operating points at 80krpm for a vehicle centrifugal compressor. The study was conducted in the following order:
-
(1) For a centrifugal compressor with a volute, a dynamic pressure experiment was carried out to obtain the dynamic pressure signals from the near choke to surge.
-
(2) Fourier transform and VMD methods were used to detect the stable and unstable conditions of the compressor, respectively. The flow instability process was investigated in detail using Hilbert marginal spectrum.
-
(3) The circumferential nonuniformity phenomenon of flow instability is discussed using unsteady simulations.
2 Brief description of VMD
In this study, VMD technique was adopted to extract the corresponding sub-signals from a pressure signal. VMD could decompose a multicomponent signal into a series of sub-signals. Every sub-signal (IMF) is an amplitude-modulated and frequency-modulated, which can be constructed as follows:
\begin{align}{u_k}(t) = {A_k}(t)\cos ( {{\phi _k}(t)}) \end{align}
where
${\phi _k}(t)$
is a non-decreasing function and
${A_k}(t)$
is the instantaneous amplitude. For each mode, it is assumed to be compacted mostly around the corresponding centre frequency
${\omega _k}$
.
Based on a Hilbert transform, the analytic signal of each IMF is acquired, and its unilateral frequency spectrum is obtained as follows:
\begin{align}\left[ {\delta (t) + \frac{j}{{\pi t}}} \right] * {u_k}(t) \end{align}
where
$\delta (t)$
is the Dirichlet function.
For each mode, an exponential term
${e^{ - j{\omega _k}t}}$
is used to shift the frequency spectrum of the analytic signal to the base band:
\begin{align}\left\{ {\left[ {\delta (t) + \frac{j}{{\pi t}}} \right] * {u_k}(t)} \right\}{e^{ - j{\omega _k}t}} \end{align}
A constrained variational problem is constructed as follows:
\begin{align}\left\{ \begin{array}{l}\mathop {\min }\limits_{\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\}} \left\{ {\sum\limits_{i = 1}^K {\bigg\| {{\partial _t}\left\{ {\left[ {\delta (t) + \dfrac{j}{{\pi t}}} \right] * {u_k}(t)} \right\}{e^{ - j{\omega _k}t}}} \bigg\|} _2^2} \right\}\\ \\[-7pt] {\rm{s}}{\rm{.t}}{\rm{. }}\sum\limits_{i = 1}^K {{u_k}} (t) = x(t)\end{array} \right. \end{align}
Where
$\left\{ {{u_k}} \right\}: = \left\{ {{u_1}, \ldots ,{u_K}} \right\}$
are the modes, and
$\left\{ {{\omega _k}} \right\}: = \left\{ {{\omega _1}, \ldots ,{\omega _K}} \right\}$
are the corresponding centre frequencies, respectively.
The solution method is to convert this constrained problem into the unconstrained variation problem; therefore, VMD introduces the penalty factor α and Lagrangian multiplication operator λ. The extended Lagrange expression is the following formula:
\begin{align}\begin{array}{l}L\left( {\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\},\lambda } \right) = \alpha \sum\limits_{i = 1}^K {\bigg\| {{\partial _t}\left\{ {\left[ {\delta (t) + \dfrac{j}{{\pi t}}} \right] * {u_k}(t)} \right\}{e^{ - j{\omega _k}t}}} \bigg\|_2^2} + \left\| {f(t) - \sum\limits_{k = 1}^K {{u_k}(t)} } \right\|_2^2 + \\[3pt] \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\langle {\lambda (t),f(t) - \sum\limits_{k = 1}^K {{u_k}(t)} } \right\rangle \end{array} \end{align}
The alternate direction method of multipliers (ADMM) is adopted to iteratively find the solution to the saddle point of the L. More details about the ADMM algorithm are provided in Ref. (Reference Dragomiretskiy and Zosso50). VMD can effectively deal with a complex, nonlinear and multicomponent signal. Each sub-signal or IMF has a very limited bandwidth and is around its corresponding characteristic centre frequency
${\omega _k}$
. In addition, the effect of the VMD method is strongly associated with the two key parameters, penalty factor α and mode number K. Therefore, one of main purpose in next section is to investigate how to obtain the optimal parameters for effective decomposition.
3 Selection of the key parameters and method validation
The mode number K represents the number of sub-signals into which the origin signal has been decomposed. The penalty factor α is related to the mode bandwidth and data fidelity, also called the balancing parameter. These two parameters affect the decomposition effect of the VMD method. Depending on whether the mode number is small or large, the results can be either under or over decomposed. These two important variables must be defined in advance when the VMD is employed to process a complex compound signal.
Before using the VMD method to decompose dynamic experiment signals, both parameters K and α must be defined. The effect of these two parameters was investigated by designing a multicomponent compressor signal x(t) at 80krpm. The input signal included typical compressor frequency components, namely the blade passing frequency (BPF), rotation frequency (RF), RI frequency, spike-wave stall frequency, modal-wave stall frequency, rotor frequency and surge frequency. The sampling frequency
$f_{\textrm{s}}$
for the simulation signal was 200kHz. The signal x(t) and each sub component are expressed in Equations (6)–(12), and illustrated in Figs. 1 and 2.
\begin{align}{x_1} = 100\cos\!(2\pi {f_{{\rm{BPF}}}}t) \end{align}
\begin{align}{x_2} = 50\sin\!(2\pi {f_{{\rm{RI}}}}t) \end{align}
\begin{align}{x_3} = 30\cos\!(2\pi {f_{{\rm{rot}}}}t) \end{align}
\begin{align}{x_4} = 150\sin\!(2\pi {f_{{\rm{spike}}}}t) \end{align}
\begin{align}{x_5} = 150\sin\!(2\pi {f_{{\rm{modal}}}}t) \end{align}
\begin{align}{x_6} = 250{\mathop{\rm sins}\nolimits} (2\pi {f_{{\rm{surge}}}}t) \end{align}
\begin{align}x(t) = {x_1} + {x_2} + {x_3} + {x_4} + {x_5} + {x_6} \nonumber\\[-10pt] \end{align}
where
$f_{\textrm{BPF}}$
= 9,333.33Hz,
$f_{\textrm{RI}}$
= 4,666.67Hz,
$f_{\textrm{rot}}$
= 1,333.33Hz,
$f_{\textrm{spike}}$
= 800Hz,
$f_{\textrm{modal}}$
= 266.67Hz, and
$f_{\textrm{surge}}$
= 66.67Hz.
Simulation signal x(t).
Frequency spectrum of simulation signal x(t).
Different parameter combinations were selected to study the effect of these two parameters on the decomposition effect. Owing to large quantities of sampled data, the time-domain diagram of each component is unclear, therefore, only the frequency spectrum of each component was displayed after VMD. Figure 3 shows the frequency spectra of the decomposition results at k = 5, which is less than the number of original signal components. When penalty factor α is very small (α = 0.01 f s), it can be observed that there is a frequency component of 266.67Hz in IMF4 and IMF5. Multiple components of the signal appear in an IMF simultaneously. When the penalty factor α is extremely large (α = 4 f s), it can be found that the component with the frequency of 4,666.67Hz is discarded.
Decomposition results of VMD at k = 5.
Figure 4 shows the frequency spectra of the decomposition results at k = 6, which is the number of original signal components. The penalty factor α is set to 0.01 f s, 0.3 f s and 10 f s. At α = 0.01 f s, it can be observed that the component with the frequency component of 266.67Hz is decomposed into IMF5 and IMF6. At α = 0.3 f s, it is obvious that the frequency component of each IMF is highly consistent with that of the original signal components. At α = 10 f s, the corresponding components of 1,333.33Hz is discarded. The corresponding components of 66.67Hz is decomposed into IMF5 and IMF6.
Decomposition results of VMD at k = 6.
Figure 5 shows the frequency spectra of the decomposition results at k = 7, which is larger than the number of original signal components. When penalty factor α is very small (α = 0.01 f s), the component with a frequency of 4,666.67Hz is decomposed into IMF2 and IMF3. When penalty factor α is extremely large (α = 10 f s), the component with the frequency of 1,333.33Hz is discarded. Meanwhile, the corresponding components of 800Hz is decomposed into IMF3 and IMF4. The corresponding components of 66.67Hz is decomposed into IMF6 and IMF7.
Decomposition results of VMD at k = 7.
According to the above analysis, an unreasonable mode number K and penalty factor α will cause a mode mixing problem. Mode mixing indicates that there are several different frequency components in the same IMF, or the same frequency component is shared by several different IMFs and the corresponding centre frequencies coincide, namely mode duplication. In addition to mode mixing, a larger α will cause important information to be missed in the signal decomposition. Therefore, it is necessary to choose the appropriate K and α.
At present, there are many indicators (entropy, energy loss coefficient, and correlation coefficient) in signal analysis; however, using a single indicator to measure signal characteristics is not optimal. Hence, a comprehensive evaluation indicator based on the energy ratio and frequency similarity factor was used to solve this problem, thereby obtaining optimal K and α.
In this study, the energy ratio
${\eta _{\rm{e}}}$
was considered as the weight of the similarity factor. The indicator is defined as a ratio of the total energies of all the IMFs to the original signal energy, which can be calculated using the following equation:
\begin{align}{\eta _{\rm{e}}}{\rm{ = }}\frac{{\sum\limits_{k = 1}^M {\sum\limits_{t = 1}^N {{{\left[ {{u_k}(t)} \right]}^2}} } }}{{\sum\limits_{t = 1}^N {{{\left[ {x(t)} \right]}^2}} }} \end{align}
where x(t), M, and N represents the original signal, the mode number, and the total length of the original signal, respectively. The closer the total energy of all IMF components is to the total energy of x(t), the better the decomposition effect is.
For over-decomposition, some spurious components with a close centre frequency will be generated; however, these modes should belong to the same component. Therefore, a similarity factor used to assess the similarity of adjacent centre frequencies is defined as follows:
\begin{align}{\eta _{\rm{s}}} = \frac{{{\omega _{k + 1}}}}{{{\omega _k}}} \end{align}
When the ratio is larger than 0.95, it can be considered that the centre frequency overlaps, which indicates that modal mixing occurs.
Eventually, the operation steps can be summarised as follows:
-
(1) The modal number K is estimated through the signal Fourier spectrogram and the K range is set.
-
(2) VMD is performed for the signal.
-
(3) The energy ratio is calculated and the minimum ratio for the decomposition result for each K is selected.
-
(4) The similarity factor for the decomposition result of the minimum ratio is calculated to obtain a reasonable decomposition.
4 Test facility and numerical method
4.1 Experimental setup
The experiments in this research were performed in the centrifugal compressor with a volute for a vehicle diesel engine. Figure 6 shows the test rig. Figures 7 and 8 show the turbocharger and the centrifugal compressor model, respectively. The compressed air from the air supply drove the turbine, thereby driving the centrifugal compressor. The compressor stage consisted of an inlet pipe, centrifugal impeller, vaneless diffuser, and asymmetric volute. The pressure ratio was 1.92 and the design mass flow rate was 0.26kg/s at the design rotating speed with 90krpm. Table 1 lists more key parameters of the investigated compressor.
Schematic diagram for test rig.
Key parameters of investigated centrifugal compressor
Turbocharger.
Centrifugal compressor model.
4.2 Data acquisition analysis
The experimental measurements included two aspects: the compressor performance and dynamic casing static pressure in the inlet of the impeller. The rotational speed and mass flow rate were collected by an electromagnetic sensor and a lemniscate flowmeter, respectively. At both the inlet and outlet, resistance sensors were used to record the temperature. The Dwyer stainless steel Pitot tube was installed at the compressor inlet and outlet pipes and equipped with the PTX5072-TC-A3-CB-H0-PA (inlet pipe) sensor and the PTX5072-TC-A3-CA-H0-PA (outlet pipe) sensor to measure the total pressure. The Pitot tube has a diameter of 3.18mm, which is 6.36% of the inlet pipe and 5.3% of the outlet pipe. In this experiment, the compressor flow instability was investigated based on the dynamic casing pressure along a circumference at the impeller inlet. At the impeller inlet, six KULITE XT-140M sensors were used to obtain the dynamic pressure, which were placed evenly and circumferentially around the casing wall. More detailed parameters of KULITE XT-140M are shown in Table 2. Figure 9 illustrates the schematic diagram for the sensor locations, marked with P1–P6.
Main parameters of pressure sensor
Schematic diagram for sensor location.
Recorded data of the dynamic experiment.
Compressor speed sweeps (from 70 to 100krpm) were performed from choke to near surge. During this period, the compressor performance characteristics and dynamic static pressure were simultaneously recorded. The compressor could operate at various mass flow rates by adjusting the valve opening. Each data recording started after the compressor operated stably for a period of time. The pressure data acquisition system records lasted 10s at a sample rate of 200kHz to capture the fully developed instability. The occurrence of low-frequency noise indicated that the compressor has entered a deep surge state. Immediately, the record was stopped and the compressor recovered from the surge by increasing the mass flow quickly to avoid damaging the test rig. The recorded data are shown in Fig. 10.
4.3 Numerical procedure
The flow instability in this centrifugal compressor was accurately captured by conducting a full-annulus 3D unsteady numerical simulation with the EURANUS solver of the commercial software NUMECA. For the unsteady simulation, the 3D Reynolds-averaged Navier–Stokes (RANS) equations were discretised with the dual time-stepping technique temporally and with the cell-centered finite volume technique spatially. The Spalart–Allmaras (SA) one-equation model was employed for turbulence closure. This turbulence model has been validated previously for the investigation on flow instability in centrifugal compressors [Reference Dodds and Vahdati53, Reference Pullan, Young, Day, Greitzer and Spakovszky54]. At the inlet, the axial flow, total pressure (101,325Pa), and total temperature (298K) were imposed while the average static pressure was applied at the outlet. The flow field results under different operating conditions were obtained by gradually increasing the outlet static pressure. In addition, during a rotating period, the physical time step was set to 280, and 60 inner iterations per time step were adopted for each physical time step.
The computational domains consist of the inlet pipe, centrifugal impeller, vaneless diffuser, volute, and outlet pipe. Significantly, the upstream inlet pipe domain extended to five times the length of the impeller inlet diameter. The downstream outlet pipe domain also extended to five times the length of the volute outlet pipe diameter. Using the module IGG/AutoGrid5, the structured mesh was generated. The rotor/stator interface connection was used between the inlet pipe and the centrifugal impeller, and between the centrifugal impeller and vaneless diffuser. Prior to performing numerical simulations, the grid independence test was verified to guarantee calculation accuracy and save calculation time. Otherwise, this study only focuses on the flow instability at 80krpm; therefore, simulations were carried out only at this speed. To ensure that the simulation results were not related to the grid distribution, three different grid configurations, coarse mesh (8,143,265), medium mesh (8,510,513) and refined mesh (9,417,769) were conducted by adjusting the node distribution in the span and streamwise directions. Considering the pressure ratio and efficiency from steady simulations, it was proven that the medium mesh was already sufficient for obtaining reliable simulation results. Otherwise, the simulation results from the medium mesh were consistent with the experiments, as illustrated in Fig. 12. Eventually, a medium mesh was used in this study. The y+ value was less than five throughout the computational domain, which meets the requirements of the turbulence model. The CFD model could provide reliable evidence for investigating the flow field in future research. Figure 11 shows the centrifugal compressor mesh.
Centrifugal compressor mesh.
Performance characteristics of centrifugal compressor.
5 Compressor performance and CFD validation
Figure 11 shows the compressor performance characteristics under four rotational speeds of 70, 80, 90 and 100krpm. The unstable operating range is marked in red for each rotational speed. The last point represents the near-surge, marked in blue. Significant differences were observed in the instability range of the compressor. At 70krpm, the compressor mass flow rate fluctuates suddenly from approximately 0.19kg/s, although the compressor can still operate normally. The instability flow range accounted for 40.34% of the overall operating range. At 80krpm, the mass flow rate fluctuated from approximately 0.21kg/s, and the instability flow range accounted for 38.46% of the overall operating range. At 90krpm, the mass flow rate fluctuated obviously from approximately 0.25kg/s, and the instability flow range accounted for 19.57% of the overall operating range. At 100krpm, the compressor entered the surge state directly from the stable flow. It can be concluded that the instability range at low speeds is significantly wider than that under high-speed conditions. As illustrated in Fig. 12, it is clear that the compressor experiences different types of work states, including stable and unstable operating conditions (RI, stall, and surge). However, it is difficult to visually distinguish between the types and the development of instability based on experimental data.
6 Compressor instability characteristics and evolution
6.1 Pre-Processing
Before data analysis, the dynamic pressure signals obtained from the experimental measurements must be pre-processed. The dynamic pressure signal can be decomposed into constant and pulsating components. The constant component is very large and produces a large amplitude at 0Hz in the spectrum analysis, obscuring other frequency components. Therefore, the average value must be subtracted from the dynamic pressure data, as shown below:
\begin{align}{p_{{\rm{pulsating}}}} = {p_{\exp }} - {p_{{\rm{average}}}} \end{align}
Generally, the BPF component and components less than the BPF are the key research objects. Therefore, it is generally necessary to filter the data to remove the high-frequency components. In this study, the Fourier transform filter is chosen. The low-pass filter blocks all high-frequency components above the BPF and allows only the low-frequency components to pass.
6.2 Instability onset detection based on Fourier transform
The work state of the compressor can be divided into two stages: stable and unstable operating conditions. Under the stable operating conditions, compressor operates well, and there are no strong nonlinear phenomena, such as stall and surge. Therefore, the Fourier transform is sufficient to detect the onset of instability and judge whether the compressor enters an unstable state. Figure 13 shows the frequency spectra of the pressure fluctuation at different circumferential positions from the high mass flow rate point to the low mass flow rate point. The horizontal, longitudinal and vertical axes represent the dimensionless frequencies obtained by dividing the rotor frequency, circumferential angle θ, and static pressure, respectively. The horizontal coordinate 1 represents the RF, whereas 7 represents the BPF. It can be observed that there are almost no significant low-frequency components below BPF under the high mass flow rate (0.36kg/s) and the medium mass flow rate (0.24kg/s) conditions. However, the spectrum characteristics were significantly different at low mass flow rates. At 0.21kg/s, the frequency components below the BPF become very noticeable. The frequency range was mainly concentrated between the 2 RF and BPF. This phenomenon indicates that new and complicated pressure fluctuations occur in the flow field. Therefore, it can be inferred that the compressor entered an unstable working state at the mass flow rate of 0.21kg/s.
Frequency spectra at different mass flow rates at 80krpm.
6.3 Instability characteristics analysis based on VMD
6.3.1 Comparison of Hilbert marginal spectrum and Fourier spectrum
Before VMD, the mode number can be estimated using the Fourier frequency spectrum of the signal x(t). The number of frequency bands in the spectrum diagram can be used as the initial value for choosing the number of decomposition modes by VMD. The amplitude of the BPF component at the 180° circumferential position clearly decreased. Based on the above analysis, the frequency spectrum envelope was extracted at the 180° circumferential position at mass flow rate of 0.21kg/s. As shown in Fig. 14, the local extreme points are marked with a red ball. There were 50 local extreme points; therefore, the initial mode number K was set to 50. The mode number K ranged from 45 to 55. Finally, the optimal K was 50 and the penalty factor was 1.55 f s. Owing to the larger number of IMF components, it is difficult to show them all in the figure. Therefore, only the first five IMF components with a larger Fourier spectrum amplitude from the 20th to 40th revolutions and their corresponding Fourier frequency spectrum results are illustrated in Figs 15 and 16. It can also be observed that these IMF components are aperiodic signals with amplitudes varying with time. These IMFs show significant nonlinear and non-stationary characteristics.
Frequency spectrum peak envelope at 0.21kg/s, at 80krpm.
Decomposition result by VMD.
Fourier spectrum of IMF components.
When the compressor operates at a low mass flow rate, the flow field becomes unstable and complex, thereby leading to flow instabilities (RI, stall or surge). Moreover, the compressor instabilities exhibits a strong nonlinear process. Therefore, the measured static pressure has strong nonlinear characteristics, such as amplitude, average and standard deviation. The traditional method cannot accurately describe this nonlinear problem. Previously, many researchers have used the Fourier transform to analyse the dynamic experimental data, which is nonlinear and non-stationary. However, Huang et al. [Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu43] argued that the Fourier transform could introduce several additional harmonic components to simulate the non-stationary data owing to the linear superposition of trigonometric functions, which may induce spurious harmonic components. Otherwise, it can be observed that the frequency characteristics of the dynamic pressure data show a multi-peak frequency band. It is not clear which important frequencies dominate during the acquisition time. The VMD method is a new, fully intrinsic and adaptive variational method. This method can adaptively decompose the nonlinear and non-stationary signals into their principal modes (IMFs) with different frequencies. Dragomiretskiy [Reference Dragomiretskiy and Zosso50] suggested that the proposed VMD method performs clearly better with less spurious components. To avoid the shortcoming of the Fourier transform, Huang introduced the Hilbert marginal spectrum to analyse nonlinear and non-stationary data, which is suitable for the measurement data of the dynamic experiment of the compressor. The Hilbert marginal spectrum approach is defined as follows:
The time-domain signal x(t) is decomposed into a series of IMFs using VMD method:
\begin{align}x(t){\rm{ = }}\sum\limits_{k = 1}^M {{u_k}} (t) \end{align}
The Hilbert transform result of each IMF
${u_k}$
can be obtained:
\begin{align}H(t) = \frac{1}{\pi }P\int\limits_{ - \infty }^{ + \infty } {\frac{{{u_k}\left( \tau \right)}}{{t - \tau }}} d\tau \end{align}
where P is Cauchy principal value.
An analytical signal can be constructed as below:
\begin{align}{z_k}(t) = {u_k}(t) + iH(t) = {a_k}(t){e^{i{\theta _k}(t)}} \end{align}
\begin{align}{a_k}(t) = \sqrt {{u_k}^2(t) + {H^2}(t)} \end{align}
\begin{align}{\theta _k}(t) = \arctan \left( {\frac{{H(t)}}{{{u_k}(t)}}} \right) \end{align}
where
${a_k}(t)$
and
${\theta _k}(t)$
represent amplitude and phase.
The instantaneous frequency of each IMF
${u_k}$
is calculated as follows:
\begin{align}{\omega _k}(t) = \frac{{d{\theta _k}(t)}}{{dt}} \end{align}
After performing the Hilbert transform on each IMF component, the signal x(t) can be expressed in the following form:
\begin{align}x(t) = {\mathop{\rm Re}\nolimits} \left( {\sum\limits_{k = 1}^M {{a_k}(t){e^{i\int {{\omega _k}(t)dt} }}} } \right) \end{align}
The equation indicates that the amplitude and instantaneous frequency are expressed as functions of time. The amplitude can be contoured on the frequency–time plane, which is called the Hilbert amplitude spectrum:
\begin{align}H\left( {\omega ,t} \right) = {\mathop{\rm Re}\nolimits} \left( {\sum\limits_{k = 1}^M {{a_k}(t){e^{i\int {{\omega _k}(t)dt} }}} } \right) \end{align}
On this basis, a marginal spectrum was introduced to investigate the instability process. The detailed definition is as follows:
\begin{align}h(\omega) = \int\limits_0^T {H({\omega ,t})} dt \end{align}
The Hilbert marginal spectrum is the time-domain integral of the Hilbert spectrum, which depicts the cumulative amplitudes of the whole data. In other words, the Hilbert marginal spectrum offers a measure of the total amplitude contribution from each frequency value. The frequency component with the large cumulative amplitude is dominant in the multicomponent signal. Huang [Reference Huang, Shen, Long, Wu, Shih, Zheng, Yen, Tung and Liu43] provides more detailed comparisons with the Fourier transform, which shows that the marginal spectrum is a method with higher resolution and is suitable for nonlinear and non-stationary signals. Yue et al. [Reference Yue, Wang, Wei, Zhang and Wang49] applied the Hilbert marginal spectrum to analyse the stall and compared it with the Fourier transform to accurately identify the frequency of the stall signal.
Figure 17 shows the comparison between the Hilbert marginal spectrum and Fourier spectrum at the 180° circumferential position at 0.21kg/s. In order to distinguish the signals, the component below RF is called the low-frequency component, and the component between RF and BPF is called the high-frequency component. Generally, the low-frequency component signal represents the stall, while the high-frequency component signal represents the RI. For the high-frequency signal, it was also called high-frequency stall [Reference Zheng and Liu55] or high-frequency pressure oscillation [Reference Zheng, Sun, Kawakubo and Tamaki56]. At 0.21kg/s, there are two distinct frequency components: 3 RF and 4.1 RF. However, there are multiple components in the Fourier spectrum, and it is difficult to determine which components are dominant in the origin signal. Compared to the Fourier frequency spectrum, the Hilbert marginal spectrum shows the dominant frequency except BPF, more clearly and has a better ability to analyse the instability signals. Although there are still many frequency components in the Hilbert marginal spectrum, they show the dominant frequencies more clearly to a great extent.
Comparison between the Hilbert marginal spectrum and Fourier spectrum at 180° circumferential position at 0.21kg/s, 80krpm.
Figure 18 illustrates the comparison of the Hilbert marginal spectra of the 180° circumferential position at 0.21 and 0.20kg/s. At 0.20kg/s, the amplitude of the high-frequency components increases, and the amplitude of the BPF decreases significantly. This phenomenon indicates that compressor has gone into the RI state at 0.21kg/s and the RI is further enhanced at 0.20kg/s. In addition, it can also be observed that there are no significant low-frequency components at 0.21 and 0.20kg/s.
Hilbert marginal spectrum at 180° circumferential position at 0.20kg/s, 80krpm.
6.3.2 Evolution and circumferential nonuniformity characteristics of flow instability
Figure 19 shows the Hilbert marginal spectrum at different circumferential positions at 0.18, 0.16, 0.14, and 0.12kg/s. At 0.18kg/s, in addition to BPF, there is a broadband frequency hump composed of high-frequency components, representing the RI. The frequency band indicates that under this operating condition, RI is the primary type of instability. The broadband frequency hump of the RI is mainly concentrated between 2 RF and 5 RF. However, there is a small difference in the frequency spectrum characteristics at different circumferential positions. At the 180° circumferential position, the amplitude of the BPF is distinctly different from that of other circumferential positions. Compared with 0.20kg/s, the amplitude of the BPF significantly decreases, reaching a value lower than the maximum amplitude of the frequency band of the RI. Meanwhile, the RI amplitude increases significantly. On the other hand, at the 60°, 120°, 180°, and 240° circumferential position, there are low-frequency components, representing the stall occurrence. However, there is no obvious stall frequency at the 0° and 300° circumferential position. It can be concluded that the instability type of the compressor has an obvious circumferential nonuniformity. At this operating point, RI is the main instability phenomenon, whereas a slight stall occurs in the partial circumferential range.
Hilbert marginal spectrum at different circumferential positions at 80krpm.
At 0.16kg/s, in addition to BPF, the broadband frequency hump of RI is mainly concentrated between RF and 3 RF. Compared to the operating condition with a mass flow rate of 0.18kg/s, the broadband frequency hump moves from high frequency to low frequency, and the Hilbert marginal spectrum amplitude decreases significantly. At all circumferential positions, low-frequency stall frequencies were observed. In addition, the stall frequency at different circumferential positions has a noticeable circumferential nonuniformity. The amplitude of the stall frequency at 240° circumferential positions is larger than that of the RI. At this operating point, the most obvious feature is that the RI and stall coexist more clearly.
When the compressor operates at the mass flow rate of 0.14kg/s, it can be observed that the low-frequency components increase significantly, whereas the high-frequency components become very small. In the entire circumferential range, the amplitude of the low-frequency stall was greater than that of the high-frequency RI. At this operating point, the compressor instability was dominated by the low-frequency stall. The compressor enters the stall state completely.
When the mass flow further decreases to 0.12kg/s near the surge point, it can be observed that the low-frequency stall completely dominates, and the high-frequency components almost disappear in the entire circumferential positions. Meanwhile, it can be seen that the amplitude of BPF decreases. At the 180° circumferential position, the amplitude of the BPF is far less than that of the stall. The amplitude of the stall is almost the same as that of the BPF at the 240° circumferential position. It can be inferred that if the mass flow further decreases, the amplitude of the BPF might be less than that of the stall. In the experiment, when the mass flow rate decreased further, the compressor quickly entered the surge state. It can be concluded that deep stall first occurs at these two positions, which is related to the impending surge. According to the experiments, the compressor surge may be triggered by a deep stall in the circumferential range of 180 to 240°.
Based on the above analysis, both the RI and stall exhibit circumferential nouniformity for the centrifugal compressor. At 0.21kg/s, the RI occurs in the entire circumferential range. At 0.18kg/s, the compressor first enters the stall condition in the circumferential range of 60° to 240°, while the other circumferential position still operates under the RI condition. As the mass flow rate decreases, the deep stall first occurs in the circumferential range of 120° to 240°. Generally, according to the above analysis, as the mass flow rate decreases, the instability process of the compressor can be summarised as follows: RI → RI & stall (coexistence stage) → stall (slight stall and deep stall) → surge. The process is clearly marked on the curve of the performance characteristics in Fig. 20.
Instability process at 80krpm.
Casing wall static pressure distributions at 80krpm.
7 Discussions
The cause of circumferential nouniformity of RI and stall was deeply explored by performing unsteady simulations at 0.26, 0.24, and 0.21kg/s (RI). The static pressure distribution of the casing wall is shown in Fig. 21. The dimensionless static pressure was obtained by dividing the static pressure by the circumferential average. The horizontal axis represents the circumferential angle, and the longitudinal axis represents the dimensionless streamwise position. In this figure, 0 and 1 represent the impeller inlet and outlet, respectively. Thus, the circumferentially nonuniform static pressure distribution is clearly indicated. A significant circumferential non-uniform distribution of the static pressure was observed around the casing wall. This nonuniform distribution phenomenon could only be caused by the asymmetric structure of the volute tongue. Previous studies [Reference Yang, Wang, Tong, Yang and Li57–Reference Zhang, Yang, Wang, Chen and Qi59] have discussed the effect of the volute tongue on the static pressure distribution in a centrifugal compressor. The airflow is divided by the volute tongue and forced to recirculate around the volute at a low flow rate. Therefore, the volute tongue acts as a diffuser at a low mass flow rate, thereby forming a high static pressure region. Eventually, the volute tongue causes an adverse pressure wave. In Fig. 21, it can be observed that there is a high static pressure strip representing the pressure wave, marked with a black curve with an arrow. The high static pressure strip extends to the impeller inlet, resulting in a circumferentially nonuniform distribution at the inlet of the impeller. A high static pressure region HB is seen at the upstream of the impeller at these three operating points. The high-pressure region HB is approximately concentrated between 120° and 300°. This region is most affected by reverse pressure wave. As the mass flow rate decreases, at 0.24kg/s, a low-pressure region LA between 120° and 270° occurs near the LE of the blade. The low-pressure region LA results in a pressure decrease caused by the blade rotation, which is manifested as a significant decrease in the amplitude corresponding to the BPF, as shown in Fig. 13. At 0.21kg/s, some high-pressure spots occur upstream of the inlet. This indicates that the degree of circumferential nonuniformity becomes more obvious. The high-pressure region HB was concentrated between approximately 120° and 240°. Based on the analysis in Section 6.3, RI first occurs at this operating point. It can be concluded that increasing the circumferential nonuniformity leads to the occurrence of RI.
Circumferential averaged streamline distribution of the compressor at 80krpm.
Figure 22 illustrates the movement of the tip leakage vortex at three different operating points. As the mass flow rate decreases, the rolling up of a shear layer between the tip leakage flow and the main flow forms a tip leakage vortex. At 0.26kg/s, the tip leakage vortex is located at the tip region near the impeller inlet. Under the action of the reverse pressure gradient, the tip leakage vortex moves to the inlet of the impeller passage. At 0.24kg/s, the vortex size increases; meanwhile, the reach of the vortex expands to the upstream of the impeller inlet, and the vortex core moves towards the blade leading edge. At this time, the low-pressure region LA is caused by the tip leakage vortex, as shown in Fig. 21. At 0.21kg/s, it can be seen that the vortex core migrates to the lower span and the reach of the vortex further extends upstream. During this period, RI first occurs at the 180° circumferential position in the compressor. This phenomenon indicates that the movement of the tip leakage vortex is closely related to the occurrence of the RI. If the mass flow rate further decreases, RI and stall just first occur in the high-pressure region HB caused by the volute tongue. With the further development of the stall, the 180° and 240° circumferential positions successively enter a deep stall, and finally, a surge occurs in the centrifugal compressor. Based on the above analysis, it can be summarised that the high static pressure region induced by the volute tongue determines the circumferential occurrence position of the RI and stall at the impeller inlet.
8 Conclusions
In this study, the evolution and circumferential nonuniformity characteristics of the flow instabilities in a centrifugal compressor with a volute were investigated based on dynamic pressure experiments. During the throttling process, six dynamic pressure sensors were installed evenly at the inlet of the compressor to obtain the dynamic pressure at the different circumferential positions. Both classical Fourier transform and advanced signal processing technique (VMD) were introduced to shed light on instability phenomena at 80krpm. The following conclusions were drawn:
-
(1) A Fourier transform was used to determine whether the compressor entered a state of instability. The VMD method and Hilbert marginal spectrum were applied to analyse the pressure signals under different instability conditions. Using this method, not only could the additional spurious harmonic components in the Fourier transform be avoided to the greatest extent, but also the frequency resolution for the nonlinear and non-stationary instability signals could be improved. This combined method can achieve better results when dealing with dynamic pressure signals.
-
(2) During the throttling process, the centrifugal compressor experienced different instability stages: RI stage, coexistence stage of RI and stall, and stall stage. In the RI stage (0.21–0.20kg/s), the Hilbert marginal spectrum amplitude gradually increased, which indicated that strength of RI gradually increased. In the coexistence stage of RI and stall (0.18–0.16kg/s), the Hilbert marginal spectrum amplitude of RI gradually decreased, and the corresponding broadband frequency hump moved from high to low frequency, which represented a gradual decrease in the RI strength. In this stage, the Hilbert marginal spectrum amplitude of the BPF at the 180° circumferential position was less than the maximum amplitude of the RI. In the stall stage (0.14–0.12kg/s), the RI almost disappeared. The compressor instability was dominated by low-frequency stall. Eventually, the compressor entered the surge state.
-
(3) During the throttling process, the instability process exhibited an obvious circumferential nonuniformity phenomenon in the centrifugal compressor. The RI first occurred at the 180° circumferential position. The stall first occurred in the circumferential range of 60°–240°. As the mass flow rate decreased, a stall occurred in the entire circumferential range. Compared to the other circumferential positions, the BPF amplitude decreased the most between 180° and 240°, which in turn triggered a surge.
-
(4) The initial occurrence position of the instability phenomena in the centrifugal compressor was determined by the asymmetric volute. The reverse pressure wave induced by the volute tongue led to a high static pressure region upstream of the impeller inlet. Under the combined action of the high static pressure region and reverse pressure gradient, the tip leakage vortex migrated to the inlet and the lower span. When the tip leakage vortex greatly exceeded the blade leading edge, it triggers instability.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant No. 51736001) and was undertaken at the Turbomachinery Institute under the Beijing Institute of Technology, China.
Declaration of competing interest
The authors have no conflict of interest declare.
