滚动轴承故障诊断(附MATLAB程序)

1、第二组实验轴承故障数据：数据打开后应采用 X105_DE_time作为分析数据，其他可作为参考， 转速 1797rpm轴承型号： 6205-2RS JEM SKF深, 沟球轴承采样频率： 12k Hz1、确定轴承各项参数并计算各部件的故障特征频率 通过以上原始数据可知次轴承的参数为： 轴承转速 r=1797r/min ；滚珠个数 n=9；滚动体直径 d=； 轴承节径 D=39mm；：滚动体接触角 =0 由以上数据计算滚动轴承不同部件故障的特征频率为：外圈故障频率 f1=r/60 * 1/2 * n(1-d/D *cos )=内圈故障频率 f2=r/60 * 1/2 * n(1+d/D *cos

2、 )= 滚动体故障频率 f3=r/60*1/2*D/d*1-(d/D)2* cos2( )= 保持架外圈故障频率 f4=r/60 * 1/2 * (1-d/D *cos )=对轴承故障数据进行时域波形分析将轴承数据 导入 MATLAB中直接做 FFT分析得到时域图如下：并求得时域信号的各项特征：（ 1）有效值：； （ 2）峰值：；（ 3）峰值因子：；（ 4）峭度：；（ 5）脉冲因子：；（ 6）裕度因子：3.包络谱分析对信号做 EMD模态分解，分解得到的每一个 IMF 信号分别和原信号 做相关分析，找出相关系数较大的 IMF 分量并对此 IMF 分量进行Hilbert 变换Empirical M

3、ode Decompositionim2fmi3fmi4fmi5fmi6fmi7fmi8fm由图中可以看出经过 EMD分解后得到的个 IMF 分量和一个残余量IMF 分量分别和原信号做相关分析后得出相关系数如下：EMD 分量IMF1IMF2IMF3IMF4IMF5相关系数EMD 分量IMF6IMF7IMF8IMF9IMF10相关系数由上表得： IMF1 的相关系数明显最大，所以选用 IMF1 做 Hilbert 包 络谱分析。所得 Hilbert 包络谱图如下：对包络谱图中幅值较大区域局部放大得到下图由以上包络图的局部放大图中可以看出包络图中前三个峰值最大也 最明显，三个峰值频率由小到大排列分

4、别为、。把这三个频率数值 和前文计算所得的理论值进行比较可知： 频率值最大为和内圈的故障 理论计算特征频率 f2=相近，说明此轴承的故障发生在轴承的内圈。 clc程序 1：原始信号时域分析及小波去噪处理clear all 轴承诊断 ');x1=(1:4096);clear z;N=4096;fs=12000;n=0:N-1;t=n/fs;f=n*fs/N;figure(1);plot(t,x1);xlabel('t');ylabel('幅值');title(' 原信号时域图 ')%小波去噪thr,sorh,keepapp=ddencmp(

5、'den','wv',x1); xd=wdencmp('gbl',x1,'db3',2,thr,sorh,keepapp);figure(2);plot(t,xd);xlabel('t');ylabel('幅值');title(' 小波去噪后时域图 ')程序 2：EMD 分解及 Hilbert 包络clcclear all轴承诊断 ');x=(1:1024);N=1024;fs=12000;n=0:N-1;f=n*fs/N;lag=N;n=0:N-1;t=n/fs;imf=e

6、md(x);m,n=size(imf); %imf 为一 m*n 阶矩阵， m 是 imf 分量， n 为数据点 emd_visu(x,1:length(x),imf,m); %实信号的信号重构及 emd 结果显示函数 for i=1:ma(i)=kurtosis(imf(i,:);%峭度b(i)=mean(imf(i,:); %均值 ;c(i)=var(imf(i,:); %方差 ;d(i)=std(imf(i,:); %均方值 e(i)=std(imf(i,:).; %均方根值f(i)=skewness(imf(i,:); %计算偏度endk,c=max(a); %k为峭度最大值， c 为

7、最大元素在数组中的位置 r,lags=xcorr(x,lag,'unbiased'); %计算序列的自相关函数 for i=1:mR,lags=xcorr(imf(i,:),lag,'unbiased'); %计算序列的自相关函数 a=corrcoef(R(1:N/2),r(1:N/2); %相关系数矩阵【对称】，主对角元素为 1 xg(i)=abs(a(1,2); %相关系数endR,C=max(xg); %R为最大值， C 为最大元素在数组中的位置figure(4);y = hilbert(imf(C,:);a = abs(y);%包络 b=fft(a);mag1=abs(b); mag=mag1*2/N; f1=(0:N-1)*fs/N;plot(f1(1: