陳曉陽 周越 沈雪瑾 王志堅 張小玲 馬純青

(1.上海大學 機電工程與自動化學院，上海200072;2.蘇州軸承廠有限公司，江蘇 蘇州215000;3.上海和錦滾子科技有限公司，上海201611)

徑向尺寸小、結構簡單、承載能力大等優(yōu)點使得滾子軸承成為機械工業(yè)領域的重要零部件. 但由于安裝誤差、熱變形、幾何形狀缺陷等原因，滾子軸承一般在偏、歪斜等工況下運行，使得軸承過早失效.Schaude［1］最先利用切片法對偏斜工況下滾子軸承接觸應力展開研究，通過數(shù)值解法求解出了偏斜工況下滾子軸承接觸應力分布. Kannel 等［2］通過切片法求解出了滾子接觸副的接觸應力分布，并對其結果進行了實驗驗證.毛月新等［3］利用影響系數(shù)法分析了偏斜工況下對數(shù)滾子與滾道的接觸應力分布，并求解了滾子在偏斜工況下的最佳凸度量. 孫殿超等［4］利用快速傅里葉變換(FFT)和共軛梯度法求解了偏斜滾子接觸副的接觸應力，并闡述了該算法具有較快的計算速度和收斂性. Ye 等［5-6］利用有限元法分析了偏斜工況下滾子的接觸問題. Harris 等［7］分析了滾子軸承中滾子偏斜和歪斜產(chǎn)生的原因及其影響.Yang 等［8］利用探針測出了圓錐滾子軸承在不同轉速、載荷及潤滑油下的歪斜角. 李偉偉等［9］采用影響系數(shù)法分析了歪斜工況滾子接觸副的接觸應力分布.但滾子軸承在偏斜和歪斜共同作用工況下的理想凸度滾子與內外圈滾道的接觸問題的研究鮮有報道，文中將采用FFT 和共軛梯度法求解上述問題，并討論對應工況下的滾子最佳凸度量.

1 偏歪斜工況滾子與滾道的接觸分析

1.1 偏歪斜工況下滾子與滾道的力學模型

在實際工況下，由于軸的變形、軸承中心不對中以及軸承的徑向游隙等原因，使得滾子接觸副處于偏斜狀態(tài)，如圖1(a)所示.偏斜的存在使?jié)L子在一端被擠壓，壓向引導擋邊產(chǎn)生摩擦力，同時由于滾子與滾道的軸向間隙等原因使得滾子處于歪斜狀態(tài)下工作，如圖1(b)所示.根據(jù)靜彈性接觸理論，半徑為R、長度為L 的滾子與滾道產(chǎn)生Ψ 度偏斜和θ 度的歪斜，M、N 分別是滾子和半徑為RI的內圈外表面上的點，其坐標表示為(x，y，z1)、(x，y，z2).在外載荷F 作用下，起始接觸點O 將延展為實際接觸區(qū)域Ω，接觸區(qū)域上的應力分布用函數(shù)P(x，y)表示，如圖2 所示.

圖1 滾子與滾道的接觸力學模型Fig.1 Mechanics model describing the contact of roller and raceway

圖2 滾子與滾道受載后的接觸區(qū)域應力分布示意圖Fig.2 Schematic diagram of stress distribution of roller and raceway in contact area after loading

根據(jù)Boussinesq 半無限體的力和變形解，在實際接觸區(qū)域Ω 內可以求解出半無限彈性體接觸問題的主導積分方程——Fredholm 方程［10］表示.為方便起見，給出其離散形式:

式中，Gij為接觸間隙矩陣，Uij為表面彈性變形矩陣，Hij為接觸體幾何形狀矩陣，Δ 為彈性趨近量.

偏斜工況下，彈性趨近量與接觸副接觸中心處的彈性趨近量δ0有以下關系:

歪斜工況下，接觸體幾何形狀矩陣可表示為

其中，Sij為滾子與平面P 之間的幾何間距，

Td為滾子凸度函數(shù).

Tij為內滾道與平面P 之間的幾何間距，

其中，

由此，即可求出彈性變形計算式及載荷和力矩平衡方程.若需求滾子與外滾道的接觸體幾何形狀矩陣，只要改變Tij函數(shù)即可，這里不再贅述.

1.2 數(shù)值解法

1.2.1 FFT 求解彈性變形

在數(shù)值計算過程中接觸副彈性變形的求解占據(jù)著主要的工作量.由彈性變形計算式是一卷積求和運算，可以利用FFT 進行快速求解，計算速度優(yōu)于多重網(wǎng)格法［11］.具體計算可參閱文獻［11-13］.

1.2.2 共軛梯度法求解接觸壓力

通過變分原理將接觸問題轉變?yōu)榍蠼鈽O值問題后，就可以通過共軛梯度法進行數(shù)值計算.該算法流程如圖3 所示，具體算法的實現(xiàn)過程可參閱文獻［4］.圖中，Q、Q′分別為初始載荷和計算載荷，pij為壓力矩陣元素，tij為迭代方向為修正量系數(shù)，ax、ay分別為網(wǎng)格的長度和寬度，ε0為迭代精度，ε1、ε2分別為壓力和矩陣的迭代誤差.

圖3 計算流程圖Fig.3 Flowchart for model solution

2 程序正確性驗證和計算時間比較

為便于比較，文獻［3，9］中的計算是文中的特例，即歪斜角為0°、偏斜角為0.02°和偏斜角為0°、歪斜角為2°. 在所用計算機為Intel Pentium Dual CPU E2160、計算環(huán)境為Visual Fortran 6.5 且網(wǎng)格數(shù)均為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］計算結果的滾子與滾道接觸區(qū)域半寬的相對誤差都小于1.0%，最大接觸應力的相對誤差都小于2.0%，說明兩種方法計算的接觸應力和接觸區(qū)域保持一致，證明了算法程序的準確性.特別是在計算時間方面，文中算法耗時大大減少，運算速度約是文獻［3，9］計算的200 倍，從而證明了該方法具有高效的計算速度和收斂性.

表2 文中計算結果與文獻［9］計算結果對比1)Table 2 Comparison of calculated results obtained in the paper and from literature［9］

3 各因素對接觸應力分布的影響

滾子軸承參數(shù):滾子副綜合彈性模量E=226GPa，滾子半徑R =5 mm，滾子長度L =25 mm，內圈滾道半徑RI=10 mm，外圈滾道半徑RO=20 mm，滾子母線為在正載工況下的最佳對數(shù)母線.

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 可知:當滾子的歪斜角一定時，隨著偏斜角的增大，滾子與內外圈滾道的“偏斜效應”逐漸明顯，重載端的接觸應力明顯大于輕載端的接觸應力，并且“偏斜效應”能明顯加大“邊緣效應”.當偏斜角增加到一定值時，輕載端的滾子與滾道將會脫離接觸，其接觸區(qū)域長度隨著偏斜角的增加而減小.對比圖4(a)和4(b)發(fā)現(xiàn)，當工況相同時，滾子與內圈的接觸區(qū)域小于與外圈的接觸區(qū)域，導致滾子與內圈的接觸應力大于滾子與外圈的接觸應力.

3.2 不同歪斜角對滾子接觸應力的影響

當外載荷F=10 kN，偏斜角Ψ =0.02°時，計算滾子歪斜角θ= 0°，1°，2°作用下滾子與滾道間的接觸應力，如圖5 所示.

圖5(a)為滾子偏斜角一定、歪斜角增加時滾子與內圈的接觸應力分布. 研究發(fā)現(xiàn):滾子發(fā)生歪斜時，滾子與內圈接觸將會產(chǎn)生“歪斜效應”，并且當歪斜角逐漸變大時，滾子與內滾道的接觸應力在接觸副的中部變大，兩端變小.“歪斜效應”逐漸明顯，且一定的偏斜角能減弱滾子的“邊緣效應”，滾子與內圈接觸區(qū)域的長度方向呈外凸的輪廓曲線. 圖5(b)為滾子偏斜角一定、歪斜角增加時的滾子與外滾道的接觸應力分布. 可知:隨著歪斜角的增大，滾子與外圈滾道的接觸應力在接觸副的兩端變大，中部變小.滾子與外圈由于偏斜角會產(chǎn)生更加嚴重的“偏斜效應”和“歪斜效應”，并且滾子與外圈接觸區(qū)域的長度方向呈外凹的輪廓曲線，這主要是由內外滾道的結構決定的，即內圈是外凸輪廓而外圈是內凹輪廓.對比圖5(a)和5(b)發(fā)現(xiàn)，隨著歪斜角的增加，當工況相同時，滾子與內圈的接觸應力由大于滾子與外圈的接觸應力逐漸轉變?yōu)樾∮跐L子與外圈的接觸應力.

3.3 滾子長徑比對接觸應力分布的影響

分析滾子長徑比對接觸應力的影響時，先保持滾子半徑恒定，令滾子在單位長度上受到的接觸壓力不變，那么接觸區(qū)域半寬b 也保持不變，因而具有相同的最大赫茲接觸壓力. 圖6 為分別計算滾子長徑比L/D 為1、3、5，且歪斜角θ =0.5°，偏斜角Ψ =0. 02°時，滾子與滾道的接觸應力分布.由圖6(a)可知:隨著長徑比的增加，由于偏斜角使得滾子與內圈滾道的重載端的接觸應力增加，“偏斜效應”增強;歪斜角的存在使得滾子與滾道的相對接觸長度減小，“歪斜效應”增強.圖6(b)表示不同長徑比滾子與外圈的接觸應力分布，可得到與圖6(a)相似的結論. 對比圖6(a)和6(b)發(fā)現(xiàn)，在歪斜角較小時，滾子與內圈的接觸應力大于滾子與外圈的接觸應力.

圖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 偏歪斜工況下滾子的凸度設計

滾子軸承一般同時運行于偏斜和歪斜工況下，針對其工況的滾子凸度設計可以通過減小滾子與滾道間的最大接觸應力，以提高滾子軸承的抗偏、歪斜能力. 根據(jù)滾子凸度加工的現(xiàn)有工藝，通過增加滾子凸度量以減小滾子接觸副的最大接觸應力是一種便于加工且較為經(jīng)濟的方法［14-16］.由分析可知，當滾子存在較小的偏、歪斜角時，滾子與內圈滾道的接觸應力要大于滾子與外圈滾道的接觸應力，因此可通過分析不同凸度量滾子與內圈滾道的接觸應力分布，提出一種應用于滾子在偏、歪斜工況下的凸度設計方法. 選取的滾子軸承參數(shù)為:滾子半徑R =5 mm，滾子長度L =12.5 mm，內圈滾道半徑RI=10 mm，滾子母線為對數(shù)母線.計算外載荷F =5 kN、歪斜角θ=1°偏斜角Ψ =0.02°時不同滾子凸度量(Td)下的滾子與內圈滾道的接觸應力分布，見圖7.

圖7 不同滾子凸度量下接觸區(qū)域中心線處應力分布Fig.7 Stress distribution on the center line of contact area with different roller crowning values

圖7 中，不同的滾子凸度量對應不同修形載荷.其中，凸度量指距離對數(shù)滾子端部2 mm 處滾子的半徑跌落量，0.68 μm 代表滾子在正載時的最佳凸度量.由圖7 可知，在偏歪斜工況下，逐漸增加滾子的凸度量，滾子與滾道的最大接觸應力先變小后變大，且接觸區(qū)域逐漸縮小. 由于當滾子凸度量逐漸增加時，滾子總存在一個最優(yōu)的凸度值使得滾子與滾道的接觸應力分布最為合理，材料利用率最高，由此可以對偏歪斜工況下滾子進行凸度設計.

對所編程序進行修改，設定一修形載荷增量ΔF，通過迭代計算出滾子在不同修形載荷下的接觸應力，直到求得的接觸應力分布曲線最為理想.由分析可得，滾子凸度量的大小與滾子所受外載荷成正比，而與滾子有效長度成反比關系. 因此，取修形載荷增量為

根據(jù)以上方法進行編程計算，可得出滾子最佳凸度曲線和修形后的中心應力分布，如圖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 結論

文中根據(jù)靜彈性接觸理論，建立了滾子相對軸承滾道在偏、歪斜工況下的力學模型，并利用FFT和共軛梯度法實現(xiàn)了對偏、歪斜工況下滾子接觸副的接觸應力問題的分析.通過比較，證明該算法具有較快的計算速度和收斂性.研究還發(fā)現(xiàn):滾子歪斜角一定時，隨著滾子偏斜角增大，滾子與內外圈滾道接觸產(chǎn)生的“偏斜效應”逐漸明顯，并且“偏斜效應”能明顯加大“邊緣效應”;滾子偏斜角一定時，隨著歪斜角的增加，“歪斜效應”逐漸明顯，且滾子與內圈接觸時“歪斜效應”更加明顯，且一定的偏斜角能減弱滾子的“邊緣效應”;而滾子與外圈滾道由于偏斜角的存在產(chǎn)生更加嚴重的“偏斜效應”和“歪斜效應”;隨著滾子長徑比的增加，由于偏斜角和歪斜角的存在，使得滾子與滾道接觸的“偏斜效應”和“歪斜效應”都明顯增強.根據(jù)滾子軸承實際運行過程中的偏、歪斜工況，文中提出了一種綜合考慮偏斜和歪斜工況下的滾子凸度設計方法，可提高滾子軸承的抗偏、歪斜能力，在一定程度上改善滾子軸承的疲勞壽命.當然，目前的研究只限于單滾子接觸副模型，而整個滾動軸承的偏歪斜工況需要進一步研究，同時可以考慮潤滑因素的影響使分析更加符合實際工況.

［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］陳曉陽，馬家駒.對數(shù)凸型滾針凸度量設計的新方法［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］王志偉，苗利軍，孟玲琴. 圓錐滾子軸承滾子凸度的優(yōu)化設計［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.