滚动轴承复合故障振动建模及Lempel-Ziv 复杂度评价

朱永生，袁 幸，张优云，朱川峰，周 智

(西安交通大学 润滑理论及轴承研究所，西安 710049)

滚动轴承被广泛应用于机械工业，是基础件也是关键核心件。由于其通常工作在高速、高温、重载的恶劣环境下，滚动轴承疲劳导致的失效往往在整机故障中占大部分比率。因此故障诊断对于滚动轴承的安全、连续运行和预知更新至关重要［1］。故障建模研究故障状态下动力学参数与响应征兆的内在联系，籍此获取服役正常与异常的证据信息［2］，其目的在于根据机械设备性能退化响应特点来识别潜在故障和构建故障库。

滚动轴承动力学建模一直为学者所重视，Fukata［3］较早分析了滚动轴承的变柔性振动。研究结论显示，滚动轴承支承下的不平衡转子具有两个激励源，即来自不平衡的转频激励和来自轴承总体刚度周期变化的内部参数激励。Tiwari等［4］采用两自由度模型研究轴承游隙造成的转子系统的非线性行为和稳定性。然后Lioulios等［5］探讨了速度波动对轴承非线性动力学的影响。Patil等［6］模拟内、外圈局部故障，研究单点缺陷所激励的时、频域特性。Liew等［7］进一步研究非线性Hertzian接触、游隙、转速等因素对动态行为的影响机制，并有效模拟了轴承元件的振动形态。Sopanen等［8－9］提出一种轴承模型用于分析局部及分布式故障，模型考虑了非线性Hertzian接触，但是缺乏外圈－轴承座的耦合效应。澳大利亚UNSW大学Sawalhi等［10］采用附加自由度的方法建立5DOF模型，从而得到了相对精确的振动响应。一些新的建模手段也被引入，Mohsen等［11］用键合图功率流构建滚动轴承振动模型，分析了轴承元件的动态接触行为。机械设备出现故障时，故障的演变存在一个由轻微到严重、由单一到复合的过程［12］。对于长期在变载工况下服役的机床轴承、机车轴承尤为显著，监测到的故障往往以多故障耦合的形态出现，研究滚动轴承复合故障建模有着重要工程意义［13］。然而，现有的动力学模型均以单点缺陷为研究对象，复合故障建模的研究鲜见报道。

为能够有效预示多耦合故障诱发的振动特征，基于轴承运动学与Hertzian接触理论，建立一种计及内圈－外圈－轴承座的6DOF振动模型。通过渐进过程描述缺陷区域，以深沟球轴承SKF6205－2RS为实例计算振动响应。与实验数据对比表明，利用该方法可深入揭示轴承局部缺陷动力特性与宏观动态性能之间的映射关系。进一步着重研究复合故障建模方法及其响应特点。在此基础上引入非线性动力系统中复杂性刻画方法——Lempel-Ziv复杂度指标用于评价不同故障状态动力系统的变化。研究表明，这种特征量是表征不同故障信息的有效参数。

1 滚动轴承接触建模

1.1 非线性弹性恢复力

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

如图1，第i个滚动体－套圈弹性趋近量δi是套圈的径向位移(xs，ys)，滚动体位置角θi和游隙γ的函数，由运动关系可知:滚动体位置角ωct，其中:ωc为公转速度，n为滚动体个数。设ω为相配轴的转速，Db为滚动体直径，Dp为轴承节圆直径，ωc设轴承座运动坐标为(xp，yp)，则内圈

－轴承座的综合弹性趋近量就为:需要说明的是，由于载荷分布区具有角度范围，δi只能取正值。

1.2 载荷 －位移关系［14］

根据Hertzian理论，点接触载荷－位移关系为:

式中:δ为弹性趋近量，K为Hertzian接触刚度。为了描述两个啮合表面之间的接触状况，定义Kinner、Kouter、∑ρ分别为滚动体与内、外圈的接触刚度与曲率和，则:

其中:E为弹性模量，ν为泊松比，μ*为变形系数。这样，一旦轴承几何尺寸确定，Kinner、Kouter就容易得到。进而求出轴承整体等效接触刚度为:

滚动轴承运动副的等效约束力可由每个滚动体所受弹性接触力求得，故总的非线性弹性恢复力可表示为:

1.3 接触阻尼

滚动轴承接触阻尼随平均负荷变化，其值很小，Kramer通过动态测试方法给出了经验计算公式［14］

2 振动系统模型

滚动轴承振动性能是其内部运动关系在外部载荷作用下与传递路径的综合作用。内圈－外圈－轴承座为非线性强耦合关系，图2表达了轴承接触界面的多自由度振动形态。其中，c0、k0、cp、kp分别为外圈、轴承座的阻尼和刚度。基于多自由度振动理论，构建多体系统环境下的轴承精细化模型如图2所示。由于滚动轴承接触问题的复杂性及各种非线性因素的影响，无疑会增加系统建模及求解的难度。鉴于此，本研究做以下基本假设。① 忽略滚动体滑移现象。② 忽略制造误差，如表面波纹度的影响。③ 假设滚动体与内、外圈之间只存在局部弹性接触变形。由第一节的分析结合图2～3，建立振动方程如下:

图2 滚动轴承振动示意图Fig.2 Schematic diagram of rolling bearing vibration

式中:ms为内圈(相配轴)质量，mo、mp分别为外圈和轴承座质量，e为不平衡距离，fn为外载荷。

3 局部故障嵌入

滚动轴承的失效大部分表现为局部缺陷，例如疲劳剥落、凹痕、点蚀等。缺陷的存在会形成一个有限的区域，具有长度l，如图3所示。当滚动体运行至缺陷位置时，弹性趋近量瞬间增加，设d为相应的实际增加值，即式(2)中的游隙γ增加d，这将导致冲击产生和能量突变。在外圈固定内圈旋转的情况下，外圈缺陷所对应的中心角α=2arcsin(l/2R)，R为外圈半径。设缺陷区域相对初始滚动体的位置角为，则滚动体通过外圈缺陷区域的渐进过程可描述为:

图3 表面缺陷区域Fig.3 The surface defects region

滚动体出现故障时，缺陷区域以转速ωb自转一周与内、外圈分别接触一次，每一次接触弹性趋近量均增加d，此时判别函数由两部分组成。设滚动体缺陷区域与外圈初始相对角位置为，则缺陷区域通过外圈的渐进过程为:

同理，滚动体缺陷区域通过内圈的渐进过程为:

复合故障指同时发生两个或多个互相关联、交叉影响的故障。以内、外圈复合故障为例，动态建模表现为式(10)～(11)的级联，即同时运行内、外圈的渐进过程。

4 计算分析及验证

式(9)所代表的状态方程组非线性强，难以解析，本文采用数值方法获取系统响应。计算方法为定步长Runge–Kutta法，步长设为1/12 000 s，初始值为10－6m。算例采用美国凯斯西储大学(CWRU)轴承研究中心的深沟球轴承(SKF6205－2RS)，转速1 772 r/min，采样频率12 000 Hz。故障轴承采用电火花加工制成，内、外圈及滚动体的故障频率依次为5.415 2、3.584 8、4.713 5乘以工频，CWRU提供了详细设计参数可供参考［15］。外圈发生故障时频谱结构相对简单，主频为故障特征频率，如图4所示。可以看出，仿真信号包络谱体现了不平衡激励的转频fr=29.66 Hz，内圈故障频率fo=106 Hz及其谐波2fo、3fo、4fo成分，这与实验信号非常吻合。与外圈故障的区别在于内圈缺陷连同转轴一起运行，故障频率必然受到工频的调制，因此包络谱的主频以 fi=159.9 Hz、fi+fr=189.5 Hz、fi－ fr=130.4 Hz及其组合谐波成分 2fi+fr、2fi－ fr、3fi+fr、3fr－ fr等为主。图5的仿真信号包络谱成功突出了主频成分，并且频率预测达到了很高的精度。仿真和实验数据的对比表明，振动模型对显著频率成分实现了良好的预测。为验证模型的鲁棒性，考察转速为1 797 r/min时滚动体故障的振动响应。从图6中清晰看出滚动体缺陷激励引起的振动在包络谱中表现为保持架转频fc=11.9 Hz对故障频率fb=141.2 Hz的调制。这是由于滚动体上缺陷位置随保持架旋转不断改变，缺陷引起的冲击响应幅值必然受到保持架旋转频率调制。以上分析充分体现了本文建模方法的可靠性。

图4 外圈故障仿真与实验信号包络谱Fig.4 Simulation and experimental envelope spectrum of outer race defect

图5 内圈故障仿真与实验信号包络谱Fig.5 Simulation and experimental envelope spectrum of inner race defect

设置相同的缺陷长度l=0.5 mm和深度d=0.1 mm，转速1 772 r/min时内、外圈复合故障振动响应如图7所示。滚动轴承复合故障的信号特征比单一故障相对复杂。从刚度损失的角度可以解释为复合故障导致轴承接触刚度发生更大变化，在不同缺陷位置引入局部柔度，使得动力学性能改变。当内、外圈复合故障出现时，局部柔度同时增加，激励系统出现两种故障频率fo和fi。包络谱低频部分fo及其谐波比较明显，fi的幅值相对较小。可以解释为复合故障由单点故障呈矢量叠加方式构成，以内、外圈复合故障为例，频域总幅值为:

图6 滚动体故障仿真与实验信号包络谱Fig.6 Simulation and experimental envelope spectrum of rolling element defect

图7 内、外圈复合故障信号及包络谱Fig.7 Vibration response of inner and outer race compound defects

图8 外圈、滚动体复合故障信号及包络谱Fig.8 Vibration response of outer race and rolling element compound defects

图9 内圈、滚动体复合故障信号及包络谱Fig.9 Vibration response of inner race and rolling element compound defects

图10为内圈、外圈、滚动体三种故障复合时的振动响应，包络谱清晰显示了三种特征频率的2倍频，fb=139.2 Hz相对微弱，并且各频率的幅值均不同。计算结果充分表明滚动轴承复合故障具有相互耦合、相互制约和强度不平衡性，也说明这些故障的产生具有内在确定性机制。同样显示，包络谱所提取的故障频率倍频分量可作为识别复合故障的有效依据。

5 复合故障Lempel-Ziv复杂度评价

滚动轴承复合故障往往使其频率组成及各成分幅、频特性复杂化，从而使时域波形发生不同程度的畸变，这种畸变也是轴承故障在时域信号中的直接反映。事实上，信号的波形、频谱结构发生变化时，信号的复杂程度也发生了变化。人们开始探索采用一些指标定量地描述轴承振动信号的复杂性［2］。姜建东等［2］用Lempel-Ziv复杂性测度定量评估大机组运行状态。Hong等［16］研究了Lempel-Ziv复杂度指标用于评估故障轴承损伤状态。实践表明，这一指标是衡量有限时间序列复杂程度的高效工具。

图10 内圈、外圈、滚动体复合故障信号及包络谱Fig.10 Vibration response of inner race，out race and rolling element compound defects

序列 SN={s1，s2，…，sN}的 Lempel-Ziv 复杂度可按CN(r)(r≤N)经过 N 次循环得到，计算过程如下［2、16］:

(1)初始化 Sv，0={}，Q0={}，CN(0)=0，r=1。令 Qr={Qr－1Sr}，由于 Qr不属于 Sv，r－1，则 CN(r)=CN(r－1)+1，Qr={}，r=r+1。

(2)令 Qr={Qr－1Sr}，判断 Qr是否属于 Sv，r－1=Sv，r－1SrSr－1;若属于，则 CN(r)=CN(r－1)，r=r+1，重复过程(2)。

(3)若不属于，则 CN(r)=CN(r－1)+1，Qr={}，r=r+1，重复过程(2)。

例如序列{0000…}的分段形式为{0·000…}，CN=2。符号序列{010101…}的分段形式为{0·1·0101…}，CN就为3。出于对数据严格的可比性考虑，Lempel和Ziv进一步提出了一种归一化基准公式，将CN值界定在［0 1］之间，详细算法为:

其中:k为SN中元素的个数，对二进制序列k=2。根据已有的经验，一般N大于3 600可保证上式成立［16］。根据上述Lempel-Ziv复杂度算法和文献［2，16］可知，一个序列的复杂度越大，说明它的添加操作越多，序列的内涵分量也越多，描述给定符号序列所需最少的、互不相同的分段形式越多，给定的序列周期性越弱，出现新内涵模式的速度也就越快。Lempel-Ziv复杂度的物理意义在于它反映了一个时间序列随着长度的增加出现新模式的速率。复杂度越大，说明数据在窗口长度时期内随时间出现的新变化越多，发生新变化的速度越快，表明这一时期的数据变化是无序而复杂的。反之复杂度越小，则说明发生新变化的速度越慢，数据变化是规则的，周期性越强。因此，振动信号的Lempel-Ziv复杂度指标能够客观反映出系统状态发生变化的情况。

为进一步完善轴承故障评估方法，本文借助上述Lempel-Ziv复杂度指标的建立思想，并将之推广到复合故障诊断领域。同样设置缺陷长度l=0.5 mm和深度d=0.1 mm的微弱故障进行数值计算。图11从左往右依次对应正常、外圈、内圈、滚动体、内－外圈复合、外圈－滚动体复合、内圈－滚动体复合、内外圈－滚动体复合故障的Lempel-Ziv复杂度变化趋势，表1为其具体数值。

图11 不同状态的Lempel-Ziv数值变化趋势Fig.11 Lempel– Ziv complexity vs.different fault

表1 不同状态的Lempel-Ziv复杂度Tab.1 Numerical value of Lempel-Ziv complexity

从结果中可以看出，正常状态的振动最为简单。滚动体故障虽然激励冲击能量不大，但是滚动轴承中滚动体的运动比较复杂，既有公转又有自转。在其运动过程中，故障区域既要与内圈接触又要与外圈接触，渐进过程由两部分组成(式(12)～(13))。不断变化的接触行为会引起强烈的调制作用，反映在频谱中会出现越来越多的谱峰及其谐波，波形中的成分越来越杂乱，振动的周期性减弱，非线性增强，相对于内、外圈更加倾向于较高的Lempel-Ziv值。复合故障破坏了单一故障的周期冲击振动结构，故障冲击加剧，直观上可以看出，信号的复杂度增加了。不同形式的复合故障其响应的复杂度也不同，这是因为不同的复合故障形成机制不同，导致信号的冲击振动结构有所差异。内外圈－滚动体三种缺陷复合故障的振动响应与两种缺陷复合故障的Lempel-Ziv指标差别不太大，体现了故障之间既相互干扰又相互制约。对复合故障的复杂度评价与我们的工程经验相吻合，再次证明了本文故障建模方法的可靠性。由上述分析可知，振动信号的复杂度较为敏感地反映了轴承状态演变过程中动力系统的变化历程，是一种有效的特征量。实际应用中整个信号的Lempel-Ziv指标单独作为诊断特征量稍微少了一些。基于波形匹配思想，兼备优良性质的多小波拥有多个尺度函数和小波函数，可以很好地匹配复合故障的多种特征波形［13］。利用多小波变换消除噪声和分解信号，从而构造不同尺度下的Lempel-Ziv复杂度向量，可以同时完成单点及复合故障的监测诊断。此方法的实验研究加之工程应用，亦是下一步的研究课题，将另文给出。

6 结论

(1)建立一种计及内圈－外圈－轴承座的滚动轴承6DOF振动模型。从轴承几何学、运动学的角度推导了局部故障的渐进过程，通过与实验数据的比较表明，该方法能有效预测振动频谱的主要成分，为后续复合故障的研究奠定了基础。

(2)进一步研究了复合故障的建模方法和振动形态。研究结果显示，复合故障具有相互耦合性、相互制约性和强度不平衡性。包络谱或共振解调谱可同时解耦复合故障特征频率及其倍频分量。一定条件下，包络谱的倍频分量可以作为诊断复合故障的依据。

(3)滚动轴承局部故障的出现相当于局部柔度的引入，进而在不同渐进过程瞬间改变系统动力学特性。不同故障信号的Lempel-Ziv复杂度指标反映了不同故障的动力学形成机制。轴承振动信号的复杂度能有效地评价轴承运行状态，灵敏地反映其状态的变迁过程。

［1］RandallR B， Antoni J. Rolling element bearing diagnostics—A tutorial［J］.Mechanical Systems and Signal Processing，2011，25(2):485－520.

［2］机械监测诊断中的理论与方法:屈梁生论文集［M］.西安:西安交通大学出版社，2009.

［3］Fukata S.On the radial vibration of ball bearings(computer simulation)［J］.Bulletin of the JSME，1985，28:899－904.

［4］Tiwari M，Gupta K.Dynamic response of an unbalanced rotor supported on ball bearings［J］.Journal of Sound and Vibration，2000，238(5):757－779.

［5］Lioulios A N，Antoniadis I A.Effect of rotational speed fluctuations on the dynamic behaviour of rolling element bearings with radial clearances［J］.International Journal of Mechanical Sciences，International Journal of Mechanical Sciences，2006，48(8):809 －829.

［6］Patil M S，Mathew J.A theoretical model to predict the effect of localized defect on vibrations associated with ball bearing［J］.International Journal of Mechanical Sciences，2010，52(9):1193－1201.

［7］Liew A，Feng N，Hahn E.Transient rotor dynamic modelling of rolling element bearing systems［J］.Transactions of the ASME，Journal of Engineering for Gas Turbines and Power，2002，124(4):984 －991.

［8］Sopanen J，Mikkola A.Dynamic model of a deep-groove ball bearing including localized and distributed defects.Part 1:theory，Proceedings of the Institution of Mechanical Engineers［J］.Part K:Journal of Multi-body Dynamics，2003(217):201－211.

［9］Sopanen J，Mikkola A.Dynamic model of a deep-groove ball bearing including localized and distributed defects.Part 2:implementation and results，Proceedings of the Institution of Mechanical Engineers［J］.Part K:Journal of Multi-body Dynamics，2003(217):213 －223.

［10］Sawalhi N，Randall R B.Simulating gear and bearing interactions in the presence of faults:Part I.The combined gear bearing dynamic model and the simulation of localized bearing faults ［J］. MechanicalSystems and Signal Processing，2008，22(8):1924－1951.

［11］Mohsen N，Michael D.Dynamic modeling of rolling element bearings with surface contact defects using bond graphs［J］.ASME，Journal of Tribology，2011，133(1):1－12.

［12］何正嘉，訾艳阳，陈雪峰.机械故障预示中的若干科学问题［J］.振动与冲击，2010，29(S):144 －147.HE Zheng-jiia， Z1 Yan-yang，CHEN Xue-feng. Some scientific issues in mechanical fault prognosis［J］.Journal of Vibration and Shock，2010，29(S):144 －147.

［13］王晓冬.基于多小波变换的机械设备早期故障与复合故障诊断研究［D］.西安:西安交通大学，2009.

［14］Harris T A.Rolling bearing analysis［M］.4th ed.New York:John Wiley and Sons，2001.

［15］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.

［16］Hong H B，Liang M.Fault severity assessment for rolling element bearingsusing the Lempel-Ziv complexity and continuous wavelet transform ［J］.Journal of Sound and Vibration，2009，320(1 －2):452 468.