歪斜状况下滚子轴承的接触应力求解与分析

李伟伟 陈晓阳 沈雪瑾 张小玲

1.上海大学，上海，200072 2.苏州轴承厂有限公司，苏州，215000

0 引言

滚子轴承具有承载能力大、径向尺寸小等优点，广泛应用在重载及径向空间受限的场合。疲劳点蚀是滚子轴承最典型的失效方式。研究发现，滚子轴承的疲劳寿命跟接触应力的若干次方成反比。赫兹于1881年给出了滚动轴承点接触下完整的理论解，而对于线接触仅给出了理想状况下的部分解。通常滚子与滚道的接触属于有限长线接触，对于此类问题的求解主要有三种方法：一是经典数值解法；二是基于变分原理的方法；三是基于有限元、边界元等现代数值方法。其中经典数值解法因具有算法简便、求解精度高等优点而应 用 最 为 广 泛。Hartnett等[1]、Chiu 等[2]将Boussinesq理论同影响系数法结合起来，求解了一般表面轮廓的接触问题，分析了正载及偏载工况下滚子与滚道的接触应力。Ahmadi等[3]改进了Hartnett的方法，将接触区域沿宽度方向划分为偶数个单元，避免了方程组的病态。文献[4-5]利用切片法分析了滚子与滚道在正载及偏载工况下的接触应力，该方法忽略了片与片之间的联系。文献[6-7]在切片法的基础上，考虑了片与片之间的联系，求解了滚子与滚道的接触问题。陈晓阳[8]将Boussinesq理论同影响系数法相结合，分析了对数轮廓滚子与滚道的接触应力。毛月新[9]分析了偏载工况下滚子与滚道的接触应力。文献[10-13]对滚子轴承歪斜工况产生的原因及影响因素进行了分析。文献[14-15]从工艺上提出了减小滚子歪斜角的方法。但关于滚子轴承在歪斜状况下滚子与滚道的接触问题鲜有报道，本文将针对此问题展开讨论。

1 歪斜工况下滚子与滚道的接触力学模型

根据弹性接触理论，在歪斜工况出现前，滚子与滚道的接触模型如图1所示。当滚子相对于滚道产生θ角的歪斜时，滚子与滚道的接触模型如图2所示。M、N分别是滚子和滚道外表面上的一点，其坐标分别为（x，y，z1）、（x，y，z2），受载前滚子与滚道在O点相接触，受载后起始接触点O扩展为一接触区域Ω，在接触区域上压力以函数p（x，y）分布，如图3所示。弹性体表面有相同坐标（x，y）的两点M和N在z轴方向上的原始距离为z1＋z2，由此可得位移平衡方程如下[8]：

图1 滚子与滚道正载下的接触模型

图2 滚子与滚道歪斜角为θ时的接触模型

图3 歪斜角为θ时，滚子与滚道受载后的应力分布示意图

式中，δ为两弹性体的弹性趋近量；ω1、ω2分别表示两弹性体在Z方向上的变形。

根据Boussinesq弹性体半空间体力—变形解，在接触区域Ω内可以得到求解一般光滑弹性体接触问题的主导方程——Fredholm积分方程：

式中，p（x′，y′）为光滑弹性体表面坐标为（x′，y′）点的压力；E1、E2分别为滚子与滚道的弹性模量；ν1、ν2分别为滚子与滚道的泊松比。

同样，在接触区域Ω内可得力的平衡方程：

式（2）和式（3）就构成了求解歪斜工况下滚子与滚道接触问题的基本方程。由于式（2）中积分核是奇异的，且积分域本身也是未知量，所以积分求解式（2）是非常困难的。计算机运算速度、精度的提高及数值计算方法的发展为数值求解有限长线接触问题提供了条件。

2 数值解法

利用结构分析中的影响系数法，首先设定一个矩形接触区域Ω1，且满足Ω1＞Ω，然后将矩形接触区域Ω1均匀划分成n个矩形单元。假设接触压力在各个矩形单元上是均匀分布的，以pj表示单元j（j＝1，2，…，n）上的接触压力，并以单元i（i＝1，2，…，n）的中心（xi，yi）处弹性体表面的原始距离作为该单元的表面函数离散值Si，则在接触区域Ω1内离散处理式（2）、式（4）得[16]

式中，Fij为影响系数，它是作用在单元j上的压力使单元i中心处产生z方向上的变形的系数；ag、bg分别为中心点坐标为（xj，yj）的网格的半长和半宽；Aj为矩形单元j的面积。

在接触区域内求解Fij，式（5）、式（6）可用矩阵形式表示为

式中，s1，s2，…，sn为编号为1，2，…，n的网格所对应的表面函数。

分析式（8）～式（11）可知：当滚子与滚道的几何形状、材料常数和外力pj已知的情况下，通过n＋1个方程求解n＋1个未知数，方程组有唯一解。据上述分析，用F ORTRAN计算机语言编程求解式（8），程序流程如图4所示。

3 程序正确性验证

3.1 与点接触赫兹解作对比

当滚子与滚道的歪斜角θ＝90°时，滚子与内圈滚道的接触可看作是点接触，可以利用赫兹理论解与程序计算结果作对比来验证程序结果的正确性。下例中，滚子半径R＝5mm，滚子长度L＝18.1mm，内圈滚道半径Ri＝5mm，外载荷Q＝600N，滚子与滚道的弹性模量和泊松比均分别为206GPa、0.3。

图4 程序求解流程图

由表1可知：假设θ＝90°，通过和点接触标准赫兹解比较可知，程序计算结果误差小于2%。

表1 赫兹理论解与程序解结果对比

滚子轴承在实际工作过程中，滚子与滚道的歪斜状况通常受到保持架及挡边的限制，歪斜角不可能达到90°，上述例子仅是为了验证程序的准确性而假设的工况。

3.2 程序解与有限元分析结果的对比

利用Abaqus有限元软件，建立滚子与滚道歪斜的三维模型，为了在提高计算效率的同时又不降低计算精度，仅取滚子与滚道的接触部分进行建模。程序求解时为了节省计算机运算时间，使网格的长度方向跟滚子长度方向相同，这样可使求解区域减小，有利于计算速度及计算精度的提高。以某厂HK1212型滚子轴承为例，滚子母线轮廓为直母线，R＝1mm，L＝7.3mm，Ri＝6mm，Ro＝8mm，θ＝3°，滚子与滚道的弹性模量均为206GPa，泊松比均为0.3。

由图5、图6可知：当滚子与内外圈滚道间的歪斜角θ＝3°时，程序计算结果和有限元分析结果应力分布趋势一致，各部分误差均在10%以内。

滚子与内圈滚道之间歪斜角θ＝90°时与赫兹点接触结果的对比及滚子与内外滚道之间歪斜角θ＝3°时与有限元结果的对比验证了程序结果的正确性。

4 分析

4.1 不同歪斜角下滚子与滚道的接触应力分布

图5 θ＝3°时滚子与内圈之间的接触应力分布

图6 θ＝3°时滚子与外圈之间的接触应力分布

滚子轴承参数为：R＝5mm，L＝25mm，Ri＝10mm，Ro＝20mm，滚子母线为最佳对数母线。当外载荷Q＝20kN时，分别计算θ＝1°和θ＝2°时滚子与滚道的接触应力分布，结果如图7、图8所示。

从图7、图8可知：歪斜工况下滚子与滚道接触时会产生“歪斜效应”，即滚子与内圈接触时，接触应力在接触副中部增大，两端减小；滚子与外圈接触时，接触应力在接触副两端增大，中部减小。随着歪斜角的增大，歪斜效应逐渐明显，且滚子与外圈接触时歪斜效应更加明显。由于随着歪斜角的增大，滚子与滚道由线接触逐渐转为点接触，因此才会出现逐渐明显的歪斜效应。

图7 θ＝1°时滚子与滚道的接触应力分布

图8 θ＝2°时滚子与滚道的接触应力分布

4.2 不同外载荷作用下滚子与滚道间的接触应力分布

滚子轴承参数为：R＝5mm，L＝25mm，Ri＝10mm，Ro＝20mm，滚子母线为最佳对数母线。分别计算Q＝5kN、Q＝20kN作用下，θ＝1°时，滚子与滚道的接触应力分布。

图9 Q＝5kN、θ＝1°时滚子与内外圈滚道的接触应力分布

图10 Q＝20kN、θ＝1°时滚子与内外圈滚道的接触应力分布

从图9、图10可知：当歪斜角相同时，外载荷越小歪斜效应越明显。由于当滚子与滚道间的歪斜角一定时，随着外载荷的减小，滚子与滚道之间的接触由线接触逐渐转为点接触，沿滚道轴线方向的接触区域长度逐渐缩短，所以歪斜效应逐渐明显。重载时，歪斜角的存在也会使滚子与滚道间最大接触应力增加，出现应力集中，这将导致滚子轴承的疲劳寿命缩短。

4.3 不同内外圈滚道半径下的接触应力分布

滚子轴承参数为：Q ＝20kN，θ＝1°，R ＝5mm，L＝25mm。当Ri＝5mm、Ro＝15mm和Ri＝20mm、Ro＝30mm时，滚子与滚道的接触应力分布如图11、图12所示。

图11 Ri＝5mm、Ro＝15mm时滚子与滚道的接触应力分布

图12 Ri＝20mm、Ro＝30mm时滚子与滚道的接触应力分布

由图11、图12可知：歪斜角、外载荷及滚子半径一定时，随着内外圈滚道半径的减小，歪斜效应逐渐明显。由于当歪斜角、外载荷及滚子半径一定时，随着内外圈滚道半径的减小，沿滚道轴线方向的接触区域长度逐渐缩短，滚子与滚道的接触由线接触逐渐转为点接触，所以歪斜效应逐渐明显。

5 结语

根据弹性接触理论，建立了滚子相对滚道歪斜时的接触模型，结合影响系数法，利用FORTRAN语言编程求解了歪斜工况下滚子与滚道的接触问题。通过与点接触赫兹解及有限元结果的对比，验证了程序结果的正确性。通过分析歪斜角、外载荷及滚道半径对滚子与滚道接触应力分布的影响可知：歪斜工况下，滚子与滚道接触时会出现歪斜效应，即滚子与内圈接触时，接触应力在接触副中部增大，两端减小；滚子与外圈接触时接触应力在接触副两端增大，中部减小；歪斜效应随着歪斜角的增大、外载荷的减小及滚道半径的减小而逐渐明显。

[1]Hartnett M J.The Analysis of Contact Stresses in Rolling Element Bearings[J].Journal of Lubrication Technology，1979，101（1）：105-109.

[2]Chiu Y P，Hartnett M J.A Numerical Solution for the Contact Problem Involving Bodies with Cylindrical Surface Considering Cylinder Effect[J].Journal of Tribology，1987，109（7）：479-486.

[3]Ahmadi N，Keer L M，Mura T.Non－Hertzian Contact Stress Analysis——Nornal and Sliding Contact[J].International Journal of Solid and Structures，1983，19（4）：357-373.

[4]Kannel J W.Comparison between Predicted and Measured Axial Pressure Distribution between Cylinders[J].Journal of Lubrication Technology，1974，96（3）：508-514.

[5]Harris T A，Kotzalas M N.滚动轴承分析[M].罗继伟，译.北京：机械工业出版社，1997.

[6]Hoeprich M R.Rolling Element Bearing Contact Geometry Analysis[J].Tribology Transaction，1995，38（4）：879-882.

[7]陈家庆，周海，吴世勤.圆柱滚子轴承载荷分布的理论研究[J].轴承，2001（6）：8-11.

[8]陈晓阳.滚针轴承优化设计及CAD技术的研究[D].杭州：浙江大学，1989.

[9]毛月新.偏载滚子轴承凸度设计及疲劳寿命研究[D].上海：上海大学，2009.

[10]崔庆武.关于滚针歪斜问题的讨论[J].轴承，1980（1）：4-9.

[11]Harris T A，Kotzalas M N，Yu W K.On the Causes and Effects of Skewing in Cylindrical Roller Bearings[J].Tribology Transactions，1998，41（4）：572-578.

[12]焦育洁，吕新生，张锡昌.航空发动机主轴圆柱滚子轴承典型故障分析[J].轴承，2003（12）：29-30.

[13]杜辉.航空发动机主轴高速圆柱滚子轴承三维瞬态拟动力学分析[D].洛阳：河南科技大学，2005.

[14]栾春刚.燃汽轮发动机主轴轴承的设计要点[J].轴承，1993（1）：7-8.

[15]张玉环，李春平，陈倩.组件M型保持架窗孔引导部件及其宽度的确定[J].轴承，1997（4）：14-17.

[16]马家驹.线接触副凸度设计的研究——（Ⅰ）计算方法[C]／／第五界全国摩擦学学术会议.武汉，1992：167-174.