歪斜狀況下滾子軸承的接觸應力求解與分析

李偉偉 陳曉陽 沈雪瑾 張小玲

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.