黄贤振, 朱会彬, 姜智元, 姜 睿
(1. 东北大学 机械工程与自动化学院, 辽宁 沈阳 110819; 2. 东北大学 航空动力装备振动及控制教育部重点实验室, 辽宁 沈阳 110819)
角接触球轴承因其寿命长、能耗低,能同时承受轴向和径向载荷,被广泛用于旋转机械中.然而,轴承长期服役于高速工况下,除了常见的疲劳[1]破坏之外,更多的是轴承出现打滑现象,严重的打滑行为会造成沟道表面摩擦磨损,同时会引起主轴振动和严重的轴承温升,降低轴承的使用寿命和可靠性.因此,为了保证机械设备的可靠性和安全性,有关轴承打滑的研究极为重要.
国内外学者对轴承打滑行为进行了大量研究.在理论分析方面,Harris[2]考虑了滚动体所受的离心力、接触应力和摩擦力,建立了预测滚动轴承打滑行为的力学分析模型.Jain等[3]建立了滚动轴承的动力学模型,研究分析了轴、径向联合载荷作用下的轴承打滑特性.Laniado-Jacome等[4]通过有限元的方法对滚动体与内外沟道之间的相对滑动进行了仿真分析.陈渭等[5]研究分析了涡动状态下的打滑问题.涂文兵等[6]建立了考虑轴承非线性因素的动态模型,研究了不同工况对滚动体打滑的影响.在实验分析方面,Hirano[7]通过监测轴承钢球磁通量的变化,确定出轴承打滑判据.Xu等[8]基于主轴热分析实验得到不同转速下轴承避免打滑的最佳预紧力,通过实验验证了Hirano[7]判据的可行性.Dong等[9]和Oktaviana等[10]通过Hirano[7]判据研究分析了避免轴承打滑的最小预紧力.王海同等[11]建立大尺寸球轴承拟静力学方程,基于Hirano[7]判据得出了轴承打滑的临界转速.Zheng等[12]通过光纤传感器监测脉冲信号对滚动体的通过频率,实现对轴承打滑的预测.Zhan等[13]采用弱磁探测器检测滚子和内沟道转速来预测轴承的打滑行为.上述理论和实验研究中,通常认为轴承结构参数和材料参数是确定的、无误差的.然而,由于制造和加工的影响,轴承的结构参数和材料参数存在一定的随机不确定性[14],这就导致球轴承打滑特性的分析产生较大的偏差.
本文综合考虑随机因素对轴承参数的影响,将轴承拟静力学分析模型和Hirano轴承打滑的判定依据相结合,以轴承是否打滑为判别条件,建立轴承打滑的极限状态方程,提出了一种球轴承打滑行为的可靠性分析模型.应用MCS方法和Kriging方法构造高精度的响应面函数进行可靠性灵敏度分析,以确定轴承结构参数对打滑现象的影响程度.为减少或避免轴承发生打滑引起早期失效提供理论依据.
图1表示角接触球轴承的几何结构及参数,图2显示了轴承中每个球的相对角位置(方位角),其中第j个滚动体方位角为φj,其表达式为
φj=2π(j-1)/Z.
(1)
图2 轴承滚动体的方位角
球轴承运动过程中,在联合载荷作用下,滚珠沿着轴向和径向相对移动,如图3所示,对轴承几何位置进行分析.
图3中,变形后任意方位角φj处,内外沟道曲率中心相对轴向距离A1j和径向距离A2j分别为
A1j=BDsinα0+δa+θRicosφj,
(2)
A2j=BDcosα0+δrcosφj.
(3)
任意方位角φj处的角度方程为
(4)
(5)
(6)
(7)
式中:fi,fo分别为内、外圈沟道曲率半径系数;X1j,X2j为滚珠中心与外圈沟道曲率中心之间的轴、径向距离;δa,δr,θ分别为变形后轴、径向位移和角位移;δij,δoj分别为内外圈的接触变形;BD为内外圈曲率中心初始距离.
图3 滚动体中心与沟道曲率中心的相对位置
根据几何位置关系,以第j个球为研究对象,确定沟道接触的变形几何相容方程为
(8)
(A1j-X1j)2+(A2j-X2j)2- [(fi-0.5)Dw+δij]2=0 .
(9)
轴承滚动体运动过程中,运动方式复杂.如图4所示,以第j个滚动体为研究对象,对其建立受力平衡方程.
图4中,λij和λoj分别为滚珠和内外圈沟道之间的控制参数,高速角接触球轴承运动中,根据外沟道控制理论,即滚珠在外沟道只发生纯滚动,因此控制参数取λij=0,λoj=2;Qij,Qoj分别为滚珠和内外圈的接触载荷,根据赫兹接触理论,接触载荷计算公式为
(10)
式中,Kij和Koj分别为内外圈的载荷变形系数.
滚珠受到离心力Fcj,陀螺力矩Mgj分别为
(11)
(12)
式中:dm为中心圆直径;ωm,ωR分别为滚珠的公转角速度、自转角速度;βj为第j个滚珠的螺旋角.其求解方式详见文献[15].
根据图4中的受力平衡关系,滚珠的受力平衡方程为
(13)
(14)
图4 滚动体受力分析
轴承在轴、径向载荷作用下,处于受力平衡状态,建立轴承内圈受力平衡方程为
(15)
(16)
(17)
式中:Ri=dm/2+(fi-0.5)Dwcosα0;Fa为轴向载荷;Fr为径向载荷;M为轴承弯矩.
应用Newton-Raphson迭代法求解上述列出的非线性方程组,可计算出内外圈的接触角αij,αoj,接触载荷Qij,Qoj等动态参数.
角接触球轴承服役于高速工况时,随着转速的增加,由于离心力和陀螺力矩的作用,滚珠与内圈沟道之间的赫兹接触面积减小.当转速提高到一定速度时,轴承内圈与滚珠之间的摩擦力小于陀螺力矩的作用,滚珠与内圈沟道失去接触时,滚动体就会发生打滑,导致摩擦力增大,产生大量的热量,容易引起轴承的早期失效.Hirano[7]确定出角接触球轴承打滑的判定依据式(18),当轴承动态参数满足此判据时,滚珠会发生打滑现象.
(18)
式中:SF为打滑系数[10],决定球轴承打滑行为的发生;Qi为内沟道与球的接触载荷;αi为内沟道与球的接触角;Fc为球的离心力.
根据上述研究,考虑随机因素对轴承结构参数的影响,通过拟静力学分析模型结合轴承打滑判据式(18),以轴承是否打滑为判别条件,即可以建立极限状态函数为
Z=g(X)=S(X)-10 .
(19)
式中:X为影响轴承打滑的参数随机变量;S(X)为数值模拟方法SF的输出值.
X=[Do,Di,Dw,ri,ro,α0]T.
函数Z>0时,表示轴承处于不发生打滑的可靠状态,称为可靠性概率或可靠度,用Pr表示.对随机变量的概率密度函数积分[16]可得到Pr为
(20)
然而,通过Newton-Raphson迭代法计算大量的轴承动态参数数据,采用Monte-Carlo方法[17]计算可靠度需要耗费很长时间.因此提出用Kriging模型来代替上述过程,提高计算效率[18].
Kriging模型是一种半参数化插值模型,可以通过已知点的数据预测出未知点的数据[19-20],在解决强非线性问题中具有显著优势,Kriging模型的具体形式为
y(x)=F(β,x)+z(x)=f(x)Tβ+z(x) .
(21)
式中:f(x)为样本点x的多项式函数;F(β,x)为线性回归模型;β为回归系数.f(x)反映模型的全局近似;z(x)为随机分布函数,反映模型的局部近似,影响整个模型的准确性,且z(x)服从正态分布,均值μ=0,方差为σ2,其协方差为
cov[z(xi),z(xj)]=σ2R(θ,xi,xj) .
(22)
式中:xi,xj为任意两个样本点;R(θ,xi,xj)为带有参数θ的相关函数,反映样本点之间的空间相关性,直接影响到模型的准确性,选用计算效果最好的高斯函数作为相关函数.
已知点的响应值Y=[y1,y2,…,ym],采用线性组合向量c来估计未知点x的响应值,即
(23)
估计值与真值之间的偏差为
(24)
误差为0时,可保证估计值的无偏性,即
FTc-f(x)=0 .
(25)
预测的均方误差为
(26)
式中:R为相关函数矩阵;r为相关向量.
为了模型的准确性,要使φ(x)为最小值,因此求线性组合向量c转化为如下最优化问题.
findc. minφ(x).s.t.FTc=f(x).
(27)
通过求解得到Kriging模型的估计值为
(28)
β*和σ2采用极大似然估计可得估计值为
β*=(FTR-1F)-1FTR-1Y,
(29)
(30)
根据上述分析,通过输入变量样本和输出响应值构建完成Kriging模型,提高计算的效率.
失效概率Pf对基本变量xi的分布参数θxi的偏导数定义为可靠性灵敏度:
(31)
为了得到轴承各结构参数对可靠度的影响程度,需要进行无量纲处理,得到失效概率对第i个变量的均值μxi和标准差σxi的灵敏度系数:
(32)
(33)
应用Kriging方法进行可靠性灵敏度分析,计算流程如图5所示.
为了验证球轴承打滑预测方法的有效性,以
图5 Kriging可靠性分析流程图
角接触球轴承B7007C为例,以转速n=5 000 r/min为研究对象,与文献中计算和实验结果进行对比,轴承的原始参数如表1所示.
表1 角接触球轴承的原始参数
表2显示了转速n=5 000 r/min时,本文计算出轴承防止打滑的预紧力,与文献[8,10]计算值与实验得到轴承避免打滑的最佳预紧力比较,理论计算值接近实验值,验证了本文的轴承打滑预测方法的有效性和准确性.
根据上文对比已经验证了轴承打滑预测方法的正确性,本文以B7005C轴承作为研究对象.
轴承生产制造过程中,由于受到大量独立因素的影响,轴承结构参数具有随机不确定性,这些随机参数一般符合正态分布,参数的标准差可以由变差系数c决定,本文假设c=0.003.轴承的随机参数对应的均值和变差系数如表3所示.
表2 轴承避免打滑的最佳预紧力
表3 轴承参数随机变量
为了研究不同工况对轴承打滑行为的影响,计算出不同转速和预紧力的打滑系数SF,计算结果如图6所示.通过SF=10建立一个临界平面区,分出打滑区域和非打滑区域.
图6 轴承打滑临界平面
轴承转速n=10 000 r/min,预紧力Fa=400 N条件下,利用拉丁超立方方法随机抽取60组轴承结构参数样本数据,作为Kriging模型的输入量.通过拟静力学模型求得打滑系数SF作为Kriging模型的输出量,建立完成Kriging模型,抽取30组随机变量测试样本数据,分别代入轴承拟静力学分析模型和Kriging模型中计算,得到30组SF的对比情况,如图7所示.图8为Kriging模型的预测误差,经过Kriging模型计算的SF与轴承拟静力学分析计算结果之间的误差非常小,验证了建立的Kriging模型的有效性.
图7 Kriging计算值与确定性计算值的对比
图8 Kriging计算值与确定性计算值的相对误差
根据本文建立的极限状态函数,分别采用Monte-Carlo方法和Kriging方法进行可靠度的计算,如图9和图10所示.两种方法具有良好的一致性,进而验证了Kriging方法的准确性.
为了探究预紧力对轴承打滑可靠度的影响,在转速n=10 000 r/min条件下进行可靠性分析,得到打滑可靠度随预紧力的变化规律.由图6可知,打滑系数SF=10时,计算得到防止打滑所需的预紧力为370 N,图9表明其可靠度为0.519,可知由于随机因素的影响,发生打滑的概率仍然较大;当预紧力达到440 N时,可靠度为0.996,即轴承发生打滑概率较低.当转速一定时,随着轴向预紧力的增大,可靠度增加.
为了探究转速对轴承打滑可靠度的影响,在预紧力Fa=400 N条件下进行可靠性分析,得到打滑可靠度随转速的变化规律.由图10可知,当转速低于9 500 r/min时,轴承打滑可靠度为1,即轴承没有发生打滑行为,当转速超过9 500 r/min,可靠度开始下降,轴承开始出现打滑现象,达到12 000 r/min时,可靠度降为0,随着转速的增加,可靠度降低.
图9 可靠度随预紧力的变化
图10 可靠度随转速的变化
为了探究轴承随机变量的变化对轴承打滑的影响程度,利用建立的Kriging模型对随机变量进行了可靠性灵敏度分析,如图11和图12所示.
图11 均值灵敏度分析
由图11和图12可知,随着内外圈沟道直径以及滚动体直径、接触角的增加,打滑可靠度降低;随轴承内外圈沟道曲率半径增加,打滑可靠度增加.
图12 标准差灵敏度分析
1) 本文建立了球轴承拟静力学分析模型,并基于Hirano试验打滑判据,考虑轴承参数的随机性更加符合轴承实际工况,提出一种角接触球轴承打滑行为的可靠性分析模型,对防止轴承发生打滑引起早期失效、提高轴承的可靠性具有重要意义.
2) 本文采用Kriging方法进行了可靠度的计算及可靠性灵敏度分析.轴承打滑可靠度会随内外圈沟道直径以及滚动体直径、接触角的增加而降低,随轴承内外圈沟道曲率半径增加而增加.其中,轴承滚动体直径的变化对轴承打滑现象影响最大,轴承内外圈沟道直径和曲率半径的变化对其影响次之,接触角的变化对其影响较小.