陈晓阳 周越 沈雪瑾 王志坚 张小玲 马纯青

(1.上海大学 机电工程与自动化学院，上海200072;2.苏州轴承厂有限公司，江苏 苏州215000;3.上海和锦滚子科技有限公司，上海201611)

径向尺寸小、结构简单、承载能力大等优点使得滚子轴承成为机械工业领域的重要零部件. 但由于安装误差、热变形、几何形状缺陷等原因，滚子轴承一般在偏、歪斜等工况下运行，使得轴承过早失效.Schaude［1］最先利用切片法对偏斜工况下滚子轴承接触应力展开研究，通过数值解法求解出了偏斜工况下滚子轴承接触应力分布. Kannel 等［2］通过切片法求解出了滚子接触副的接触应力分布，并对其结果进行了实验验证.毛月新等［3］利用影响系数法分析了偏斜工况下对数滚子与滚道的接触应力分布，并求解了滚子在偏斜工况下的最佳凸度量. 孙殿超等［4］利用快速傅里叶变换(FFT)和共轭梯度法求解了偏斜滚子接触副的接触应力，并阐述了该算法具有较快的计算速度和收敛性. Ye 等［5-6］利用有限元法分析了偏斜工况下滚子的接触问题. Harris 等［7］分析了滚子轴承中滚子偏斜和歪斜产生的原因及其影响.Yang 等［8］利用探针测出了圆锥滚子轴承在不同转速、载荷及润滑油下的歪斜角. 李伟伟等［9］采用影响系数法分析了歪斜工况滚子接触副的接触应力分布.但滚子轴承在偏斜和歪斜共同作用工况下的理想凸度滚子与内外圈滚道的接触问题的研究鲜有报道，文中将采用FFT 和共轭梯度法求解上述问题，并讨论对应工况下的滚子最佳凸度量.

1 偏歪斜工况滚子与滚道的接触分析

1.1 偏歪斜工况下滚子与滚道的力学模型

在实际工况下，由于轴的变形、轴承中心不对中以及轴承的径向游隙等原因，使得滚子接触副处于偏斜状态，如图1(a)所示.偏斜的存在使滚子在一端被挤压，压向引导挡边产生摩擦力，同时由于滚子与滚道的轴向间隙等原因使得滚子处于歪斜状态下工作，如图1(b)所示.根据静弹性接触理论，半径为R、长度为L 的滚子与滚道产生Ψ 度偏斜和θ 度的歪斜，M、N 分别是滚子和半径为RI的内圈外表面上的点，其坐标表示为(x，y，z1)、(x，y，z2).在外载荷F 作用下，起始接触点O 将延展为实际接触区域Ω，接触区域上的应力分布用函数P(x，y)表示，如图2 所示.

图1 滚子与滚道的接触力学模型Fig.1 Mechanics model describing the contact of roller and raceway

图2 滚子与滚道受载后的接触区域应力分布示意图Fig.2 Schematic diagram of stress distribution of roller and raceway in contact area after loading

根据Boussinesq 半无限体的力和变形解，在实际接触区域Ω 内可以求解出半无限弹性体接触问题的主导积分方程——Fredholm 方程［10］表示.为方便起见，给出其离散形式:

式中，Gij为接触间隙矩阵，Uij为表面弹性变形矩阵，Hij为接触体几何形状矩阵，Δ 为弹性趋近量.

偏斜工况下，弹性趋近量与接触副接触中心处的弹性趋近量δ0有以下关系:

歪斜工况下，接触体几何形状矩阵可表示为

其中，Sij为滚子与平面P 之间的几何间距，

Td为滚子凸度函数.

Tij为内滚道与平面P 之间的几何间距，

其中，

由此，即可求出弹性变形计算式及载荷和力矩平衡方程.若需求滚子与外滚道的接触体几何形状矩阵，只要改变Tij函数即可，这里不再赘述.

1.2 数值解法

1.2.1 FFT 求解弹性变形

在数值计算过程中接触副弹性变形的求解占据着主要的工作量.由弹性变形计算式是一卷积求和运算，可以利用FFT 进行快速求解，计算速度优于多重网格法［11］.具体计算可参阅文献［11-13］.

1.2.2 共轭梯度法求解接触压力

通过变分原理将接触问题转变为求解极值问题后，就可以通过共轭梯度法进行数值计算.该算法流程如图3 所示，具体算法的实现过程可参阅文献［4］.图中，Q、Q′分别为初始载荷和计算载荷，pij为压力矩阵元素，tij为迭代方向为修正量系数，ax、ay分别为网格的长度和宽度，ε0为迭代精度，ε1、ε2分别为压力和矩阵的迭代误差.

图3 计算流程图Fig.3 Flowchart for model solution

2 程序正确性验证和计算时间比较

为便于比较，文献［3，9］中的计算是文中的特例，即歪斜角为0°、偏斜角为0.02°和偏斜角为0°、歪斜角为2°. 在所用计算机为Intel Pentium Dual CPU E2160、计算环境为Visual Fortran 6.5 且网格数均为256 ×128 的条件下，分别采用两种方法计算，其结果如表1 和2 所示.

表1 文中计算结果与文献［3］计算结果对比1)Table 1 Comparison of calculated results obtained in the paper and from literature［3］

由表1、2 可知，文中计算结果与文献［3，9］计算结果的滚子与滚道接触区域半宽的相对误差都小于1.0%，最大接触应力的相对误差都小于2.0%，说明两种方法计算的接触应力和接触区域保持一致，证明了算法程序的准确性.特别是在计算时间方面，文中算法耗时大大减少，运算速度约是文献［3，9］计算的200 倍，从而证明了该方法具有高效的计算速度和收敛性.

表2 文中计算结果与文献［9］计算结果对比1)Table 2 Comparison of calculated results obtained in the paper and from literature［9］

3 各因素对接触应力分布的影响

滚子轴承参数:滚子副综合弹性模量E=226GPa，滚子半径R =5 mm，滚子长度L =25 mm，内圈滚道半径RI=10 mm，外圈滚道半径RO=20 mm，滚子母线为在正载工况下的最佳对数母线.

3.1 不同偏斜角对滚子接触应力的影响

当外载荷F=10 kN、歪斜角θ=0.5°时，计算在偏斜角Ψ=0°，0.02°，0.05°的作用下滚子与滚道的接触应力分布，如图4 所示.

图4 不同偏斜角滚子与滚道的接触应力分布Fig.4 Contact stress distribution at different misaligned angles between roller and raceway

由图4 可知:当滚子的歪斜角一定时，随着偏斜角的增大，滚子与内外圈滚道的“偏斜效应”逐渐明显，重载端的接触应力明显大于轻载端的接触应力，并且“偏斜效应”能明显加大“边缘效应”.当偏斜角增加到一定值时，轻载端的滚子与滚道将会脱离接触，其接触区域长度随着偏斜角的增加而减小.对比图4(a)和4(b)发现，当工况相同时，滚子与内圈的接触区域小于与外圈的接触区域，导致滚子与内圈的接触应力大于滚子与外圈的接触应力.

3.2 不同歪斜角对滚子接触应力的影响

当外载荷F=10 kN，偏斜角Ψ =0.02°时，计算滚子歪斜角θ= 0°，1°，2°作用下滚子与滚道间的接触应力，如图5 所示.

图5(a)为滚子偏斜角一定、歪斜角增加时滚子与内圈的接触应力分布. 研究发现:滚子发生歪斜时，滚子与内圈接触将会产生“歪斜效应”，并且当歪斜角逐渐变大时，滚子与内滚道的接触应力在接触副的中部变大，两端变小.“歪斜效应”逐渐明显，且一定的偏斜角能减弱滚子的“边缘效应”，滚子与内圈接触区域的长度方向呈外凸的轮廓曲线. 图5(b)为滚子偏斜角一定、歪斜角增加时的滚子与外滚道的接触应力分布. 可知:随着歪斜角的增大，滚子与外圈滚道的接触应力在接触副的两端变大，中部变小.滚子与外圈由于偏斜角会产生更加严重的“偏斜效应”和“歪斜效应”，并且滚子与外圈接触区域的长度方向呈外凹的轮廓曲线，这主要是由内外滚道的结构决定的，即内圈是外凸轮廓而外圈是内凹轮廓.对比图5(a)和5(b)发现，随着歪斜角的增加，当工况相同时，滚子与内圈的接触应力由大于滚子与外圈的接触应力逐渐转变为小于滚子与外圈的接触应力.

3.3 滚子长径比对接触应力分布的影响

分析滚子长径比对接触应力的影响时，先保持滚子半径恒定，令滚子在单位长度上受到的接触压力不变，那么接触区域半宽b 也保持不变，因而具有相同的最大赫兹接触压力. 图6 为分别计算滚子长径比L/D 为1、3、5，且歪斜角θ =0.5°，偏斜角Ψ =0. 02°时，滚子与滚道的接触应力分布.由图6(a)可知:随着长径比的增加，由于偏斜角使得滚子与内圈滚道的重载端的接触应力增加，“偏斜效应”增强;歪斜角的存在使得滚子与滚道的相对接触长度减小，“歪斜效应”增强.图6(b)表示不同长径比滚子与外圈的接触应力分布，可得到与图6(a)相似的结论. 对比图6(a)和6(b)发现，在歪斜角较小时，滚子与内圈的接触应力大于滚子与外圈的接触应力.

图5 不同歪斜角滚子与滚道的接触应力分布Fig.5 Contact stress distribution at different skewed angles between roller and raceway

图6 不同长径比滚子与滚道的接触应力分布Fig.6 Contact stress distribution at different length-to-diameter ratios between roller and raceway

4 偏歪斜工况下滚子的凸度设计

滚子轴承一般同时运行于偏斜和歪斜工况下，针对其工况的滚子凸度设计可以通过减小滚子与滚道间的最大接触应力，以提高滚子轴承的抗偏、歪斜能力. 根据滚子凸度加工的现有工艺，通过增加滚子凸度量以减小滚子接触副的最大接触应力是一种便于加工且较为经济的方法［14-16］.由分析可知，当滚子存在较小的偏、歪斜角时，滚子与内圈滚道的接触应力要大于滚子与外圈滚道的接触应力，因此可通过分析不同凸度量滚子与内圈滚道的接触应力分布，提出一种应用于滚子在偏、歪斜工况下的凸度设计方法. 选取的滚子轴承参数为:滚子半径R =5 mm，滚子长度L =12.5 mm，内圈滚道半径RI=10 mm，滚子母线为对数母线.计算外载荷F =5 kN、歪斜角θ=1°偏斜角Ψ =0.02°时不同滚子凸度量(Td)下的滚子与内圈滚道的接触应力分布，见图7.

图7 不同滚子凸度量下接触区域中心线处应力分布Fig.7 Stress distribution on the center line of contact area with different roller crowning values

图7 中，不同的滚子凸度量对应不同修形载荷.其中，凸度量指距离对数滚子端部2 mm 处滚子的半径跌落量，0.68 μm 代表滚子在正载时的最佳凸度量.由图7 可知，在偏歪斜工况下，逐渐增加滚子的凸度量，滚子与滚道的最大接触应力先变小后变大，且接触区域逐渐缩小. 由于当滚子凸度量逐渐增加时，滚子总存在一个最优的凸度值使得滚子与滚道的接触应力分布最为合理，材料利用率最高，由此可以对偏歪斜工况下滚子进行凸度设计.

对所编程序进行修改，设定一修形载荷增量ΔF，通过迭代计算出滚子在不同修形载荷下的接触应力，直到求得的接触应力分布曲线最为理想.由分析可得，滚子凸度量的大小与滚子所受外载荷成正比，而与滚子有效长度成反比关系. 因此，取修形载荷增量为

根据以上方法进行编程计算，可得出滚子最佳凸度曲线和修形后的中心应力分布，如图8 所示.

图8 滚子修形凸度曲线和修形前后中心应力分布曲线Fig.8 Modification crown curves of roller and stress distribution curves of roller center before and after modification

由图8 可知，当滚子歪斜角θ =1°、偏斜角Ψ =0.02°时，对应的最佳凸度量为1.13 μm，滚子与滚道的最大接触应力为2.23 GPa，而滚子在正载工况下的凸度量为0.68 μm，此凸度量下的最大接触应力为2.77 GPa，修形后的滚子与滚道的最大接触应力减小了19.5%，可见滚子的凸度修形可以提高轴承的抗偏、歪斜能力，从而在一定程度上提高滚子轴承的疲劳寿命.

5 结论

文中根据静弹性接触理论，建立了滚子相对轴承滚道在偏、歪斜工况下的力学模型，并利用FFT和共轭梯度法实现了对偏、歪斜工况下滚子接触副的接触应力问题的分析.通过比较，证明该算法具有较快的计算速度和收敛性.研究还发现:滚子歪斜角一定时，随着滚子偏斜角增大，滚子与内外圈滚道接触产生的“偏斜效应”逐渐明显，并且“偏斜效应”能明显加大“边缘效应”;滚子偏斜角一定时，随着歪斜角的增加，“歪斜效应”逐渐明显，且滚子与内圈接触时“歪斜效应”更加明显，且一定的偏斜角能减弱滚子的“边缘效应”;而滚子与外圈滚道由于偏斜角的存在产生更加严重的“偏斜效应”和“歪斜效应”;随着滚子长径比的增加，由于偏斜角和歪斜角的存在，使得滚子与滚道接触的“偏斜效应”和“歪斜效应”都明显增强.根据滚子轴承实际运行过程中的偏、歪斜工况，文中提出了一种综合考虑偏斜和歪斜工况下的滚子凸度设计方法，可提高滚子轴承的抗偏、歪斜能力，在一定程度上改善滚子轴承的疲劳寿命.当然，目前的研究只限于单滚子接触副模型，而整个滚动轴承的偏歪斜工况需要进一步研究，同时可以考虑润滑因素的影响使分析更加符合实际工况.

［1］Schaude B. Optimale wlzk rperpro filierung von zylinder rollen lagern［D］.Karlsruhe:Karlsruhe University，1978:36-44.

［2］Kannel J W，Hatnett M J. Theoretical and experimental evaluation of edge stresses under sever edge loads［J］.ASLE Transactions，1981，26(1):25-30.

［3］毛月新，沈雪瑾，陈晓阳，等. 偏载滚子轴承的接触应力分布计算及其滚子凸度设计［J］. 中国机械工程，2009，20(16):1918-1922.Mao Yue-xin，Shen Xue-jin，Chen Xiao-yang，et al. Study on contact stress distribution of roller bearings and roller profile design under misaligned loads［J］.China Mechanical Engineering，2009，20(16):1918-1922.

［4］孙殿超，陈晓阳，王志坚，等.利用FFT 和共轭梯度法求解偏载滚子副的接触应力［J］. 机械设计与制造，2014，11(11):43-46.Sun Dian-chao，Chen Xiao-yang，Wang Zhi-jian，et al.Numerical solution of the contact stress for load-offsetting roller pairs based on fast Fourier transform and conjugate method［J］. Machinery Design ＆ Manufacture，2014，11(11):43-46.

［5］Ye Zhen-Huan，Wang Li-qin，Gu Le，et al.Effects of tilted misalignment on loading characteristics of cylindrical roller bearings［J］.Mechanism and Machine Theory，2013，69:153-167.

［6］Tong Bao-hong，Liu Ying，Cheng Xin-ming，et al. Finiteelement analysis on dynamic contact characteristics of cylindrical roller bearing considering misalignment with journal［J］.Transactions of the Chinese Society of Agricultural Engineering，2012，28(13):42-47.

［7］Harris T A，Kotzalas M N，Yu W K. On The causes and effects of roller skewing in cylindrical roller bearing［J］.Tribology Transactions，1998，41(4):572-578.

［8］Yang Y，Danyluk S，Hoeprich M. A study on rolling element skew measurement in a tapered roller bearing with a specialized capacitance probe［J］. Journal of Tribology，2000，122(3):534-538.

［9］李伟伟，陈晓阳，沈雪瑾，等. 歪斜状况下滚子轴承的接触应力求解及分析［J］. 中国机械工程，2011，22(17):2034-2039.Li Wei-wei，Chen Xiao-yang，Shen Xue-jin，et al. Contact pressure calculation and analysis of roller bearings in skewing condition［J］.China Mechanical Engineering，2011，22(17):2034-2039.

［10］Love A E H. A treatise on the mathematical theory of elasticity［D］.Cambridge:Cambridge University Press，1927:21-27.

［11］王文中，胡元中，王慧. 基于快速傅里叶变换和共轭梯度法求解干接触问题［J］.机械工程学报，2006，42(7):14-18.Wang Wen-zhong，Hu Yuan-zhong，Wang Hui. Numerical solution of dry contact problem based on fast Fourier transform and conjugate gradient method ［J］. Chinese Journal of Mechanical Engineering，2006，42(7):14-18.

［12］Liu S，Wang Q，Liu G. A versatile method of discrete convolution and FFT (DC-FFT)for contact analyses［J］.Wear，2000，243(1):101-111.

［13］Wang W Z，Wang H，Liu Y C.A comparative study of the methods for calculation of surface elastic deformation［J］.Journal of Engineering Tribology，2003，217(2):145-153.

［14］陈晓阳，马家驹.对数凸型滚针凸度量设计的新方法［J］.摩擦学学报，1993，13(1):25-32.Chen Xiao-yang，Ma Jia-ju. A new numerical evaluation of crown drop for a needle roller with Lundberg’s profile［J］.Tribology，1993，13(1):25-32.

［15］王志伟，苗利军，孟玲琴. 圆锥滚子轴承滚子凸度的优化设计［J］.机械设计与制造，2014，6(6):55-61.Wang Zhi-wei，Miao Li-jun，Meng Ling-qin.Optimal design of roller convexity of tapered roller bearing［J］.Machinery Design ＆ Manufacture，2014，6(6):55-61.

［16］Fujiwara H，Kobayashi T，Kawase T，et al. Optimized logarithmic roller crowning design of cylindrical roller bearings and its experimental demonstration［J］.Tribology Transaction，2010，53:909-916.