基于接触模型和峭度最优Laplace小波的滚动轴承量化诊断✳

袁 幸,朱永生,张优云,洪 军

(1.西安交通大学润滑理论及轴承研究所,陕西 西安 710049;2.西安交通大学机械制造系统工程国家重点实验室,陕西西安 710049)

引 言

飞机、电力机车、高档机床等重大装备的核心技术问题之一就是轴承,它对整个装备制造业的发展水平有着举足轻重的作用,其运行状态往往直接影响整台机器的精度、可靠性及寿命[1]。 Sawalhi N将最小熵盲反卷积和谱峭度应用于滚动轴承故障检测中[2]。Karthik K应用互信息进行了故障特征选取[3]。谭继勇针对滚动轴承冲击信号的检测问题,构造了随机共振自适应检测算法,提高了信噪比,避免了漏峰现象[4]。文献[5]提出了基于经验模态分解和奇异值差分谱的轴承故障诊断方法。从已有研究结果来看,以先进信号分析为基础的滚动轴承诊断技术研究较多,人们大都针对损伤信号处理与分类,研究相关的算法,却很少关注服役轴承具体的损伤尺寸。对于滚动轴承这样的基础部件,应用广泛但规格种类繁多,工况差别大,不同损伤程度的样本获取困难。以机车为例,轮对轴承的故障样本少,一个轮对检修流水线在3个月时间检测的1 200多个轮对轴承中仅有2例故障轴承[6]。再者,单纯的基于信号处理的传统方法较少考虑轴承系统的具体结构特点,侧重于对故障存在性的检测,不易描述损伤演化过程,因而难以准确评价轴承可以持续运行的安全余量。

本文系统地提出根据正、反问题进行损伤量化诊断的原理及方法。建立多体接触振动模型来描述滚动轴承空间运动特性,涵盖轴承界面复杂相互作用及摩擦学系统行为,以接触表面局部区域尺寸变化来刻画损伤发展趋势,实现复杂工况下轴承的损伤建模。通过精细划分损伤区域由正问题求解标准模式数据库,将反问题计算转化为多维优化问题,运用灰色关联度辨识轴承损伤位置和具体尺寸。

1 滚动轴承损伤接触建模

1.1 接触刚度系数

由 Hertzian理论,点接触负荷-变形关系为[7]

式中 W为弹性趋近量,Eschmann给出了两固体点接触弹性变形表达式

式中 Δ,kb,_是由材料属性决定的 Hertzian系数;d分别为弹性模量、泊松比与接触表面曲率和。点接触的接触力可表示为

因此轴承接触刚度系数为

∑ d的计算方法由 Harris给出,详见文献 [7],一旦轴承设计尺寸和材料确定,K可由式(2)～(4)求出。

1.2 接触关系的坐标表达

图1是滚动轴承坐标示意图,第i个钢球-套圈接触变形Wi是内圈在X,Y方向位移(xs,ys),钢球位置角θi和游隙c的函数,如下式

设(xb,yb)为钢球坐标,由于振动传递作用,考虑到钢球自身的振动,局部接触变形变为

θi表示第i个钢球的位置角

式中 kc为轴承公转速度,N为钢球个数。设转速为k,则,其中Db为钢球直径,Dp为轴承节圆直径。

图1 滚动轴承示意图Fig.1 Schematic diagram of rolling bearing

设(xo,yo)为外圈处的运动坐标,接触变形同样可表示为

1.3 动态接触建模

球轴承接触阻尼随平均负荷变化,其值很小,Kramer通过频域测试方法给出了阻尼近似计算公式[8]

故xb方向包含阻尼的动态接触力可写成

1.4 拖动力模型

点线接触中两个滚动元件之间的相对滑动速度引起的润滑剂与滚动元件之间的摩擦力被称为拖动力。拖动特性是决定轴承元件动态特性的关键因素之一,是滚动轴承动力学分析中不可缺少的基本参数。引入摩擦系数_,yb和xb方向拖动力的函数关系为

摩擦系数_由弹流润滑模型决定,设Λ为油膜厚度与表面粗糙度的比值,则[7]

式中 _bd受润滑参数及温度的影响,在本文使用的润滑模型中,摩擦系数_bd是滑 /滚率s的函数。

通过以上分析,对深沟球轴承的平面运动,考虑到刚度、阻尼、拖动力和滑 /滚率的动态接触力表达式为

1.5 运动方程式

滚动轴承的振动存在着振动源和振动传递路径的交互,如图2所示,系统响应最终为传感器接收。

图2 滚动轴承系统振动模型Fig.2 Vibration model of rolling bearing system

轴承表面损伤形态决定了激振力的频谱,界面动态接触决定了振动系统的传递特性,最终的振动频谱由上述二者共同决定。正问题分析中,外圈质量和刚度将激励系统较低的固有频率。考虑到轴承的高频谐振频带,设置mp,cp,kp为附加自由度来模拟轴承-轴承座结合部的约束以激励高频响应。将传递路径的质量、刚度、阻尼及几何位置作为系统参数,根据牛顿定律可建立内圈-钢球-外圈接触振动微分方程组如下

式中ms,mb,mo分别为内圈与轴的质量,钢球质量,外圈质量;ko,co为外圈刚度与阻尼;fx,fy为径向载荷与不平衡力的总和。由于工程环境、安装条件等因素的影响,使得轴承-轴承座结合部的材料特性、几何尺寸和边界条件具有微小的偏差。为了与实际工况等效,mp,cp,kp的具体数值可通过健康信号的频谱加以修正。模型参数导致的误差做了误差置零处理[9],获得数值计算与实验结果一致的数值模型。

1.6 表面损伤建模

滚动轴承由疲劳诱发的表面损伤主要包括裂纹、凹点、剥落等。最早研究轴承损伤振动模型的是McFadden和 Smith,该模型通过轴承运转过程中钢球经过故障位置产生的一系列冲击与冲击衰减函数的卷积来模拟振动[10]。Choudhury A提出外加冲击力的形式模拟损伤所激励的动态效应,建立了滚动轴承表面损伤的随机振动数学模型[11]。Behzad M等通过基于模型的振动响应与传递函数相结合的方法较好地模拟了轴承的振动行为[12]。

到目前为止,提出的大多数为滚动轴承损伤的数字信号仿真与模拟冲击模型,未能有效地将损伤融入到滚动轴承自身的几何学和动力学方程中。本文通过找出滚动轴承在运转过程中钢球通过损伤区的运动变化规律,在宏观上建模并与接触振动相结合为分析滚动轴承疲劳累积损伤提供一种可用的方法,损伤模型如图 3所示。

损伤长度为L,深度为h,当钢球运行至损伤区域时,游隙值将随之增加h,使得W=W+h,瞬间通过长度L后,游隙值恢复为W。当钢球处于损伤区域时,导致接触面间的 Hertzian接触力突然降低或变为零,相当于滚动轴承的非线性时变系统的激励发生突变使得振动能量出现波动,损伤程度可用L×h或L/h来定量衡量。

图3 表面损伤实物与几何模型Fig.3 Physical and geometric model of surface damage

2 数值计算与分析

式(15)所代表的状态方程组非线性强,难以解析,本文采用数值方法获取系统响应。在MATLAB中编制了相应程序,计算方法为定步长 Runge-Kutta法。算例采用美国凯斯西储大学(CWRU)轴承研究中心的深沟球轴承 (6205-2RS SKF),转速1 797 r/min,采样频率 12 000 Hz,CWRU提供了轴承详细设计参数[13]。轴承受到来自转子不平衡激励所产生的强迫振动,其振动频率为转子的旋转频率。轴承也将产生由于总刚度连续周期变化而形成的VC(Varying compliance)振动及其谐波组合[14,15]。这是因为,在径向载荷的作用下,各钢球的受力情况是不一样的,且随着钢球上的某一点的运动位置不同受力情况亦不一样,随着钢球相对于径向载荷作用线的移动,轴承刚度以数倍于钢球沿静止套圈转动的频率呈周期性变化。这种振动为参数激励的通过振动,VC振动频率与外圈损伤频率数值上相等属于轴承自身本体振动成分,VC=kc×nb=106.7 Hz。图4的功率谱清晰地显示了VC振动及其高次谐波,初步表明文中计算具有一定的可靠性。

图4 正常状态波形及功率谱图Fig.4 Waveform and spectrum of the normal bearing

图5～7分别是损伤信号波形及功率谱图,内、外圈和钢球损伤频率依次为fi=161.3 Hz,fo=106.7 Hz,fb=140.4 Hz。与图 4比较,损伤时波形具有明显的周期冲击特征,这是游隙瞬间变化引起的,波形重复出现的时间间隔为损伤出现的周期0.006 2,0.009 4,0.007 1 s。

由于接触振动模型全面考虑了振动的传递路径,诱发轴承各阶固有频率,使得高频谐振频率与损伤频率调制,在高频处谱图出现以损伤频率为主的边频带,幅值远离载波中心程逐渐减小的趋势。由于内圈损伤位置随轴旋转不断改变,损伤引起的冲击响应幅值必然受到转频调制,会在损伤频率各次谐波处出现边频谱线,造成内圈损伤功率谱图更加复杂,如图5所示。试验信号与仿真结果取得较好的一致性,也说明了本文建模方法的有效性。

图5 内圈损伤状态波形及功率谱图Fig.5 Waveform and spectrum of the bearing with inner race damage

图6 外圈损伤状态波形及功率谱图Fig.6 Waveform and spectrum of the bearing with outer racedamage

图7 钢球损伤状态波形及功率谱图Fig.7 Waveform and spectrum of the bearing with rolling element damage

3 故障量化特征提取

3.1 故障演变过程描述

将L×h设为轻度、中度、严重 3个阶段,分别取值为(0.1,0.1),(0.6,0.8),(1.2,1.5)mm。图8中(a)～(c)为 3个阶段的外圈损伤频谱,可见谱图结构差异明显。轻微损伤时损伤频率被谐波分量所干扰,低频调制边频能量较弱,随着损伤程度加深低频调制谱线逐渐清晰,并且谱图能量分布变化显著。图8(d)为损伤发展时归一化振动烈度变化曲线,损伤加剧冲击越明显,振动烈度值越大;冲击越不明显,甚至被VC及其谐波频率所淹没时,振动烈度越小。不难看出,通过剖分L×h的数值可以有效区分振动响应。以L×h的大小量化损伤区域等效于局部不同柔度的植入,从而诱发振动响应变化,形成模式的多样性。

图8 外圈不同损伤程度的频谱及振动烈度值Fig.8 Spectrum and vibration intensity of the bearing with outer race damagefor different fault severities

3.2 峭度最优Laplace小波特征提取

机械诊断属于模式识别问题,为了增强特征向量的敏感性,达到精确提取轴承损伤时冲击响应信号的目的,本文研究了一种峭度最优 Laplace小波。R Lind等人构造的 Laplace小波是一种复指数小波,解析表达式为[16]

式中 参数矢量γ=[a,kn]决定了小波的性能;成员变量a,kn为模态动力学参数;A用来归一化小波函数。峭度对冲击信息非常敏感,又是区分非高斯分布的指标,高的峭度值意味着信号含有丰富的冲击成分。基于波形匹配的思想,本文用峭度最大化准则定量描述参数矢量γ。在一定范围内变化a和kn,选择使得峭度最大的a,kn值为参数,经计算获得最优值a=0.85,kn=16,这就是峭度最优 Laplace小波变换。本质上,基于峭度最优的 Laplace小波变换也就是控制小波滤波器的通带带宽,参数优化相当于调整品质因子Q,最终结果使小波冲击特征更加显著,得到与信号最相似的基元函数。

时域峭度、脉冲等 5个指标已被证明可较好区分多重状态[17]。仅依靠时域信息对轴承损伤尺寸的精确识别还不够理想,因为信息含量太少,频谱结构变化得不到体现。因此,结合频域无量纲指标组建标准模式,频域指标定义为平均频率、频谱偏斜度、频谱峭度、均方根比[18]。由于工况变动的随机性和采集信息的噪声影响,量化诊断具有一定的不确定性。模糊分析是解决识别过程中不确定性问题的有力工具。以 [0,1]之间的模糊隶属度表示特征向量有助于提升特征敏感度,优化特征空间布局,增强信息可靠性。根据经验可选择升半岭型函数表示隶属度[17]。至此,构造了量化诊断标准模式向量

以外圈损伤尺寸 0.001,0.002,0.003 mm为例,根据模型计算特征向量,结果如图9所示。可见,经峭度最优 Laplace小波提取的 9个无量纲指标对差值为 0.001 mm的损伤尺寸实现了良好的分类,也说明动力学响应参数可以较好映射局部损伤尺寸,为量化诊断奠定了基础。

图9 外圈不同损伤尺寸的特征向量Fig.9 Feature vectors of the bearing with outer race damage for different fault size

4 局部损伤的定量检测

4.1 定量识别的灰色关联度分析

量化诊断正问题实际上是通过对含有不同损伤位置和尺寸的滚动轴承进行动力学求解,获得损伤特征并转化为模糊向量,由计算得到的离散值构建损伤数据库。本文将反问题转化为轴承的实际振动输出与接触振动数值模型的数值输出间差异最小化为优化目标的优化问题。利用关联寻优求解出与输入值最相似的样本点,进而测出损伤位置及损伤长度L,深度h。正问题计算得到标准模式数据库

式中i为模式类别,k=9表示特征向量个数。令待检模式为ςj(k)={xj(1),xj(2),… ,xj(9)},其中j为模式类别,k的意义同上。记Si对Sj的关联系数为Yij(k),Yij(k)由下式确定[19]

式中 Δij=|Si(k)-Sj(k)|称为绝对差,d为分辨系数,量化诊断时可取d=0.5,则Si与Sj的关联度为

关联度从物理意义讲是对灰色过程的白化处理,从数学角度讲是空间几何曲线的比较。反问题辨识损伤尺寸可用如下数学关系式描述

正、反问题进行滚动轴承量化诊断的流程如图10。

图10 量化诊断流程Fig.10 Flow chart for the quantitative identification

4.2 案例验证

案例同样源自 CWRU轴承中心,采用电火花机床分别在测试轴承的内、外圈和钢球上加工深度为11 mils(0.279 4 mm)的盲孔,盲孔的直径为 7 mils(0.177 8 mm),14 mils(0.355 6 mm),21 mils(0.533 4 mm),这样就可以通过实验数据完成量化诊断。测试轴承分别在 4种不同载荷(载荷分别为0～ 3马力)以及相对应的 1 797,1 772,1 750,1 730 r/min四种转速下测试。振动信号由16通道DAT数据记录仪采集,采样频率为12 k Hz。 CWRU未提供第4种工况下的外圈数据,因此本文用内、外圈和钢球在前 3种工况下的样本各 20个按照图 10的流程进行量化诊断。

正问题计算区间定为L=0.1～1 mm,h=0.1～1 mm,搜索步长0.001 mm。表1给出了任意一种工况下h的反问题分析结果,表 2为 3种工况下L的分析结果。计算出损伤尺寸预测值的相对误差,结果表明,损伤深度预测相对误差不超过 4%,损伤长度的相对误差不超过 7%。

表1 量化诊断结果(h/mm)Tab.1 Results for the quantitative identification(h/mm)

表 2 量化诊断结果(L/mm)Tab.2 Results f or thequantitative identif ication(L/mm)

5 结 论

(1)滚动轴承的振动响应是轴承内部几何关系、运动关系在外部载荷作用下的综合体现。结合轴承界面复合约束作用及摩擦学行为,建立了计及传递路径的滚动轴承动态接触振动模型,获取了轴承不同损伤状态下的本体振动信号。该模型可以敏锐地捕捉到轴承损伤时振动信号的冲击信息且涵盖了频谱的中、高频成分,理论结果与案例数据分析相吻合。

(2)局部损伤的出现引起内部游隙的变化,进而影响到轴承的振动特性。不同损伤位置和尺寸形成了不同的损伤模式;反之,运用模式向量也可反求损伤位置和具体尺寸。

(3)峭度对信号中的冲击成分非常敏感,同时峭度又是一种无量纲指标,它不受信号绝对水平的影响。根据特征波形匹配原理,研究峭度最优 Laplace小波复合特征提取方法,全面挖掘了振动信号时、频域动态信息。

(4)将量化诊断中的反问题转化为优化辨识,以实测特征向量作为输入,运用灰色关联度辨识具体损伤尺寸。诊断实践表明,识别误差不超过7%,这为方法的进一步实际应用提供了有力的依据。

[1] Robert B Randall,Jére Antoni.Rolling element bearing diagnostics—— A tutorial[J].Mechanical Systems and Signal Processing,2011,25(2):485—520.

[2] Sawalhi N,Randall R B,Endo H.The enhancement of fault detection and diagnosis in rolling element bearings using minimum entropy deconvolution combined with spectral kurtosis[J].Mechanical Systems and Signal Processing,2007,21(6):2 616—2 633.

[3] Karthik K,Nataraj C.Feature selection for fault detection in rolling element bearings using mutual information[J].Journal of Vibration and Acoustics,2011,133(6):061001—11.

[4] 谭继勇,陈雪峰,何正嘉.冲击信号的随机共振自适应检测方法 [J].机械工程学报,2010,46(23):61— 67.Tan Jiyong,Chen Xuefeng,He Zhengjia.Impact signal detection method with adaptive stochastic resonance[J].Journal of Mechanical Engineering,2010,46(23):61—67.

[5] 张超,陈建军,徐亚兰.基于 EMD分解和奇异值差分谱理论的轴承故障诊断方法[J].振动工程学报,2011,24(5):539—545.Zhang Chao,Chen Jianjun,Xu Yalan.A bearing fault diagnosis method based on EMD and difference spectrum theory of singular value[J].Journal of Vibration Engineering,2011,24(5):539—545.

[6] 屈梁生.机械监测诊断中的理论与方法 [M].西安:西安交通大学出版社,2009.Qu Liangsheng.Theory and Methodology for Mechanical Condition Monitoring and Fault Diagnosis[M].Xi'an:Xi'an Jiao Tong University Press,2009.

[7] Harris T A.Rolling Bearing Analysis[M].4th ed.New York:John Wiley and Sons,2001.

[8] Kramer E.Dynamic of Rotors and Foundations[M].Berlin:Springer,1993.

[9] 陈雪峰.小波有限元理论与裂纹故障诊断的研究[D].西安:西安交通大学,2004.Chen Xuefeng.Research on the Theory of Wavelet Finite Element Method and Cracks Fault Diagnosis[D].Xi'an:Xi'an Jiao Tong University,2004.

[10]Mc Fadden P D,Smith J D.Model for the vibration produced by a single point defect in a rolling element bearing[J].Journal of Sound and Vibration,1984,96(1),69—82.

[11]Choudhury A, Tandon N.Vibration response of rolling element bearings in a rotor bearing system to a local defect under radial load[J].ASME,Journal of Tribology,2006,128(2):252—261.

[12]Behzad M,Bastami A R,Mba D.A new model for estimating vibrations generated in the defective rolling element bearings[J].ASME,Journal of Vibration and Acoustics,2011,133(4):041011—8.

[13]The Case Western Reserve University Bearing Data Center Website.Bearing data center seeded fault test data[EB/OL].http:http://www.eecs/cwru/edu/laboratory/bearing.

[14]TIWARI M,GUPTA K,PRAKASH O.Effect of radial internal clearance of a ball bearing on the dynamics of a balanced horizontal rotor[J].Journal of Sound and Vibration,2000,238(5):723—756.

[15]TIWARIM,GUPTA K,PRAKASH O.Dynamic response of an unbalanced rotor supported on ball bearings[J].Journal of Sound and Vibration,2000,238(5):757—779.

[16]何正嘉,訾艳阳,张西宁.现代信号处理及工程应用[M].西安:西安交通大学出版社 ,2007.He Zhengjia,Zi Yanyang,Zhang Xining.Modern Signal Processing and Engineering Application[M].Xi'an:Xi'an Jiao Tong University Press,2007.

[17]徐敏,虞和济,张瑞林.设备故障诊断手册 [M].西安:西安交通大学出版社,1998.Xu Min,Yu Heji,Zhang Ruilin.Handbook of Machinery Fault Diagnosis[M].Xi'an:Xi'an Jiao Tong University Press,1998.

[18]Chen P,Masatoshi T,Toshio T.Fault diagnosis method for machinery in unsteady operating condition by instantaneous power spectrum and genetic programming[J].Mechanical Systems and Signal Processing,2005,19(1):175—194.

[19]邓聚龙.灰色系统基本方法 [M].武汉:华中科技大学出版社,2005.Deng Julong.The Primary Methods of Grey System Theory[M].Wuhan: Huazhong University of Science& Technology Press,2005.