CN115711738B - Method for extracting fault characteristic frequency of main bearing of aeroengine in service environment - Google Patents
Method for extracting fault characteristic frequency of main bearing of aeroengine in service environment Download PDFInfo
- Publication number
- CN115711738B CN115711738B CN202211301763.0A CN202211301763A CN115711738B CN 115711738 B CN115711738 B CN 115711738B CN 202211301763 A CN202211301763 A CN 202211301763A CN 115711738 B CN115711738 B CN 115711738B
- Authority
- CN
- China
- Prior art keywords
- frequency
- feature
- fault
- signal
- signals
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 34
- 238000000605 extraction Methods 0.000 claims abstract description 23
- 238000010183 spectrum analysis Methods 0.000 claims abstract description 11
- 230000001133 acceleration Effects 0.000 claims abstract description 9
- 238000001228 spectrum Methods 0.000 claims description 33
- 238000001914 filtration Methods 0.000 claims description 19
- 230000000737 periodic effect Effects 0.000 claims description 19
- 238000005316 response function Methods 0.000 claims description 9
- 230000009467 reduction Effects 0.000 claims description 8
- 230000003044 adaptive effect Effects 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 7
- 230000003595 spectral effect Effects 0.000 claims description 7
- 238000004458 analytical method Methods 0.000 claims description 6
- 238000010606 normalization Methods 0.000 claims description 6
- 238000012937 correction Methods 0.000 claims description 3
- 230000009466 transformation Effects 0.000 claims description 3
- 238000003745 diagnosis Methods 0.000 abstract description 7
- 238000005516 engineering process Methods 0.000 abstract description 2
- 238000005096 rolling process Methods 0.000 description 14
- 238000004901 spalling Methods 0.000 description 7
- 238000000926 separation method Methods 0.000 description 6
- 238000012360 testing method Methods 0.000 description 5
- 230000005540 biological transmission Effects 0.000 description 4
- 238000013461 design Methods 0.000 description 3
- 238000012544 monitoring process Methods 0.000 description 3
- 230000008859 change Effects 0.000 description 2
- PMGQWSIVQFOFOQ-YKVZVUFRSA-N clemastine fumarate Chemical compound OC(=O)\C=C\C(O)=O.CN1CCC[C@@H]1CCO[C@@](C)(C=1C=CC(Cl)=CC=1)C1=CC=CC=C1 PMGQWSIVQFOFOQ-YKVZVUFRSA-N 0.000 description 2
- 230000003111 delayed effect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 238000002485 combustion reaction Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 230000011664 signaling Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
A method for extracting the fault characteristic frequency of main bearing of aeroengine in service environment. The invention discloses a novel method for extracting characteristic frequency of weak fault signals of a main bearing of a military aircraft engine in a strong noise service environment. The method mainly comprises the following steps: aiming at the collected vibration acceleration signals of the aeroengine case measuring points, the weak random signal extraction method technology is adopted to separate the signals, discrete components and random components of the signals are obtained, and dimensionless characteristic quantities reflecting the bearing fault characteristic frequency local/overall characteristic significance level are constructed. The method can effectively extract random components in the signals, filter interference components of discrete parts such as rotating speed, gear meshing frequency and the like in square envelope spectrum analysis, improve the corresponding amplitude of a frequency band where the fault frequency is located, inhibit the influence of clutter on fault diagnosis, improve the accuracy of main bearing fault feature extraction, and provide an important method and an effective path for main bearing fault feature extraction of the aeroengine.
Description
Technical Field
The invention belongs to the field of fault diagnosis of mechanical systems of aeroengines, and particularly relates to a novel method for extracting characteristic frequencies of weak fault signals of main bearings of military aeroengines in a strong noise service environment.
Background
Under the actual working condition, the main bearing of the aeroengine is influenced by the radial action of centrifugal load, the contact angles of the rolling bodies and the inner ring and the outer ring deviate from the initial values, so that the actual fault characteristic frequency deviates from the theoretical fault characteristic frequency, and the extraction of the fault characteristics of the main bearing is not facilitated. In addition, after the vibration signal of the main bearing of the aero-engine passes through a complex transmission path, the fault characteristic signal becomes extremely weak, and a plurality of rotor rotating speed frequency and frequency multiplication thereof, gear engagement frequency and frequency multiplication thereof of an accessory transmission case, rotor blade passing frequency and frequency multiplication thereof, engine air-driven and combustion and other interference frequency components are doped in the signal, so that the fault characteristic frequency components of the main bearing are easily confused and are difficult to identify. Therefore, the periodic signals with discrete frequency spectrum characteristics are required to be removed, weak main bearing fault characteristic signals are stripped from complex modulated original signals, the actual contact angle of the rolling body and the inner ring/outer ring is judged, and the fault characteristics are effectively extracted.
The basic principle of the weak random signal extraction method is that a Frequency Response Function (FRF) between a main vibration sequence and a delay vibration sequence is searched, and then a self-reference self-adaptive separation filter is obtained through inverse Fourier transform, so that effective separation of bearing fault characteristic signals and rest signals is realized based on the filter. Therefore, the invention provides a main bearing fault feature extraction method suitable for aeroengine case measuring point signals. The weak random signal extraction technology is utilized to strip the part irrelevant to the main bearing fault diagnosis in the case measuring point signal layer by layer, and the part is combined with the autocorrelation analysis and the square envelope spectrum to construct the fault feature significance, and the fault feature significance is applied to the main bearing fault complete machine test of the engine, and the verification result shows that the main bearing fault feature extraction can be effectively realized, so that an important method and an effective path are provided for the main bearing fault diagnosis of the aeroengine.
Disclosure of Invention
The invention aims at solving the technical problems in the background technology, and provides a method for extracting the main bearing fault characteristics of an aeroengine, which can provide an important method and an effective way for continuously improving the use safety level of a military aeroengine.
In order to achieve the above purpose, the invention is realized by the following technical scheme:
s1, aiming at the collected vibration acceleration signals of the measuring points of the engine case of the aeroengine, adopting a weak random signal extraction technology to separate the vibration acceleration signals of the aeroengine, and obtaining discrete components and random components of the signals. The discrete components comprise periodic signals such as rotating speed frequency and frequency multiplication thereof, gear meshing frequency and frequency multiplication thereof, and the random part contains bearing fault characteristic frequency, noise signals and the like.
S2, square enveloping, autocorrelation noise reduction and spectrum analysis are carried out on the obtained random components, and bearing fault characteristic frequency and corresponding contact angles are extracted.
S3, based on the feature frequency extracted by the S2, constructing dimensionless feature quantity reflecting the local/whole feature significance level of the bearing fault feature frequency, and providing a feature significance calculating method based on an S-shaped function, so that the normalization of the dimensionless feature is realized, and the bearing fault feature is accurately identified.
The step S1 specifically includes:
s1-1: dividing a main vibration sequence and a delay vibration sequence: x is an original vibration sequence, and is divided into sub-signals with lengths of 2L+τ, wherein L is the order of a self-reference adaptive filter to be generated, and τ is the delay between a main vibration sequence and a delay vibration sequence. Delayed vibration sequenceAnd a main vibration sequence x k (n) can be expressed as:
x k (n)＝w L (n)x[n-L+k(2L+τ)]，0＜n≤L，0＜k≤K (2)
wherein ,wL (n) is a window function of length L, a Parzen window function is used, K is the number of sub-signal segments, and K segments are used.
S1-2: constructing a self-reference adaptive filter: the purpose of the filter is to achieve a slaveTo x k (n) optimal prediction, such that +.>On the premise of keeping the periodic part of the X-ray detector unchanged, the X-ray detector is used for detecting the X-ray k (n) optimal prediction, non-periodic part and wideband noise will be from->Filtering out. The self-referenced adaptive optimal filter frequency response function H (f) that achieves the minimum mean square prediction error can be expressed as:
wherein ,SUV (f) Is the cross power spectrum of signals U and V, S UU (f) For the self-power spectrum of the signal U, p denotes the periodic component and r denotes the aperiodic component. Whilex k (n)＝p k (n)+r k (n), H (f) can be expressed in simplified terms:
wherein ,delay vibration sequence for the kth sub-signal>Performing Fourier transform with length of M, X k，M (f) Main vibration sequence x for k-th segment sub-signal k (n) performing a fourier transform of length M.And M is greater than or equal to L, the effective length of the frequency response function is equal to the length of the signal participating in Fourier transform, when M is greater than L, the excess part is +.>Or x k After (n) 0 is added, the frequency spectrum resolution of the bearing fault is reduced, so that the invention adopts M=L to obtain a filter frequency response function H (f) with the length of L, and the influence of different data lengths on the characteristic frequency is eliminated.
S1-3: frequency domain filtering and phase correction: dividing x into sub-signals x [ n+ (K ' -1) ] L with the lengths of L, wherein n is more than 0 and less than or equal to L, and K ' is the number of sub-signal segments and is K ' in total. The frequency domain filtering of the sub-signals is sequentially performed by H (f), which can be expressed as:
X k′，L (f)＝H(f)·X k′,L (f) (5)
wherein ,Xk′，L (f) For the k' th segment filtering result, X k′，L (f) Is the fourier transform of the k' th segment length L sub-signal. In order to ensure that the phase of the obtained filtered signal is the same as that of the original signal, the filtered frequency spectrum needs to be reconstructed according to the amplitude-frequency characteristic after filtering and the phase-frequency characteristic before filtering, and the method can be expressed as follows:
Re＝|X k′，L (f)|cos{Ψ[X k′，L (f)]} (7)
Im＝|X k′，L (f)|sin{Ψ[X k′，L (f)]} (7)
X k′，L (f)＝Re+i·Im (8)
wherein ,Xk′，L (f) Is X k′，L (f) The reconstructed spectrum, |is the magnitude of the phase, and ψ is the phase of the phase.
S1-4: obtaining the periodic component x p And an aperiodic component x r : for the obtained X k′，L (f) Performing inverse Fourier transform to obtain periodic component x with length L p，k′ And with x [ n+ (k' -1) L]Subtracting to obtain the aperiodic component x with length L r，k′ . Repeating the step S1-3 to obtain x corresponding to 0 < K' < K p，k′ and xr，k′ Head-to-tail phaseIs connected to obtain x p and xr 。
It should be noted that the choice of the time delay factor τ and the filter order L has a great influence on the separation performance of the filter. Wherein the time delay factor tau should be slightly larger than the correlation length of the random component, and the filter order L is at least 10-20 times the minimum frequency selected to remove the discrete component.
The step S2 specifically includes:
s2-1: square enveloping is carried out on the random component obtained in the step S1;
s2-2: performing autocorrelation noise reduction on a signal obtained by square envelope;
s2-3: carrying out spectrum analysis on the noise-reduced signal;
s2-4: extracting bearing fault characteristic frequency, specifically comprising:
assuming that the contact angle has a transformation range of alpha min To alpha max Between them, calculate the square envelope spectrum Y SE The searching range of the failure characteristic frequency of the middle bearing is [ f ] min ，f max ]. Wherein Z is the number of rolling bodies, D is the diameter of the rolling bodies, D is the pitch diameter of the bearing, alpha is the contact angle and f i For the rotation frequency of the inner ring, f o For the rotation frequency of the outer ring, the outer ring of the rolling bearing is generally fixed, and the inner ring rotates, i.e. |f i -f o |＝f i 。
In the square envelope spectrum Y SE In the method, the maximum value of spectral lines in the searching range of the fault characteristic frequency of the bearing is the corresponding amplitude value of the fault characteristic frequencyThe corresponding frequency is the fault characteristic frequency f oc The corresponding angle is the contact angle alpha, < + >>And α is represented as:
the step S3 specifically includes:
s3-1: solving the bearing fault characteristic frequency f oc Average value of nearby spectral lines, assuming f oc The corresponding bandwidth is delta f, the average value calculation range is delta f, and the square envelope spectrum interval is f ε ，f oc The average of the nearby spectral lines is expressed as:
f i ∈[f oc -Δf/2-δf/2，f oc -Δf/2]，f j ∈[f oc +Δf/2，f oc +Δf/2+δf/2]
s3-2: constructing a reflection bearing fault characteristic frequency f oc Dimensionless feature quantity S of local feature saliency loc ：
S3-3: constructing a reflection bearing fault characteristic frequency f oc Dimensionless feature quantity S of overall feature significance gen Assume that the frequency range of the desired analysis is [ f ]' min ，f′ max ]S is then gen Can be expressed as:
wherein ,indicating that the first 5 maxima are taken.
S3-4: obtaining a characteristic frequency f capable of reflecting the bearing faults at the same time oc Is characterized by (1)Dimensionless feature quantity S of the degree of the writing and the overall feature significance:
S＝S loc ·S gen (15)
s3-5: in order to achieve normalization of dimensionless features and accurately represent the prominence and the significance degree of vibration amplitude values at the fault feature frequency of the bearing relative to surrounding points, a feature significance value calculation method based on an S-shaped function is provided. The function expression is as follows:
normalizing the obtained dimensionless characteristic quantity S to obtain the final bearing fault characteristic significance S o ，S o ∈[O，1]。
Compared with the prior art, the invention has the following technical effects:
1) The method can effectively extract random components in the signals, filter interference components of discrete parts such as rotating speed, gear meshing frequency and the like in square envelope spectrum analysis, improve the corresponding amplitude of the frequency band of the fault frequency and inhibit the influence of clutter on fault diagnosis.
2) The method can effectively improve the accuracy of extracting the main bearing fault characteristics, accurately identify the characteristic frequency of the early spalling fault of the main bearing of the aeroengine under the influence of a complex transmission path and strong noise interference, and determine the contact angle between the main bearing rolling body and the inner/outer ring under the actual working condition.
3) The bearing fault feature normalized dimensionless feature quantity constructed by the method can effectively represent the fault feature of the bearing, and can be used as an effective index and an advantageous criterion for the fault diagnosis and on-line monitoring of the main bearing of the aeroengine.
Drawings
FIG. 1 is a diagram showing the construction method of a delay vibration sequence and a main vibration sequence in a weak random signal extraction method
FIG. 2 is a technical flow chart of the present invention
FIG. 3 is a physical diagram of an initial spalling failure of a bearing
FIG. 4 is a time domain waveform and spectrogram of the original signal
FIG. 5 is a partial spectrogram of the original signal
FIG. 6 is a graph of the square envelope of the original signal
FIG. 7 shows a discrete component time domain waveform and spectrogram
FIG. 8 is a graph of a time domain waveform and a spectrum of the separated random component
FIG. 9 is a square envelope spectrum after auto-correlation of random components
FIG. 10 is an S-shaped graph
FIG. 11 is a graph of the saliency of the original signal versus the random component features
FIG. 12 is a plot of the contact angle distribution of the original signal versus the random component
FIG. 13 is a graph showing the statistical comparison of contact angles of the original signals and the random components
Detailed Description
The invention further provides a method for extracting fault characteristics of the whole vibration main bearing of the aeroengine based on the weak random signal extraction technology by referring to the accompanying drawings.
S1, aiming at the collected vibration acceleration signals of the measuring points of the engine case of the aeroengine, separating the vibration acceleration signals of the aeroengine by adopting a weak random signal extraction technology to obtain discrete components and random components of the signals. The discrete components comprise periodic signals such as rotating speed frequency and frequency multiplication thereof, gear meshing frequency and frequency multiplication thereof, and the random part contains bearing fault characteristic frequency, noise signals and the like.
S2, square enveloping, autocorrelation noise reduction and spectrum analysis are carried out on the obtained random components, and bearing fault characteristic frequency and corresponding contact angles are extracted.
S3, based on the feature frequency extracted by the S2, constructing dimensionless feature quantity reflecting the local/whole feature significance level of the bearing fault feature frequency, and providing a feature significance calculating method based on an S-shaped function, so that the normalization of the dimensionless feature is realized, and the bearing fault feature is accurately identified.
The step S1 specifically includes:
s1-1: dividing a main vibration sequence and a delay vibration sequence: x is an original vibration sequence, and is divided into sub-signals with lengths of 2L+τ, wherein L is the order of a self-reference adaptive filter to be generated, and τ is the delay between a main vibration sequence and a delay vibration sequence. Delayed vibration sequenceAnd a main vibration sequence x k (n) can be expressed as:
x k (n)＝w L (n)x[n-L+k(2L+τ)]，0＜n≤L，0＜k≤K (2)
wherein ,wL (n) is a window function of length L, a Parzen window function is used, K is the number of sub-signal segments, and K segments are used. The dividing method of the main vibration sequence and the delay vibration sequence is as shown in fig. 1.
S1-2: constructing a self-reference adaptive filter: the purpose of the filter is to achieve a slaveTo x k (n) optimal prediction, such that +.>On the premise of keeping the periodic part of the X-ray detector unchanged, the X-ray detector is used for detecting the X-ray k (n) optimal prediction, non-periodic part and wideband noise will be from->Filtering out. The self-referenced adaptive optimal filter frequency response function H (f) that achieves the minimum mean square prediction error can be expressed as:
wherein ,SUV (f) Is the cross power spectrum of signals U and V, S UU (f) For the self-power spectrum of signal U, p represents the periodic component and r represents the aperiodic component. Whilex k (n)＝p k (n)+r k (n), H (f) can be expressed in simplified terms:
wherein ,delay vibration sequence for the kth sub-signal>Performing Fourier transform with length of M, X k，M (f) Main vibration sequence x for k-th segment sub-signal k (n) performing a fourier transform of length M. And M is greater than or equal to L, the effective length of the frequency response function is equal to the length of the signal participating in Fourier transform, when M is greater than L, the greater part should be +.>Or x k After (n), 0 is complemented, and the filter frequency response function H (f) with the length of L is obtained by adopting M=L.
S1-3: frequency domain filtering and phase correction: dividing x into sub-signals x (n+ (K ' -1) L) with the lengths of L, wherein n is more than 0 and less than or equal to L, and K ' is the number of sub-signal segments and is K ' in total. The frequency domain filtering of the sub-signals is sequentially performed by H (f), which can be expressed as:
X k′，L (f)＝H(f)·X k′，L (f) (5)
wherein ,Xk′，L (f) For the k' th segment filtering result, X k′，L (f) Is the fourier transform of the k' th segment length L sub-signal. In order to ensure that the phase of the obtained filtered signal is the same as that of the original signal, the filtered frequency spectrum needs to be reconstructed according to the amplitude-frequency characteristic after filtering and the phase-frequency characteristic before filtering, and the method can be expressed as follows:
Re＝|X k′，L (f)|cos{Ψ[X k′，L (f)]} (6)
Im＝|X k′，L (f)|sin{Ψ[X k′，L (f)]} (7)
X k′，L (f)＝Re+i·Im (8)
wherein ,Xk′，L (f) Is X k′，L (f) The reconstructed spectrum, |is the magnitude of the phase, and ψ is the phase of the phase.
S1-4: obtaining the periodic component x p And an aperiodic component x r : for the obtained X k′，L (f) Performing inverse Fourier transform to obtain periodic component x with length L p，k′ And subtracting from x (n+ (k' -1) L) to obtain an aperiodic component x with length L r，k′ . Repeating the step S1-3 to obtain x corresponding to 0 < K' < K p，k′ and xr，k′ End to obtain x p and xr 。
It should be noted that the choice of the time delay factor τ and the filter order L has a great influence on the separation performance of the filter. Wherein the time delay factor tau should be slightly larger than the correlation length of the random component, and the filter order L is at least 10-20 times the minimum frequency selected to remove the discrete component.
The step S2 specifically includes:
s2-1: and (3) square enveloping the random component obtained in the step S1.
S2-2: and carrying out autocorrelation noise reduction on the signal obtained by square envelope.
S2-3: and carrying out spectrum analysis on the noise-reduced signal.
S2-4: and extracting the fault characteristic frequency of the bearing. The specific operation is as follows:
assuming that the contact angle has a transformation range of alpha min To alpha max Between them, calculate the square envelope spectrum Y SE The searching range of the failure characteristic frequency of the middle bearing is [ f ] min ，f max ]. Wherein Z is the number of rolling bodies, D is the diameter of the rolling bodies, D is the pitch diameter of the bearing, alpha is the contact angle and f i For the rotation frequency of the inner ring, f o For the rotation frequency of the outer ring, the outer ring of the rolling bearing is generally fixed, and the inner ring rotates, i.e. |f i -f o |＝f i 。
In the square envelope spectrum Y SE In the method, the maximum value of spectral lines in the searching range of the fault characteristic frequency of the bearing is the corresponding amplitude value of the fault characteristic frequencyThe corresponding frequency is the fault characteristic frequency f oc The corresponding angle is the contact angle alpha, < + >>And α is represented as:
the step S3 specifically includes:
s3-1: solving the bearing fault characteristic frequency f oc Average value of nearby spectral lines, assuming f oc The corresponding bandwidth is delta f, the average value calculation range is delta f, and the square envelope spectrum interval is f ε ，f oc The average of the nearby spectral lines is expressed as:
f i ∈[f oc -Δf/2-δf/2，f oc -Δf/2]，f j ∈[f oc +Δf/2，f oc +Δf/2+δf/2]
s3-2: constructing a dimensionless feature quantity S reflecting the local feature significance of the bearing fault feature frequency foc loc ：
S3-3: constructing a reflection bearing fault characteristic frequency f oc Dimensionless feature quantity S of overall feature significance gen Assume that the frequency range of the desired analysis is [ f ]' min ，f′ max ]S is then gen Can be expressed as:
wherein ,indicating that the first 5 maxima are taken.
S3-4: obtaining a characteristic frequency f capable of reflecting the bearing faults at the same time oc Dimensionless feature quantity S of local feature saliency and global feature saliency:
S＝S loc ·S gen (15)
s3-5: in order to achieve normalization of dimensionless features and accurately represent the prominence and the significance degree of vibration amplitude values at the fault feature frequency of the bearing relative to surrounding points, a feature significance value calculation method based on an S-shaped function is provided. The function expression is as follows:
normalizing the obtained dimensionless characteristic quantity S to obtain the final outer ring fault characteristic significance S o ，S o ∈[0，1]。
In order to illustrate the applicability of the method of the invention to the use environment of a real engine, the test data of a complete machine rack of the main bearing spalling fault of a certain domestic military aviation engine is selected, the three-pivot bearing is an angular contact ball bearing with an outer ring spalling fault, the outer ring spalling size is about 3 mm by 5 mm, and the method belongs to early spalling, as shown in figure 3. The test was run for 155 hours and 34 minutes and eventually stopped by vibration and swarf signaling alarms due to rapid increases in the oil temperature. The vibration acceleration sensor is arranged at an airborne vibration monitoring measuring point of the intermediate casing, and the sampling frequency is 200kHz.
According to the method for monitoring the fault evolution state of the rolling bearing of the aeroengine newly proposed by the inventor, the bearing is in early spalling fault in 108.9 hours before the test, so that the vibration signal in 108.9 hours before the test is selected, the engine is limited to be in a middle stable state, and the characteristic frequency f of the outer ring fault can be calculated at the moment corresponding to N1 rotating speed of 8051rpm, N2 rotating speed of 14417rpm and initial contact angle of 39 DEG of the bearing o 2155Hz. The original waveform of a set of vibration signals and the spectrum characteristics thereof are shown in fig. 4, the local spectrum (1800 Hz-2300 Hz) is shown in fig. 5, and the spectrum after square envelope and autocorrelation is shown in fig. 6. It can be seen that the characteristic frequency component of the outer ring fault of the main bearing is difficult to directly judge from the influence of noise such as a transmission path, strong noise, a periodic signal and the like due to frequency multiplication. In addition, theoretical analysis shows that the bearing roller is influenced by centrifugal radial load, so that the contact angle between the bearing roller and the outer ring is reduced, the contact angle between the bearing roller and the inner ring is increased, and the contact angle is changed along with the change of working conditions, so that the extraction of the fault characteristic frequency of the outer ring is more difficult.
Firstly, a weak random signal extraction method is carried out on the signals, discrete interference components in vibration acceleration signals of the whole machine of the aeroengine are filtered, random components containing main bearing fault characteristic frequencies are obtained, and on the basis, the main bearing outer ring fault characteristics are identified and extracted according to the contact angle change rule of the rolling bodies and the outer ring. According to the method described in step S1, τ=100 (slightly greater than the random component correlation length f s /f o Approximately 92.81), l=2700 (greater than 20 times the N1 rotational frequency 8051/60 approximately 134.18 Hz), the discrete components and spectra obtained by separation are shown in fig. 7, and the random components and spectra obtained by separation are shown in fig. 8. From the above, a large amount of frequency multiplication in the original frequency spectrum is effectively separated into discrete components by a weak random signal extraction technology, and discrete frequency components, such as rotor rotation speed frequency, gear engagement frequency, frequency multiplication and the like, which have influence on the subsequent main bearing outer ring fault feature extraction are filtered, so that the weak random signal extraction technology effectively separates the discrete components and the random components of the signals, and the interference is reduced for the subsequent main bearing outer ring fault feature extraction.
Subsequently, referring to the method described in step S2, the characteristic frequency of the bearing outer ring fault is extracted, and the square envelope spectrum of the random signal is obtained through the methods of square envelope, autocorrelation noise reduction and spectrum analysis, as shown in fig. 9. Comparing the spectrum of the original signal of fig. 7, there are energy concentrations of 2133Hz, 2146Hz, 2164Hz frequency components on both sides of the theoretical outer ring failure frequency 2155Hz, corresponding to contact angles of 32.3 °,36.7 °,41.7 °, respectively. From the theory of bearing dynamics, it is known that the contact angle of the rolling elements with the outer ring will be smaller than the design angle, i.e. should be smaller than the design value of 39 ° due to the centrifugal force of the rolling elements. Thus the possibility of 2164Hz being the actual outer ring failure frequency and 41.7 ° being the actual contact angle can be excluded. The other 2133Hz and 2146Hz cannot be accurately judged, the 2146Hz is estimated to be the actual outer ring fault frequency more likely only by the fact that the corresponding amplitude of the 2146Hz frequency component is larger than that of the 2133Hz frequency component, the possibility that 2133Hz is the theoretical outer ring fault characteristic frequency cannot be eliminated, and the actual outer ring fault frequency cannot be effectively and accurately extracted. As can be seen from the spectrogram of the random signal in fig. 9, 2133Hz (16 times the frequency of the N1 rotation speed) and 2164Hz (9 times the frequency of the N2 rotation speed) frequency components are effectively filtered, and 2146Hz can be judged to be the actual outer ring fault characteristic frequency, and 36.7 ° is the actual contact angle. Therefore, the method can effectively extract random components in the signals, filter interference components of discrete parts in square envelope spectrum analysis, improve the corresponding amplitude of the frequency band of the fault frequency and inhibit the influence of clutter on fault diagnosis.
It is worth noting that 2146Hz is closer to 16 times of the N1 rotation speed than 2133Hz, but after being processed by the weak random signal extraction technology, the 2146Hz frequency component is reserved, which means that in this state, the outer ring fault frequency is close to 16 times of the N1 rotation speed, so that the outer ring fault frequency is effectively excited, but the frequency spectrum cannot be effectively identified and extracted due to harmonic interference of the rotation speed frequency and the meshing frequency. After the processing of the weak random signal extraction technology, harmonic components are effectively filtered, the bearing fault frequency is reserved, and finally the actual outer ring fault frequency can be effectively and accurately extracted, so that the rule that the contact angle of the rolling body and the outer ring is smaller than the design angle under the actual working condition is met and verified.
In order to further verify the bearing fault feature normalized dimensionless feature quantity provided by the invention, vibration data of 61 hours 43 minutes to 68 hours 58 hours are selected, large-state steady-state vibration acceleration data with the N2 rotating speed being more than 13900rpm is limited to perform feature significance and contact angle calculation, the parameter k of the S-shaped function is 0.1, the parameter a is 50, and the S-shaped curve is shown in figure 10. The calculated original signal and random component feature saliency points are compared with each other as shown in fig. 11, the outer ring feature saliency of which the saliency is greater than 0.8 corresponds to the outer ring fault feature frequency, the corresponding contact angle is calculated, the contact angle distribution points are compared with each other as shown in fig. 12, and the contact angle statistics pair is as shown in fig. 13.
As can be seen from fig. 11 (a) and fig. 12 (a), the significance of the failure feature of the outer ring calculated from the original signal is relatively concentrated, which means that the energy distribution of the failure frequency of the outer ring extracted is more concentrated, and the contact angle is mainly concentrated at about 41.7 °. However, the above analysis has been performed, the calculated contact angle does not conform to the actual law, and the frequency is concentrated at 9 times the frequency of the N2 rotation speed, which is an interference component. The random signal after the DRS processing is effectively screened, and as can be seen from fig. 12 (b) and fig. 13 (b), the contact angle is mainly concentrated between 28 ° and 36 °, and interference components such as rotation speed frequency multiplication are effectively filtered, so that the method provided by the invention can effectively extract the characteristic frequency of the outer ring fault and calculate the corresponding contact angle.
Finally, it should be noted that: the embodiments described above are only for illustrating the technical solution of the present invention, and are not limiting; although the invention has been described in detail with reference to the foregoing embodiments, it will be understood by those of ordinary skill in the art that: the technical scheme described in the foregoing embodiments can be modified or some or all of the technical features thereof can be replaced with equivalents; such modifications and substitutions do not depart from the spirit of the invention.
Claims (5)
1. The characteristic frequency extraction method for the weak fault signal of the main bearing of the military aircraft engine in the strong noise service environment is characterized by comprising the following steps:
s1, measuring vibration acceleration signals aiming at collected aeroengine case measuring points, and separating the signals by adopting a weak random signal extraction technology to obtain discrete components and random components of the signals; the discrete component comprises a rotating speed frequency and a frequency multiplication thereof, a gear meshing frequency and a frequency multiplication signal thereof, and the random component comprises a bearing fault characteristic frequency and a noise signal;
s2, square envelope, autocorrelation noise reduction and spectrum analysis are carried out on the obtained random components, and bearing fault characteristic frequency is extracted;
s3, constructing dimensionless feature quantities reflecting the local/overall feature significance level of the bearing fault feature frequency based on the feature frequency extracted by the S2, and adopting a feature significance calculation method based on an S-shaped function to realize normalization of the dimensionless feature so as to accurately identify the bearing fault feature; the specific operation of the step S1 is as follows:
s1-1: dividing a main vibration sequence and a delay vibration sequence: x is an original vibration sequence, and is divided into sub-signals with lengths of 2L+τ, wherein L is the order of a self-reference self-adaptive filter to be generated, τ is the delay between a main vibration sequence and a delay vibration sequence, and the delay vibration sequenceAnd a main vibration sequence x k (n) can be expressed as:
wherein ,wL (n) is a window function with the length of L, a Parzen window function is adopted, K is the number of sub-signal segments, and K segments are used;
s1-2: constructing a self-reference self-adaptive optimal filter:the purpose of the filter is to achieve a slaveTo x k (n) optimal prediction, such that +.>On the premise of keeping the periodic part of the X-ray detector unchanged, the X-ray detector is used for detecting the X-ray k (n) optimal prediction, non-periodic part and wideband noise will be from->The self-reference adaptive optimal filter frequency response function H (f) that achieves the minimum mean square prediction error can be expressed as:
wherein ,SUV (f) Is the cross power spectrum of signals U and V, S UU (f) For the self-power spectrum of signal U, p represents the periodic component and r represents the non-periodic component;
s1-3: frequency domain filtering and phase correction: dividing x into sub-signals x (n+ (K ' -1) L) with the lengths of L, wherein n is more than 0 and less than or equal to L, K ' is the number of sub-signal segments, K ' is the total number of segments, and H (f) is utilized to sequentially carry out frequency domain filtering on the sub-signals, and the frequency domain filtering is expressed as follows:
wherein ,for the k' th segment filtering result, X k′,L (f) For fourier transformation of the k' th sub-signal with length L, to ensure that the phase of the obtained filtered signal is the same as that of the original signal, the spectrum after filtering needs to be reconstructed according to the amplitude-frequency characteristic after filtering and the phase-frequency characteristic before filtering, and the method is expressed as follows:
wherein ,is->The reconstructed frequency spectrum, |is the amplitude value, and ψ is the phase value;
s1-4: obtaining deterministic component x p And a random component x r : for the obtainedPerforming inverse Fourier transform to obtain periodic component x with length L p,k′ And do nothing to->Subtracting to obtain the aperiodic component x with length L r,k′ Repeating the step S1-3 to obtain x corresponding to K' being more than 0 and less than or equal to K p,k′ and xr,k′ End to obtain x p and xr ；
The step S3 specifically comprises the following steps:
s3-1: solving the bearing fault characteristic frequency f oc Average value of nearby spectral lines, assuming f oc The corresponding bandwidth is delta f, the average value calculation range is delta f, and the square envelope spectrum interval is f ε ，f oc Nearby spectrumThe average value of the line is expressed as:
f i ∈[f oc -Δf/2-δf/2,f oc -Δf/2],f j ∈[f oc +Δf/2,f oc +Δf/2+δf/2]
s3-2: constructing a reflection bearing fault characteristic frequency f oc Dimensionless feature quantity S of local feature saliency loc ：
S3-3: constructing a reflection bearing fault characteristic frequency f oc Dimensionless feature quantity S of overall feature significance gen Assume that the frequency range of the desired analysis is [ f ]' min ,f′ max ]S is then gen Can be expressed as:
wherein ,representing taking the first 5 maxima;
s3-4: obtaining a characteristic frequency f capable of reflecting the bearing faults at the same time oc Dimensionless feature quantity S of local feature saliency and global feature saliency:
S＝S loc ·S gen (15)
s3-5: in order to achieve normalization of dimensionless features and accurately characterize the prominence and significance of vibration amplitude at the bearing fault feature frequency relative to surrounding points, a feature significance value based on an S-shaped function is calculated by using the following formula:
normalizing the obtained dimensionless characteristic quantity S to obtain the final bearing fault characteristic significance S o ,S o ∈[0,1]。
2. The method for extracting feature frequency of fault signal according to claim 1, wherein in step S1-2,x k (n)＝p k (n)+r k (n), H (f) may be represented as:
wherein ,delay vibration sequence for the kth sub-signal>Performing Fourier transform with length of M, X k,M (f) Main vibration sequence x for k-th segment sub-signal k (n) performing a fourier transform of length M.
3. The method for extracting feature frequency of fault signal according to claim 2, wherein,
and obtaining a filter frequency response function H (f) with the length L by adopting M=L so as to eliminate the influence of different data lengths on the characteristic frequency.
4. The method for extracting feature frequency of fault signal according to claim 1, wherein,
the time delay factor tau should be greater than the correlation length of the random component and the filter order L should be at least 10-20 times the minimum frequency selected to remove the discrete component.
5. The fault signature frequency extraction method according to claim 1, wherein the specific operation of step 2 is:
s2-1: square enveloping is carried out on the random component obtained in the step S1;
s2-2: performing autocorrelation noise reduction on a signal obtained by square envelope;
s2-3: carrying out spectrum analysis on the signal after autocorrelation and noise reduction;
s2-4: bearing failure feature frequencies are extracted from the spectral analysis.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211301763.0A CN115711738B (en) | 2022-10-24 | 2022-10-24 | Method for extracting fault characteristic frequency of main bearing of aeroengine in service environment |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202211301763.0A CN115711738B (en) | 2022-10-24 | 2022-10-24 | Method for extracting fault characteristic frequency of main bearing of aeroengine in service environment |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN115711738A CN115711738A (en) | 2023-02-24 |
| CN115711738B true CN115711738B (en) | 2023-10-31 |
Family
ID=85231539
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202211301763.0A Active CN115711738B (en) | 2022-10-24 | 2022-10-24 | Method for extracting fault characteristic frequency of main bearing of aeroengine in service environment |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN115711738B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117851873B (en) * | 2024-03-07 | 2024-05-28 | 唐智科技湖南发展有限公司 | Bearing running state evaluation method and system based on dynamic contact angle |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107631877A (en) * | 2017-08-11 | 2018-01-26 | 南京航空航天大学 | A kind of rolling bearing fault collaborative diagnosis method for casing vibration signal |
| CN110470475A (en) * | 2019-09-04 | 2019-11-19 | 中国人民解放军空军工程大学航空机务士官学校 | A kind of aero-engine intershaft bearing early-stage weak fault diagnostic method |
| CN113326782A (en) * | 2021-06-01 | 2021-08-31 | 西安交通大学 | Rolling bearing fault feature automatic extraction method based on envelope spectrum form fitting |
| CN114486263A (en) * | 2022-02-15 | 2022-05-13 | 浙江大学 | Noise reduction and demodulation method for vibration signal of rolling bearing of rotary machine |
| CN114705426A (en) * | 2022-01-06 | 2022-07-05 | 杭州爱华仪器有限公司 | Early fault diagnosis method for rolling bearing |
-
2022
- 2022-10-24 CN CN202211301763.0A patent/CN115711738B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107631877A (en) * | 2017-08-11 | 2018-01-26 | 南京航空航天大学 | A kind of rolling bearing fault collaborative diagnosis method for casing vibration signal |
| CN110470475A (en) * | 2019-09-04 | 2019-11-19 | 中国人民解放军空军工程大学航空机务士官学校 | A kind of aero-engine intershaft bearing early-stage weak fault diagnostic method |
| CN113326782A (en) * | 2021-06-01 | 2021-08-31 | 西安交通大学 | Rolling bearing fault feature automatic extraction method based on envelope spectrum form fitting |
| CN114705426A (en) * | 2022-01-06 | 2022-07-05 | 杭州爱华仪器有限公司 | Early fault diagnosis method for rolling bearing |
| CN114486263A (en) * | 2022-02-15 | 2022-05-13 | 浙江大学 | Noise reduction and demodulation method for vibration signal of rolling bearing of rotary machine |
Also Published As
| Publication number | Publication date |
|---|---|
| CN115711738A (en) | 2023-02-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110470475B (en) | Early weak fault diagnosis method for intermediate bearing of aircraft engine | |
| CN106769033B (en) | Variable speed rolling bearing fault recognition methods based on order envelope time-frequency energy spectrum | |
| CN109682601A (en) | The initial failure recognition methods of rolling bearing under a kind of variable speed operating condition | |
| US7027953B2 (en) | Method and system for diagnostics and prognostics of a mechanical system | |
| US7409854B2 (en) | Method and apparatus for determining an operating status of a turbine engine | |
| CN104236908B (en) | Combined slicing bearing fault diagnosis method on basis of MID (modulation intensity distribution) algorithm | |
| CN115711738B (en) | Method for extracting fault characteristic frequency of main bearing of aeroengine in service environment | |
| CN108151869A (en) | A kind of mechanical oscillation characteristic index extracting method, system and device | |
| Guo et al. | Tooth root crack detection of planet and sun gears based on resonance demodulation and vibration separation | |
| Zhao et al. | Rolling element bearing instantaneous rotational frequency estimation based on EMD soft-thresholding denoising and instantaneous fault characteristic frequency | |
| CN116659860B (en) | Novel method for monitoring main bearing fault evolution of aeroengine in service environment | |
| Lin et al. | A review and strategy for the diagnosis of speed-varying machinery | |
| Chen et al. | Proportional selection scheme: A frequency band division tool for rolling element bearing diagnostics | |
| Li et al. | Application of a Method of Identifiying Instantaneous Shaft Speed from Spectrum in Aeroengine Vibration Analysis | |
| Babiker et al. | Initial fault time estimation of rolling element bearing by backtracking strategy, improved VMD and infogram | |
| Lu et al. | Detection of weak fault using sparse empirical wavelet transform for cyclic fault | |
| Loughlin et al. | Conditional moments analysis of transients with application to helicopter fault data | |
| Chen et al. | Improvement on IESFOgram for demodulation band determination in the rolling element bearings diagnosis | |
| Qin et al. | Adaptive fast chirplet transform and its application into rolling bearing fault diagnosis under time-varying speed condition | |
| Zhu et al. | Adaptive combined HOEO based fault feature extraction method for rolling element bearing under variable speed condition | |
| He et al. | Weak fault detection method of rolling bearing based on testing signal far away from fault source | |
| Mauricio et al. | Cyclostationary-based tools for bearing diagnostics | |
| CN116644280A (en) | High-frequency vibration diagnosis method for failure of main bearing of engine | |
| CN117109923A (en) | Rolling bearing fault diagnosis method and system | |
| CN116150585A (en) | Rotary machine fault diagnosis method based on product envelope spectrum |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |