Welcome to my official account "fault diagnosis and python learning"
1. Data introduction
Including inner ring, outer ring, rolling element and normal data, which are one respectively.
1730_12k_0.007-InnerRace: inner ring fault
1730_12k_0.007-OuterRace3: outer ring fault
1730_12k_0.014-Ball: rolling element failure
1730_48k_Normal: normal data
For students who do not understand CWRU bearing data set, see here:
CWRU dataset introduction phase 1
CWRU dataset introduction phase 2
CWRU dataset introduction phase 3
CWRU dataset introduction phase 3
2. Load CWRU inner ring fault data
Next, we will analyze the fault data of bearing inner race. The original data is a mat file, which is a matlab file. Define a function to read the data
def data_acquision(FilePath):
"""
fun: from cwru mat File reading acceleration data
param file_path: mat File absolute path
return accl_data: Acceleration data, array type
"""
data = scio.loadmat(file_path) # Load mat data
data_key_list = list(data.keys()) # The mat file is of dictionary type. Get all the keys of the dictionary and convert them to list type
accl_key = data_key_list[3] # Get'X108_DE_time'
accl_data = data[accl_key].flatten() # Get'x108_ DE_ The value corresponding to time 'is the vibration acceleration signal, and the two-dimensional array is expanded into a one-dimensional array
return accl_data
data_ Acquisition (FilePath) input the parameter FilePath and output a 1-dimensional array data. Here is a demonstration
file_path = r'E:/graduate student/pytorch/CSDN code/fault_diagnosis_signal_processing/Part 4-Envelope spectrum/1730_12k_0.007-InnerRace.mat' xt = data_acquision(file_path) plt.figure(figsize=(12,3)) plt.plot(xt) print(xt)
Output results:
[ 0.22269856 0.09323776 -0.14651649 ... -0.36125573 0.31138814 0.17055689]
3. Hilbert Demodulation (envelope spectrum) analysis
Hilbert Demodulation, also known as envelope spectrum analysis.
3.1 Hilbert Huang Transform
set up
x
(
t
)
x(t)
x(t) is a real-time domain signal, and its Hilbert transform is defined as:
h
(
t
)
=
1
π
∫
−
∞
+
∞
x
(
τ
)
t
−
τ
d
τ
=
x
(
t
)
∗
1
π
t
h(t)=\frac{1}{\pi} \int_{-\infty}^{+\infty} \frac{x(\tau)}{t-\tau} \mathrm{d} \tau=x(t) * \frac{1}{\pi t}
h(t)=π1∫−∞+∞t−τx(τ)dτ=x(t)∗πt1
Then the original signal
x
(
t
)
x(t)
x(t) and its Hilbert transform signal
h
(
t
)
h(t)
h(t) can construct a new analytic signal
z
(
t
)
z(t)
z(t):
z
(
t
)
=
x
(
t
)
+
j
h
(
t
)
=
a
(
t
)
e
j
φ
t
z(t)=x(t)+j h(t)=a(t) e^{j \varphi t}
z(t)=x(t)+jh(t)=a(t)ejφt
# step1: do Hilbert transform ht = fftpack.hilbert(xt) print(ht)
Output results:
[-0.02520403 -0.28707983 -0.00610516 ... 0.1100125 0.22821944 -0.11203138]
Output at this time h ( t ) h(t) h(t) is the analytic signal a ( t ) a(t) Imaginary part coefficient of a(t)
yes z ( t ) z(t) Take the mode of z(t) and get its amplitude a ( t ) : a(t): a(t):
a ( t ) = ∣ z ( t ) ∣ = x 2 ( t ) + h 2 ( t ) {a(t)=|z(t)|=\sqrt{x^{2}(t)+h^{2}(t)}} a(t)=∣z(t)∣=x2(t)+h2(t)
Note: a ( t ) a(t) a(t) is envelope signal, also known as analytic signal
3.2 obtaining envelope spectrum
sampling_rate = 12000 am = np.fft.fft(at) # fft is used to obtain the amplitude of at after Hilbert transform am = np.abs(am) # Calculate the absolute value of the amplitude (at this time, the absolute value is very large) am = am/len(am)*2 am = am[0: int(len(am)/2)] freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # Obtain fft frequency, which includes positive frequency and negative frequency freq = freq[0:int(len(freq)/2)] # Get positive frequency plt.plot(freq, am)
Observations:
(1) Where the frequency is 0Hz, the amplitude is relatively large
(2) In the low-frequency part, it seems to see the first and second harmonic
3.3 remove DC component
The amplitude at 0Hz is relatively high, which makes the amplitude of other frequencies low and inconvenient to observe. This phenomenon is called DC component, the method of removing DC component, y = y-mean(y)
sampling_rate = 12000 at = at - np.mean(at) # To DC component am = np.fft.fft(at) # fft is used to obtain the amplitude of at after Hilbert transform am = np.abs(am) # Calculate the absolute value of the amplitude (at this time, the absolute value is very large) am = am/len(am)*2 am = am[0: int(len(am)/2)] freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # Obtain fft frequency, which includes positive frequency and negative frequency freq = freq[0:int(len(freq)/2)] # Get positive frequency plt.plot(freq, am)
Output results:
sampling_rate = 12000 at = at - np.mean(at) # To DC component am = np.fft.fft(at) # fft is used to obtain the amplitude of at after Hilbert transform am = np.abs(am) # Calculate the absolute value of the amplitude (at this time, the absolute value is very large) am = am/len(am)*2 am = am[0: int(len(am)/2)] freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # Obtain fft frequency, which includes positive frequency and negative frequency freq = freq[0:int(len(freq)/2)] # Get positive frequency plt.figure(figsize=(12,3)) plt.plot(freq, am) plt.xlim(0,500)
Output results:
4. Calculate fault characteristic frequency
Inner ring fault characteristic frequency: F B P F I = n f r 2 ( 1 + d D cos α ) F_{\mathrm{BPFI}}=\frac{n f_{r}}{2}\left(1+\frac{d}{D} \cos \alpha\right) FBPFI=2nfr(1+Ddcosα)
Outer ring fault characteristic frequency: F B P F O = n f r 2 ( 1 − d D cos α ) F_{\mathrm{BPFO}}=\frac{n f_{r}}{2}\left(1-\frac{d}{D} \cos \alpha\right) FBPFO=2nfr(1−Ddcosα)
Rolling element fault characteristic frequency: F B S F = D f r 2 d [ 1 − ( d D cos α ) 2 ] F_{\mathrm{BSF}}=\frac{D f_{r}}{2 d}\left[1-\left(\frac{d}{D} \cos \alpha\right)^{2}\right] FBSF=2dDfr[1−(Ddcosα)2]
n n n: Number of scrolls, f r f_{r} fr: shaft speed d d d: Diameter of ball (sub) D D D: Bearing pitch diameter
Bearing model: 6205-2RSL JME SKF deep groove ball bearing
d
d
d=7.94mm,
D
D
D=39.04mm,
α
\alpha
α=0,
n
n
n=9
4.1 define a bearing fault characteristic frequency calculation function
For convenience, a calculation function of bearing fault characteristic frequency is defined, which only needs to input parameters d d d, D D D, α \alpha α, n n n and f r f_{r} fr is enough
def bearing_fault_freq_cal(n, d, D, alpha, fr=None):
'''
Basic description:
Calculate the fault characteristic frequency of rolling bearings
Detailed description:
Enter 4 parameters n, fr, d, D, alpha
return C_bpfi, C_bpfo, C_bsf, C_ftf, fr
Inner ring Outer ring Speed of needle cage
Parameters
----------
n: integer
The number of roller element
fr: float(r/min)
Rotational speed
d: float(mm)
roller element diameter
D: float(mm)
pitch diameter of bearing
alpha: float(°)
contact angle
fr:: float(r/min)
rotational speed
Returns
-------
BPFI: float(Hz)
Inner race-way fault frequency
BPFO: float(Hz)
Outer race-way fault frequency
BSF: float(Hz)
Ball fault frequency
FTF: float(Hz)
Cage frequency
'''
C_bpfi = n*(1/2)*(1+d/D*np.math.cos(alpha))
C_bpfo = n*(1/2)*(1-(d/D)*np.math.cos(alpha))
C_bsf = D*(1/(2*d))*(1-np.square(d/D*np.math.cos(alpha)))
C_ftf = (1/2)*(1-(d/D)*np.math.cos(alpha))
if fr!=None:
return C_bpfi*fr/60, C_bpfo*fr/60, C_bsf*fr/60, C_ftf*fr/60, fr/60
else:
return C_bpfi, C_bpfo, C_bsf, C_ftf, fr
Next, calculate the fault characteristic frequency of CWRU at 1730rpm
bpfi, bpfo, bsf, ftf, fr = bearing_fault_freq_cal(n=9, alpha=0, d=7.94, D=39.04, fr=1730)
print('Inner ring fault characteristic frequency',bpfi)
print('Outer ring fault characteristic frequency',bpfo)
print('Rolling element fault characteristic frequency',bsf)
print(ftf)
print(fr)
5. Verification of theoretical fault characteristic frequency and actual fault characteristic frequency
sampling_rate = 12000 at = at - np.mean(at) # To DC component am = np.fft.fft(at) # fft is used to obtain the amplitude of at after Hilbert transform am = np.abs(am) # Calculate the absolute value of the amplitude (at this time, the absolute value is very large) am = am/len(am)*2 am = am[0: int(len(am)/2)] freq = np.fft.fftfreq(len(at), d=1 / sampling_rate) # Obtain fft frequency, which includes positive frequency and negative frequency freq = freq[0:int(len(freq)/2)] # Get positive frequency plt.figure(figsize=(12,3)) plt.plot(freq, am) plt.xlim(0,500) plt.vlines(x=156.13, ymin=0, ymax=0.2, colors='r') # Frequency doubling plt.vlines(x=156.13*2, ymin=0, ymax=0.2, colors='r') # Double frequency
Output results:
Red is the theoretical inner circle fault characteristic frequency, and the blue line is the actual fault characteristic frequency. Although it is not completely coincident, this is within the allowable range. Therefore, in the actual situation, if the envelope spectrum appears, the inner ring fault of the bearing can occur.
6. Compared with fft
Can we see the fault characteristic frequency without Hilbert transform demodulation and direct fft? The following is a comparative analysis
sampling_rate = 12000 am = np.fft.fft(xt) # fft is used to obtain the amplitude of at after Hilbert transform am = np.abs(am) # Calculate the absolute value of the amplitude (at this time, the absolute value is very large) am = am/len(am)*2 am = am[0: int(len(am)/2)] freq = np.fft.fftfreq(len(xt), d=1 / sampling_rate) # Obtain fft frequency, which includes positive frequency and negative frequency freq = freq[0:int(len(freq)/2)] # Get positive frequency plt.plot(freq, am)
Output results:
It can be seen that there are more fft high frequencies of the original inner ring fault signal
plt.plot(freq, am) plt.xlim(0, 500) plt.vlines(x=156.13, ymin=0, ymax=0.05, colors='r') # Frequency doubling plt.vlines(x=156.13*2, ymin=0, ymax=0.05, colors='r') # Double frequency
It can be seen that the fault characteristic frequency amplitude is low and is covered by other frequency amplitudes in the case of direct fft. In turn, Hilbert Demodulation can make it easier to observe the low frequency of fault characteristic frequency.
7. Envelope spectral function
In order to use envelope spectrum more conveniently, an envelope spectrum function PLT is encapsulated here_ envelope_ spectrum
def plt_envelope_spectrum(data, fs, xlim=None, vline= None):
'''
fun: Draw envelope spectrum
param data: Input data, 1D array
param fs: sampling frequency
param xlim: Abscissa of picture xlim，default = None
param vline: Picture vertical line, default = None
'''
#----To DC component----#
data = data - np.mean(data)
#----Do Hilbert transform----#
xt = data
ht = fftpack.hilbert(xt)
at = np.sqrt(xt**2+ht**2) # Obtain the analytic signal at = sqrt(xt^2 + ht^2)
am = np.fft.fft(at) # Obtain the amplitude by fft transformation of analytical signal at
am = np.abs(am) # Calculate the absolute value of the amplitude (at this time, the absolute value is very large)
am = am/len(am)*2
am = am[0: int(len(am)/2)] # Take positive frequency amplitude
freq = np.fft.fftfreq(len(at), d=1 / fs) # Obtain fft frequency, which includes positive frequency and negative frequency
freq = freq[0:int(len(freq)/2)] # Get positive frequency
plt.plot(freq, am)
if vline: # Whether to draw a vertical line
plt.vlines(x=vline, ymax=0.2, ymin=0, colors='r') # Height y 0-0.2, color red
if xlim: # Whether xlim is set for the abscissa of the picture
plt.xlim(0, xlim)
plt.xlabel('freq(Hz)') # Abscissa label
plt.ylabel('amp(m/s2)') # Ordinate label
7.1 outer ring fault data test
file_path = r'E:/graduate student/pytorch/CSDN code/fault_diagnosis_signal_processing/Part 4-Envelope spectrum/1730_12k_0.007-OuterRace3.mat' data = data_acquision(file_path) plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bpfo)
Output results:
It can be seen that the actual outer ring fault characteristic frequency is relatively obvious
7.2 test and analysis of rolling element fault data
file_path = r'E:/graduate student/pytorch/CSDN code/fault_diagnosis_signal_processing/Part 4-Envelope spectrum/1730_12k_0.014-Ball.mat' data = data_acquision(file_path) plt_envelope_spectrum(data = data, fs=12000, xlim=300, vline=bsf)
It can be seen that the actual rolling element fault characteristic frequency is not obvious
8. Summary
(1) Envelope spectrum can detect inner ring and outer ring faults, which is difficult for rolling elements
(2) Using fft directly is difficult to detect the fault characteristic frequency, and the fault characteristic frequency is easy to be covered up by high frequency