A method for monitoring and early warning of rotating stall in an axial flow compressor

Abstract

The application provides a rotating stall monitoring and early warning method for an axial flow compressor. Through a dynamic pressure sensor in a dynamic data acquisition system, pressure pulsation signals of different positions on the compressor wall surface in the circumferential and axial directions from a steady state to a stall process are collected. The collected signals are decomposed by VMD to obtain intrinsic mode function components. IMF components containing stall precursors and pre-stall precursor characteristic frequencies in a frequency band range are selected. The energy of the relevant component signals IMF is calculated respectively. Meanwhile, the energy of the original signal in the frequency band from 0 to the blade passing frequency is calculated in the frequency domain. The energy proportion is calculated. The abnormal threshold of the energy proportion distribution at each time is calculated to early warn the rotating stall of the compressor. The method can achieve the effect of timely early warning of stall precursors and pre-stall precursors.

Description

Technical Field

[0001] This invention relates to the field of active control of compressor stability, and more particularly to a method for monitoring and early warning of rotating stall in an axial compressor. Background Technology

[0002] With the development of high-performance aero-engine and gas turbine technologies, compressors are required to have high pressure ratios, high efficiency, and high reliability. However, these three requirements are mutually restrictive; higher pressure ratios and efficiencies generally mean that compressors are more prone to instability. Currently, existing control strategies for compressor instability mainly include active and passive control technologies. The main idea of ​​active control technology is to detect stall precursors and pre-stall precursors from signals before instability occurs and apply reverse disturbances to suppress instability and broaden the stable operating range. Stall precursors can be divided into modal stall precursors and spike stall precursors, while pre-stall precursors refer to minute disturbances that occur before the stall precursors appear.

[0003] Currently, the main detection techniques for stall precursors and pre-stall precursors include correlation analysis and wavelet transform. Correlation analysis only analyzes pressure signals in the time domain, and the correlation coefficient results are easily affected by random interference. Wavelet transform can perform multi-resolution analysis of signals in the time and frequency domains, but its effectiveness depends on the selection of wavelet basis functions and it cannot adaptively decompose the signal. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention proposes a method for monitoring and early warning of rotating stall in axial compressors based on Variational Mode Decomposition (VMD). The specific technical solution adopted in this invention is as follows:

[0005] Step 1: Collect pressure pulsation signals from steady state to stall process at different positions of the compressor wall in the circumferential and axial directions using the dynamic pressure sensor in the dynamic data acquisition system.

[0006] Step 2: Perform VMD decomposition on the acquired signal to obtain the intrinsic mode function components IMF1, IMF2, ... IMF k , where k is the number of components to be decomposed.

[0007] Step 3: Select the IMF component within the frequency band that contains the characteristic frequencies of stall precursor and pre-stall precursor, and denote it as IMF. s IMF sp .

[0008] Step 4: Calculate the IMF of each component signal. s IMF sp Energy E s and Esp Simultaneously, the energy E of the original signal (the signal before decomposition) within the frequency domain, from 0 to the Blade Pass Frequency (BPF), is calculated. bpf .

[0009] Step 5: Calculate the energy percentage K s K sp :

[0010]

[0011]

[0012] Step 6: Calculate K at each time point using the box plot method. s K sp The outlier threshold of the distribution, if K s Exceeding the upper limit threshold UL s , or K sp Exceeding the upper limit threshold UL sp This will provide an early warning of compressor rotational stall.

[0013] In step 1 above, to avoid frequency aliasing and signal amplitude distortion, the sampling frequency should be greater than 10 times the BPF.

[0014] Optionally, in step 1 above, the circumferential measuring points of the pulsating pressure are arranged on the compressor casing wall, and are evenly distributed circumferentially at the rotor inlet section or at a section 15% of the distance from the blade tip to the leading edge. The axial measuring points of the pulsating pressure are arranged along the chord length direction on the casing wall above the rotor, and the total length of the measuring points covers the chord length of the rotor blade tip.

[0015] In step 3 above, the stall precursor characteristic frequency is generally lower than the rotor rotation frequency F. r A certain characteristic frequency; the characteristic frequency of the pre-stall precursor is generally a certain characteristic frequency lower than that of the BPF.

[0016] In step 3 above, if a suitable IMF component cannot be found, the number of components k and the penalty factor α in the VMD algorithm are changed, and the signal is decomposed again.

[0017] In step 6 above, the upper limit threshold UL s UL sp Determined by the following formula:

[0018] UL = Q3 + 1.5IQR

[0019] In the formula, Q3 represents K under normal conditions. s K sp The upper quartile of the distribution; Q1 is the normal state of K s K spThe lower quartile of the distribution; IQR is the interquartile range, calculated using the following formula:

[0020] IQR = Q3 - Q1.

[0021] The beneficial effects that this invention can bring are as follows:

[0022] 1) The technical solution provided by this invention can be used to monitor and warn of axial compressor rotational stall.

[0023] 2) The technical solution provided by the present invention can adaptively decompose the dynamic pressure signal of the compressor.

[0024] 3) By performing VMD decomposition on the pressure pulsation signal, the energy ratio of the IMF components containing the characteristic frequencies of stall precursors and pre-stall precursors to the original signal is calculated. Combined with statistical methods, the effect of timely early warning of stall precursors and pre-stall precursors is achieved. Attached Figure Description

[0025] Figure 1 A flowchart of the axial compressor rotational stall early warning method of the present invention is provided.

[0026] Figure 2 A schematic diagram of the structure of a three-stage axial compressor with inlet and outlet guide vanes and a schematic diagram of the stall warning measurement point arrangement are given in an embodiment of the present invention. Detailed Implementation

[0027] To more clearly illustrate the technical solution of this invention, this embodiment takes a three-stage axial flow compressor with inlet and outlet guide vanes as an example to introduce the specific implementation process of the technical solution of this invention. It should be understood that the embodiments described herein are merely for explaining the invention and are not intended to limit the invention; the invention can also be applied to other embodiments.

[0028] This invention proposes a method for monitoring and early warning of rotating stall in an axial compressor. The specific process is as follows (see attached diagram). Figure 1 ;

[0029] Step 1: Collect pressure pulsation signals from steady state to stall process at different positions of the compressor wall in the circumferential and axial directions using the dynamic pressure sensor in the dynamic data acquisition system.

[0030] like Figure 2 As shown, four dynamic pressure sensors are evenly arranged circumferentially at the inlet section AA of the first row of moving blades. In addition, ten dynamic pressure sensors are equally spaced along the axial direction at the tip of the first row of moving blades, with the measuring point range covering 1.5 times the chord length of the tip.

[0031] In this embodiment, compressor stall can be induced by gradually closing the intake throttle valve to the stall boundary.

[0032] In step 1 above, to avoid frequency aliasing and signal amplitude distortion, the sampling frequency should be greater than 10 times the blade passage frequency. In this embodiment, the design speed of the 3-stage axial compressor is 9600 rpm, and the number of blades in the first row is 16; therefore, the sampling frequency of the dynamic pressure sensor should be greater than 25.6 kHz. While ensuring a sufficient sampling frequency, the amount of data processed can be reduced by downsampling.

[0033] Optionally, if a phase-locked signal is acquired during the experiment, the time-domain signal can be converted to the angular domain using resampling technology, and then order tracking can be implemented to perform subsequent processing and analysis of the signal from the order perspective.

[0034] Step 2: Perform VMD decomposition on the acquired signal to obtain the intrinsic mode function components IMF1, IMF2, ... IMF k , where k is the number of components to be decomposed.

[0035] The VMD algorithm decomposes a signal into multiple intrinsic mode components (AM-FM signals) with different center frequencies and finite bandwidths. It iteratively updates the center frequency and bandwidth by constructing a variational optimization problem to minimize the sum of the estimated bandwidths of each mode component, thereby determining the frequency center and bandwidth of each mode component.

[0036] Specifically, the constrained variational optimization problem in the VMD algorithm can be expressed as:

[0037]

[0038]

[0039] In the formula, For gradient calculation, t is time, j is the imaginary unit, δ(t) is the Dirac function, and u k (t) is the intrinsic mode function, ω k denoted as , where is the center frequency of each intrinsic mode function component, and f(t) is the original dynamic pressure signal before decomposition.

[0040] By introducing a quadratic penalty factor α and the Lagrange operator λ(t), the constrained variational problem is transformed into an unconstrained variational problem, and the constructed Lagrange expression is:

[0041]

[0042] Solve the above equation using the alternating direction multiplier method, and iteratively update {u}. k}、{ω k}, λ until the convergence condition is met.

[0043] Step 3: Select the IMF component within the frequency band that contains the characteristic frequencies of stall precursor and pre-stall precursor, and denote it as IMF. s IMF sp .

[0044] Stall precursor frequency is generally lower than the rotor rotation frequency F. r A certain characteristic frequency, mainly corresponding to the circumferential propagation of the stall cluster, can be initially selected as 0.3F. r ~0.7F r The characteristic frequency of pre-stall precursors is generally a certain characteristic frequency below BPF, which is related to the scale of small disturbances, and can be initially selected as 0.2BPF to 0.8BPF.

[0045] If a suitable IMF component cannot be found, change the number of modal components k and the penalty factor α in the VMD algorithm, and re-decompose the signal.

[0046] Optionally, the adaptive values ​​of the number of modal components k and the penalty factor α can be achieved through optimization algorithms, correlation coefficients, information entropy, and other methods.

[0047] Step 4: Calculate the IMF of each component signal. s IMF sp Energy E s and E sp And in the frequency domain, the energy E of the original signal (the signal before decomposition) within the band from 0 to the Blade Pass Frequency (BPF) is calculated. bpf .

[0048] According to Passevar's theorem, the energy of a signal is equal in the time and frequency domains. Therefore, the energy of the component signal can be directly calculated in the time domain, and the energy of the original signal in the (0, BPF) band can be calculated in the frequency domain. When the signal energy is characterized by its effective value, E... s E sp and E bpf It can be obtained through the following formula:

[0049]

[0050]

[0051]

[0052] In the formula, f s (k), f sp (k) represents the component signal IMF. s IMF spX(k) represents the amplitude of the kth data point in the time domain, X(k) represents the amplitude of the kth data point in the frequency domain of the original signal, N1 and N2 represent the number of data points of the component signal in the entire time domain, and N3 represents the number of data points of the original signal in the frequency domain (0, BPF).

[0053] Step 5: Calculate the energy percentage K s K sp :

[0054]

[0055]

[0056] Step 6: Calculate K at each time point using the box plot method. s K sp The outlier threshold of the distribution, if K s Exceeding the upper limit threshold UL s , or K sp Exceeding the upper limit threshold UL sp This provides an early warning for compressor rotating stall. The upper limit threshold UL is used for this purpose. s UL sp Determined by the following formula:

[0057] UL = Q3 + 1.5IQR

[0058] In the formula, Q3 represents K under normal conditions. s K sp The upper quartile of the distribution; Q1 is the normal state of K s K sp The lower quartile of the distribution; IQR is the interquartile range, calculated using the following formula:

[0059] IQR = Q3 - Q1

[0060] Furthermore, after determining the warning, the warning signal can be output in the form of a switch signal and connected to the control system to perform audible and visual alarms and activate the actuators of the stabilization device (such as tip jet, adjustable guide vanes and stationary vanes), thereby achieving active control of compressor rotational stall.

[0061] It should be noted that, without departing from the technical solution of this invention, those skilled in the art can make equivalent substitutions or modifications to the relevant technical features, and the embodiments after such substitutions or modifications still fall within the protection scope of this invention.

Claims

1. A method for monitoring and early warning of rotating stall in an axial compressor, characterized in that, The method includes the following steps: Step 1: Collect pressure pulsation signals from steady state to stall process at different positions of the compressor wall in the circumferential and axial directions using the dynamic pressure sensor in the dynamic data acquisition system; Step 2: Perform VMD decomposition on the acquired signal to obtain the intrinsic mode function components IMF1, IMF2, ... IMF k , where k is the number of components to be decomposed; Step 3: Select the IMF component within the frequency band that contains the characteristic frequencies of stall precursor and pre-stall precursor, and denote it as IMF. s IMF sp ; Step 4: Calculate the IMF of each component signal. s IMF sp Energy E s and E sp Simultaneously, the original signal, i.e. the signal before decomposition, is calculated in the frequency domain, from 0 to the energy E within the frequency band through which the blade passes. bpf ; Step 5: Calculate the energy percentage K s K sp : Step 6: Calculate K at each time point using the box plot method. s K sp The outlier threshold of the distribution, if K s Exceeding the upper limit threshold UL s , or K sp Exceeding the upper limit threshold UL sp This will provide an early warning of compressor rotational stall.

2. The method according to claim 1, characterized in that, In step 1 above, to avoid frequency aliasing and signal amplitude distortion, the sampling frequency should be greater than 10 times the blade passage frequency.

3. The method according to claim 1, characterized in that, In step 1 above, the circumferential measuring points of the pulsating pressure are arranged on the compressor casing wall, and are evenly distributed circumferentially at the rotor inlet section or at a section 15% of the distance from the blade tip to the leading edge; the axial measuring points of the pulsating pressure are arranged along the chord length direction on the casing wall above the rotor, and the total length of the measuring points covers the chord length of the rotor blade tip.

4. The method according to claim 1, characterized in that, In step 2, the constrained variational optimization problem in the VMD algorithm is expressed as: In the formula, For gradient calculation, t is time, j is the imaginary unit, δ(t) is the Dirac function, and u k (t) is the intrinsic mode function, ω k Here, f(t) represents the center frequency of each intrinsic mode function component, and f(t) represents the original dynamic pressure signal before decomposition. By introducing a quadratic penalty factor α and the Lagrange operator λ(t), the constrained variational problem is transformed into an unconstrained variational problem, and the constructed Lagrange expression is: Solve the above equation using the alternating direction multiplier method, and iteratively update {u}. k }、{ω k }, λ until the convergence condition is met.

5. The method according to claim 1, characterized in that, In step 3 above, the stall precursor characteristic frequency is generally lower than the rotor rotation frequency F. r A certain characteristic frequency; the characteristic frequency of the pre-stall precursor is generally a certain characteristic frequency lower than the blade passing frequency.

6. The method according to claim 4, characterized in that, In step 3 above, if a suitable IMF component cannot be found, the number of components k and the penalty factor α in the VMD algorithm are changed, and the signal is decomposed again.

7. The method according to claim 1, characterized in that, In step 6 above, the upper limit threshold UL s UL sp Determined by the following formula: UL = Q3 + 1.5IQR; In the formula, Q3 represents K under normal conditions. s K sp The upper quartile of the distribution; Q1 is the normal state of K s K sp The lower quartile of the distribution; IQR is the interquartile range, calculated using the following formula: IQR = Q3 - Q1.

Images