考虑滑动和热膨胀的球轴承局部故障动力学建模

雷宏卫， 成建联

(长安大学 道路施工技术与装备教育部重点实验室，西安 710064)

滚动轴承是各种旋转机械的重要零件，若其发生故障，会导致设备不能正常运转。轴承的动力学特性可以反映轴承是否存在故障[1]。磨损、剥落等常见局部故障会影响轴承寿命，当发生这些故障时，轴承的动力学特性也会发生变化。因此对局部故障轴承的动力学特性进行研究很有意义[2]。Mcfadden等[3]建立了恒定径向载荷下考虑轴承形状、转速、载荷分布以及振动衰减的内圈局部故障滚动轴承动力学模型。Patel等[4]考虑多种因素包括各组件质量和不同缺陷状态，建立球轴承动力学二自由度模型，模型结果表明具有多个局部故障的轴承相比单一故障轴承的振动谱具有更高的幅值。Patil等[5]将轴承简化为弹簧-阻尼系统，利用赫兹接触理论计算滚动体与滚道的接触力，对轴承的动力学特性进行研究，发现外圈局部故障所引起的振动幅值高于内圈局部故障和滚动体故障引起的振动幅值。Nakhaeinejad等[6]基于多体动力学建立了考虑陀螺效应、离心力、弹性挠度、滑动的局部故障动力学模型。曹宏瑞等[7]研究了外圈滚道表面状态和滚动体与滚道间隙的变化对轴承动力学特征的影响，发现轴承出现磨损故障后，其振动响应时域信号波形无明显规律；若间隙增加，振动的幅值以及磨损也会增加。牛蔺楷等[8]建立六自由度轴承局部表面损伤动力学模型，充分考虑了滚动体三维运动和相对滑动，研究了滚动体在通过局部故障时的动力学特性。徐可君等[9]研究了滚动体存在缺陷的中介轴承采用不同支承形式时的振幅特征。田晶等[10]忽略滚动体滑动，建立了包含弹流润滑和时变位移因素的动力学模型，结果表明缺陷尺寸增大会增加中介轴承的振动响应。

轴承运转时内部各组件会由于相互摩擦而产生热量，这些热量使得轴承产生不均匀的膨胀，影响轴承性能。Ma等[11]建立精确计算球轴承发热速率的模型，分析了多种因素对轴承温度的影响。Hannon等[12-13]建立滚动轴承传热模型，研究了轴承以及轴承座、转轴的温度场。李欢锋等[14]使用拟静力学的方法建立角接触球轴承生热模型，研究了内沟曲率和外沟曲率对轴承生热的影响。但多数对轴承温度场的研究没有包含轴承温度对其动力学特性的影响。

综上所述，对局部故障球轴承的动力学特性和轴承温度场已有很多研究成果，但大都缺少对轴承温度变化引起的组件热膨胀对其动力学特性影响的考虑。本文综合考虑滚动体经过缺陷区域不同状态的时变刚度、时变位移，以及轴承转动导致各组件摩擦产生热量引发的热膨胀对局部故障轴承动力学特性的影响，建立二自由度动力学方程，并利用试验对其进行验证。

1 轴承弹簧-阻尼模型

为简化研究，本文对研究对象做出如下假设[15]：

(1) 轴承的故障为外滚道单一局部剥落，剥落形式为与滚道等宽的矩形凹槽，且凹槽具有一定深度；

(2) 滚动体均匀分布在滚道中且互不接触；

(3) 滚动体与滚道之间满足赫兹接触理论；

(4) 接触力只在径向发生作用；

(5) 轴承的外圈与轴承座、内圈与转轴皆为刚性连接，且转轴质量包含内圈。

在上述条件下，将球轴承系统简化为弹簧-阻尼模型，其中外加载荷Fy作用在转轴上。

2 滚动体的滑动

如图1所示，将第i个滚动体的角度位置记为Ψi。在轴向载荷Fr作用下，图中θr范围内为滚动体承载区。Ψi可以表示为[16]

(1)

图1 轴承受载及位置关系Fig.1 Loading and angle relationship of bearing

式中:ωs表示转轴的角速度;ωc为保持架角速度;Z为滚动体数量;Db是滚动体直径;Dp为轴承节圆直径；α表示轴承接触角，深沟球轴承的α=0。Ψ0表示第1个滚动体的角度位置。

当滚动体刚受载时会发生相对滑动，此时Ψi表示为

(2)

式中,γ为在[-θslip,θslip]上的均匀分布函数，其中0.01 rad<θslip<0.02 rad。

3 时变刚度模型

通过Hertz接触理论可得滚动体与滚道的接触力[17]

F=Kδn

(3)

式中,K为滚动体与内、外圈等效接触刚度；δ表示径向位移;n是载荷-变形系数，球轴承取n=1.5。

等效接触刚度K如式(4)所示

(4)

其可以由滚动体与内圈滚道的接触刚度Ki和滚动体与外圈滚道的接触刚度Ko计算得出，如式(5)所示

(5)

式中:∑ρi为滚动体与内圈接触的曲率之和;∑ρo为滚动体与外圈接触的曲率之和；Eeq为等效弹性模量；ξ2为第二类椭圆积分；ξ1为第一类椭圆积分；e为椭圆度参数。外圈滚道的局部剥落会导致剥落区附近一定范围内的接触刚度发生变化[18]，记θch为剥落区前后接触刚度发生变化的角度范围，βen、βex分别为滚动体进入和退出外圈滚道剥落区角度位置，βmid为剥落区中心在滚道上的角度位置，2θsp为βen、βex之间的角度，如图2(a)。

(a)

(b)图2 剥落区附近接触刚度Fig.2 Contact stiffness values near the defect

当滚动体运转至βen-θch<Ψi<βen范围时，其与外圈滚道之间的接触刚度Ko进行修正后得到新的接触刚度Ke1为

Ke1=Ko[1-R[1+(Ψi-βen)/θch]S]

(6)

当滚动体逐渐远离退出边缘，即βex<Ψi<βex+θch时，修正后的接触刚度Ke2可表示为

Ke2=Ko[1-R[1+(βex-Ψi)/θch]S]

(7)

式中：R为在θch范围内的载荷-变形接触比例数;S为形状系数。随着滚动体在θch内滚动，不同位置处对应的Ke值不同。

当滚动体处于βen<Ψi<βex范围内时，如图2(b)所示，有如下关系

Ke=(1-R)Ko

(8)

4 时变位移模型

滚动体经过剥落区时其位移会随时间发生变化[19]，如图2(a)所示。设H为剥落区深度，Hr为滚动体通过剥落区时的最大径向位移，其值如式(9)所示。L为剥落区长度，Hd为滚动体时变位移,Do为轴承的外圈滚道直径。图2(a)所示情况剥落区的深度大于滚动体通过此区域时的最大径向位移，图2(b)所示情况剥落区深度等于滚动体最大径向位移。

(9)

图2(a)情况下滚动体经过剥落区域时会与凹槽两侧边缘接触。

根据几何关系，其时变位移Hd表示为

(10)

图2(b)所示情况下，滚动体先后经过接触剥落区左侧、底部、右侧，再次进入正常滚道。Hd为

(11)

式中,θen为滚动体与剥落区域底部刚接触时角度位置与βen之间的夹角，如图2(b)所示，其值为

(12)

5 动力学模型的建立

5.1 轴承热膨胀分析

球轴承在高速运转时，滚动体的公转搅油、自旋以及其与滚道的差速滑动等都会产生热量[20]。由于轴承各组成部分的结构、散热系数有所不同，内部产生的热量会引起各部分的不均匀膨胀，导致热诱导载荷的产生。因此，在对轴承进行动力学研究时热膨胀因素不容忽视。

轴承常用的热分析方法有整体法和局部法两种。使用整体法得出的计算结果精度低，且难以对实际工况准确预测，本文采用精度高的局部法对球轴承进行热分析。球轴承的严重故障通常发生在热诱导载荷上升的初期阶段。在轴承刚开始运转时，轴承内部产生的摩擦热有很大比例转移到滚动体上，此时滚动体的温度会明显高于内滚道和外滚道。轴承工作一段时间之后，一方面热量通过轴承座和转轴散失，另一方面轴承座吸收了轴承的热量而发生膨胀，轴承中的热诱导载荷逐渐减小[21]。

由于内、外圈的散热条件比滚动体好，而且滚动体的体积、质量小，因此在热诱导预紧力上升的初期阶段，滚动体的热膨胀程度最大。综上，本文在建立动力学模型时主要考虑滚动体热膨胀对轴承动力学特性的影响。使用局部法计算轴承生热率的方法如下所示[22-23]

(13)

式中:Hfb表示滚动体与滚道的摩擦生热率，其中包含了差动滑动摩擦、陀螺转动摩擦和自旋摩擦三种类型。τybi为滚动体在接触点的剪切力，Vybi表示滚动体相对滑动速度，ωti、Fxbi为滚动体的自转角速度和摩擦力，Mbi、ωri为滚动体的自旋摩擦力矩和角速度。其中下标y,x分别表示沿接触椭圆的短轴和长轴方向，b取值0时表示内圈，取值1时表示外圈。

(14)

式中:Hli为滚动体公转搅油生热率;Hc、Hi分别表示保持架与引导面和滚动体的摩擦生热率;Flui为润滑油对滚动体摩擦力;ωsi、ωo分别为滚动体公转和外圈的角速度。Hc表示保持架与引导面的摩擦生热率，Dc为引导面直径，Fsl是滑动摩擦力，csl为滑动系数，ωb为套圈角速度。Hi为滚动体与保持架的摩擦生热率，μ、Qi分别为滚动体与保持架的摩擦因数和接触载荷，ωi为滚动体角速度。忽略热辐射的情况下，热传导热阻Rh和对流热阻Rc有如下关系式

(15)

式中:Lh、Ah为导热特征长度和换热面积;Kh为材料的热导率;he为对流系数。he可表示为

(16)

式中:de为对流换热特征长度;Nu为努塞尔数;λa为流体热导率。轴承的瞬态热平衡方程为

(17)

式中:T表示待求节点温度;Tj是待求节点相关的各节点温度;Rjp为节点j，p间热阻;Hp是待求节点的生热量;mp为节点相关的质量;Cp表示材料比热容。如图3为热网络模型示意图。根据热网络模型和瞬态热平衡方程，建立18个节点的微分方程组。采用四阶龙格库塔法求解微分方程组，得到轴承瞬态温度特性。图4为热膨胀位移求解流程。

图3 轴系的热网络节点分布Fig.3 Thermal network diagram of the shaft-bearing system

图4 热膨胀求解流程Fig.4 Thermal expansion solution process

5.2 轴承热膨胀位移计算模型

球轴承滚动体径向热膨胀位移可表示为[24]

ub=βbTbDb

(18)

式中:ub代表滚动体径向热位移;βb是滚动体热膨胀系数;Tb表示滚动体温度。

5.3 轴承动力学模型

χi为滚动体与滚道的接触参数，定义如式(19)

(19)

考虑滑动和热膨胀的球轴承局部故障动力学模型如式(20)

(20)

6 试验验证

试验轴承型号SKF6205，参数如表1所示。设定步长为Δt=3×10-5s，径向载荷Fx=0，Fy=400 N。初始位移x0=10-6m，y0=10-6m。为了验证本文轴承温度和热膨胀计算方法的正确性，试验测量轴承外圈和内圈的温度，在轴承座上沿径向加工通孔以安装温度传感器测量外圈温度，传感器量程0～750 ℃，分辨率0.05 ℃。利用红外测温法测量内圈温度，所使用传感器量程为-50 ℃～800 ℃，分辨率为0.1 ℃。

表1 SKF6205轴承主要参数Tab.1 The parameters of SKF6205 bearing

将试验测得温度结果与使用本文方法模拟所得轴承外圈和内圈温度进行对比。如图5所示，为试验和模拟外圈温度结果。图6所示为试验和模拟的内圈温度结果。对比可知两者误差小于5%。图7所示为模拟热诱导载荷预测曲线，与文献[20]中的预测趋势一致。上述结果证明了本文轴承热计算方法的可行性。

图5 轴承外圈温度Fig.5 The temperature of outer ring

图6 轴承内圈温度Fig.6 The temperature of inner ring

图7 热诱导载荷预测Fig.7 Prediction of thermally induced load

使用电火花在两个轴承外圈上分别加工出与滚道等宽的长度为0.5 mm和4.0 mm，深度为0.2 mm的矩形凹槽。电机主轴转速为2 500 r/min，轴承的特征频率理论值计算公式如式(21)所示[25]，表2为轴承特征频率理论值。

(21)

表2 轴承的特征频率Tab.2 Characteristic frequency of bearing

式中:fo为保持架通过外圈的频率;fr为转频;外圈剥落时的振动特征频率为Zfo。

图8、图9所对应的轴承剥落凹槽长度为0.5 mm。图8(a)为本文模型得到的振动信号时域图，图8(b)为将时域信号经过EEMD法处理过后的频谱图，从中可清晰识别故障频率及其二倍频。图9为试验所得时域图和频谱图，可以发现试验所得的频谱中噪声较大，这是由轴承的加工、安装误差等因素导致的。

(a)

(b)图8 0.5 mm剥落轴承仿真时域与频域Fig.8 Simulation of 0.5 mm defect bearing in time domain and frequency domain

(a)

(b)图9 0.5 mm剥落轴承试验时域与频域Fig.9 0.5 mm defect bearing experiment time domain and frequency domain

图10、图11对应凹槽长度为4.0 mm，可以其从频谱图中识别故障频率点及其二倍频点。随着凹槽尺寸

(a)

(b)图10 4.0 mm剥落轴承仿真时域与频域Fig.10 Simulation of 4.0 mm defect bearing in time domain and frequency domain

(a)

(b)图11 4.0 mm剥落轴承试验时域与频域Fig.11 4.0 mm defect bearing experiment time domain and frequency domain

的变大，其信号幅值也增大。这表明轴承局部剥落尺寸对轴承寿命有重要影响，剥落尺寸越大尺轴承受到的冲击越大，寿命越短。模拟信号的故障特征频率与试验值误差小于3%。

7 结 论

(1) 考虑了滚动体相对滑动、时变接触刚度和热膨胀因素，建立了能更准确反映实际振动特性的二自由度轴承动力学方程。

(2) 所建立模型包含与剥落尺寸相关的滚动体时变位移。分析结果表明轴承滚道剥落尺寸越大，所引起的冲击振幅越大。说明剥落尺寸的大小对轴、承寿命有重要影响。

(3) 本文所建立动力学模型的故障频率与试验相比，误差小于3%。验证了所建立模型的可行性。