彈性接觸問(wèn)題的快速算法及其在圓錐滾子接觸應(yīng)力分析中的應(yīng)用

米俊，夏伯乾

(鄭州大學(xué) 機(jī)械工程學(xué)院，鄭州 450001)

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

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

1 彈性接觸問(wèn)題的基本方程

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

圖1 兩個(gè)彈性體接觸示意圖

ω1+ω2+z1+z2=δ(接觸區(qū)域以?xún)?nèi))，

(1)

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

(2)

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

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

將(3)式代入(1)式中，可得變形協(xié)調(diào)方程為

(4)

接觸應(yīng)力還應(yīng)滿足平衡方程

(5)

在接觸區(qū)s內(nèi)p(x,y)≥0；在接觸區(qū)s外p(x,y)=0。

(4)和(5)式構(gòu)成了彈性接觸問(wèn)題的基本方程。

2 接觸問(wèn)題的快速直接迭代算法

在求解應(yīng)力之前，要將求解區(qū)域分為若干個(gè)單元，假定每個(gè)單元結(jié)點(diǎn)上的應(yīng)力為定值pij，則變形協(xié)調(diào)方程(4)離散后變?yōu)?/p>

(6)

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

(7)

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

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

(8)

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

通過(guò)Gauss-Seidel迭代法求解以上方程組，雖然降低了方程組的維數(shù)，也解決了方程組的病態(tài)問(wèn)題，但(8)式中每個(gè)子方程組的系數(shù)矩陣都是滿元矩陣，當(dāng)接觸區(qū)域的網(wǎng)格劃分比較稠密時(shí)，需要耗費(fèi)較多計(jì)算時(shí)間。

(9)

這樣一來(lái)， 方程組(9)中每個(gè)等式展開(kāi)后皆為方程組,且此方程組的系數(shù)矩陣都由滿元矩陣變成了帶寬為2β+1的稀疏帶狀矩陣，其中β為矩陣的半帶寬，變換后的帶狀矩陣與前面的滿元矩陣相比，運(yùn)算速度要快得多，節(jié)省了大量計(jì)算時(shí)間[1]。并不是半帶寬β越大，計(jì)算精度越高，當(dāng)半帶寬β取5～7時(shí)，即當(dāng)?shù)仃嚸芏葹闈M陣的15%～20%時(shí)，計(jì)算精度最高。半帶寬β取值過(guò)大或者說(shuō)迭代矩陣密度過(guò)高，反而不利于提高計(jì)算精度。則(9)式和離散后的 (7) 式便構(gòu)成求解任意彈性表面接觸問(wèn)題的快速直接迭代算法。對(duì)于這種直接迭代算法的求解過(guò)程，其程序設(shè)計(jì)框圖如圖2所示。

圖2 程序設(shè)計(jì)框圖

3 圓錐滾子接觸問(wèn)題的幾何模型

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

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

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

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

(10)

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

(11)

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

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

bi=0.5di(1+sin2θ)。

(12)

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

(13)

4 數(shù)值算例

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

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

圖4 圓錐滾子的三維應(yīng)力分布

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

圖5 不同載荷下中心平面接觸應(yīng)力分布曲線

5 結(jié)論

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

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