彈性接觸問題的快速算法及其在圓錐滾子接觸應力分析中的應用

米俊，夏伯乾

(鄭州大學 機械工程學院，鄭州 450001)

接觸問題是滾動軸承壽命設計中的基本問題，對滾子軸承而言，由于滾子和滾道之間的接觸是有限長線接觸，而且滾子端部有應力集中，因此經典的求解滾子軸承接觸問題的Hertz線接觸理論就無法得到理想的結果。對非Hertz接觸問題，目前廣泛地采用數值解法。數值分析的手段有2類：一是使用有限元軟件，如ANSYS，Nastran，I-DEAS等；二是使用經典的數值分析方法，即在二維區域內對接觸問題的積分方程進行數值求解。用有限元軟件分析滾動軸承的接觸問題，過程復雜、占用計算資源較多，很不經濟，所以現實中大量關于滾動軸承接觸問題的分析與研究仍然廣泛應用經典的數值分析方法。

圓錐滾子由于幾何形狀的復雜性，關于其接觸應力計算與分析方面的研究文獻還比較少。文中將一種求解彈性表面接觸問題的快速直接迭代算法[1]應用于圓錐滾子接觸應力的求解，得到了不同載荷下圓錐滾子的中心接觸應力分布，為圓錐滾子軸承的優化設計提供參考。

1 彈性接觸問題的基本方程

彈性體A和B在外加載荷F作用下相互接觸，如圖1所示。在這里先作如下假設：(1)接觸體之間的摩擦可以忽略不計；(2)接觸體幾何尺寸遠大于接觸區域的幾何尺寸；(3)接觸區域視為平面。則有

圖1 兩個彈性體接觸示意圖

ω1+ω2+z1+z2=δ(接觸區域以內)，

(1)

ω1+ω2+z1+z2>δ(接觸區域以外)，

(2)

式中：δ為載荷作用下兩彈性體之間產生的彈性趨近量；ω1，ω2分別為彈性體A和B在接觸點的彈性變形量。其中[2]，

式中：p為接觸應力；υ為泊松比；E為彈性模量；s為接觸應力作用面積；x′,y′為積分變量；i=1,2。

將(3)式代入(1)式中，可得變形協調方程為

(4)

接觸應力還應滿足平衡方程

(5)

在接觸區s內p(x,y)≥0；在接觸區s外p(x,y)=0。

(4)和(5)式構成了彈性接觸問題的基本方程。

2 接觸問題的快速直接迭代算法

在求解應力之前，要將求解區域分為若干個單元，假定每個單元結點上的應力為定值pij，則變形協調方程(4)離散后變為

(6)

(5)式離散后的表達式為

(7)

式中：α′為矩形單元面積，α′=4ab;pij≥0。

(6)式是m×n維的大型方程組，而且已經被證明是1個病態方程組[3-4]，求解起來很困難。因此需要設法將 (6)式降階，并克服其病態。具體方法是將整個區域分為若干部分分別求解，在離散后的求解區域上，利用Gauss-Seidel迭代法對方程組進行求解，構造迭代格式如下

(8)

Pj=(pj1,pj2,…,pjn),

通過Gauss-Seidel迭代法求解以上方程組，雖然降低了方程組的維數，也解決了方程組的病態問題，但(8)式中每個子方程組的系數矩陣都是滿元矩陣，當接觸區域的網格劃分比較稠密時，需要耗費較多計算時間。

(9)

這樣一來， 方程組(9)中每個等式展開后皆為方程組,且此方程組的系數矩陣都由滿元矩陣變成了帶寬為2β+1的稀疏帶狀矩陣，其中β為矩陣的半帶寬，變換后的帶狀矩陣與前面的滿元矩陣相比，運算速度要快得多，節省了大量計算時間[1]。并不是半帶寬β越大，計算精度越高，當半帶寬β取5～7時，即當迭代矩陣密度為滿陣的15%～20%時，計算精度最高。半帶寬β取值過大或者說迭代矩陣密度過高，反而不利于提高計算精度。則(9)式和離散后的 (7) 式便構成求解任意彈性表面接觸問題的快速直接迭代算法。對于這種直接迭代算法的求解過程，其程序設計框圖如圖2所示。

圖2 程序設計框圖

3 圓錐滾子接觸問題的幾何模型

圖3給出了圓錐滾子與平面接觸的幾何模型。圓錐滾子被垂直于y軸的平面所截的面是大小不同的橢圓，這些截取的橢圓在接觸點附近的曲率是不相同的,因此這些橢圓的接觸變形和接觸應力也就不同。

根據圖3的幾何關系，距離圓錐滾子小端Li處，圓錐體橫截面的直徑di和垂直于y軸的橢圓截面的長半軸ai為[6]

圖3 圓錐滾子與平面接觸的幾何模型

di=(1-Li/L)d1+d2Li/L，

(10)

ai=0.5di(cosθ+sinθtan 2θ)，

(11)

式中：d1為圓錐滾子小端直徑；d2為圓錐滾子大端直徑；L為滾子素線有限長度；θ為圓錐滾子的半錐角。

因此,橢圓短半軸可以近似地表示為

bi=0.5di(1+sin2θ)。

(12)

通過計算可知，當θ=10°時，(12)式的相對誤差小于0.2%，而實際的圓錐滾子半錐角在2°左右，因此采用(12)式可以滿足計算精度要求。故橢圓截面在接觸點附近的曲率半徑為[7]

(13)

4 數值算例

算例中圓錐滾子的小端和大端直徑分別為d1=10.0 mm，d2=12.0 mm，滾子的有限長度L=15.0 mm，半錐角θ=3°，圓錐滾子和平面的彈性模量及泊松比分別為E1=E2=210 GPa，υ1=υ2=0.3。應力迭代矩陣半帶寬β=6。

圖4是外加載荷為3 000 N，網格節點數為35×35時圓錐滾子與平面接觸的三維應力分布圖。從圖4中可以看出圓錐滾子兩端的應力集中效應，而且小端的接觸應力要比大端的稍微大一些，其他部分應力分布比較均勻。

圖4 圓錐滾子的三維應力分布

圖5為圓錐滾子承受5種不同載荷時，滾子中心平面的接觸應力分布情況。從圖5中可以看出，同一曲線上滾子中心應力值右端要比左端大一些，即圓錐滾子小端的中心接觸應力比大端的中心接觸應力要大一些。這是因為圓錐滾子的小端曲率比大端大，應力集中明顯一些。在其他的部分，接觸應力分布比較均勻。

圖5 不同載荷下中心平面接觸應力分布曲線

5 結論

(1)給出了圓錐滾子接觸問題的簡化計算模型，將一種求解任意彈性表面接觸問題的快速算法應用到圓錐滾子接觸應力的求解中，并得到了滿意的結果。

(2)計算出了不同載荷下圓錐滾子中心平面的接觸應力分布。通過分析可知，在一定外載荷作用下，圓錐滾子兩端存在應力集中，而且小端的接觸應力比大端的接觸應力要大，在沿滾子素線的其他部分，接觸應力分布則比較均勻。