圓柱滾子軸承數(shù)值求解方法的研究與分析

李猛帥，鄭甲紅

（陜西科技大學(xué)，陜西 西安 710021）

1 確定相關(guān)參數(shù)

圓柱滾子軸承的接觸問題屬于有限長線接觸問題，此類問題是在基于點(diǎn)接觸問題的基礎(chǔ)上展開討論的，因此點(diǎn)接觸問題的一些結(jié)論對線接觸是有用的。包括以下參數(shù)：主曲率半徑R、變形前后滾子表面對應(yīng)點(diǎn)與起始接觸點(diǎn)之間的距離Z、接觸應(yīng)力與外載荷之間的平衡方程。

主曲率半徑是根據(jù)圓柱滾子軸承的實(shí)際滾子狀態(tài)來確定的，變形前后滾子表面對應(yīng)點(diǎn)與起始接觸點(diǎn)之間的距離與點(diǎn)接觸不同，它在計(jì)算程序中作為初始參數(shù)輸入，而且需要將其在計(jì)算程序中轉(zhuǎn)化成向量的形式，它的確定方式如下。

采用文獻(xiàn)[1]中的式1.21—式1.22b得出：

式（2）（3）中：R11、R12、R21、R22分別為2個(gè)接觸物體在一點(diǎn)上的主曲率半徑。

對于圓柱滾子軸承而言，接觸角α=0°，R12、R21、R22都為無窮大，代入式（2）與式（3）并將2式相加可得：

這樣就確定了Z的表達(dá)式，在MATLAB軟件中輸入時(shí)，需要將其以向量的形式輸入，而且它的個(gè)數(shù)也要和2種數(shù)值求解方法中的網(wǎng)格數(shù)一一對應(yīng)，因?yàn)檫@樣才能滿足矩陣之間的基本運(yùn)算，從而進(jìn)一步求解基本平衡方程。

線接觸問題的基本平衡方程與點(diǎn)接觸問題的基本方程在原理上是類似的，可理解為是將點(diǎn)接觸問題的基本方程轉(zhuǎn)化成求和的形式。

2 非Hertz接觸問題的求解方法

有限長線接觸問題中短圓柱滾子的接觸狀態(tài)一般是在非理想狀態(tài)下的，實(shí)際接觸狀態(tài)比較復(fù)雜。一般在求解時(shí)，是將所有可能接觸區(qū)域劃分成多個(gè)矩形小單元格，而且要注意所求解區(qū)域要稍微大于實(shí)際接觸區(qū)域，這樣也符合實(shí)際工程要求。

劃分矩形單元如圖1所示，將可能接觸的區(qū)域劃分成多個(gè)矩形小單元，其中任意一個(gè)單元j的半長和半寬分別為aj和bj，假設(shè)單元上的接觸應(yīng)力為pj[1]。

圖1 劃分矩形單元

則單元j上的接觸應(yīng)力在任意一個(gè)單元i的中心上產(chǎn)生的位移為：

式（6）中：Dij為柔度系數(shù)。

柔度系數(shù)Dij通過MATLAB軟件進(jìn)行計(jì)算，而且它的下標(biāo)在MATLAB中代表著任意一個(gè)矩陣，也就是單元矩陣劃分多少個(gè)，那么對應(yīng)的柔度系數(shù)矩陣在MATLAB中就有多少個(gè)。那么線接觸的基本方程就可以表示為：

式（7）（8）中：m、n為單元格數(shù)目，它們與柔度系數(shù)矩陣下標(biāo)表示的意義是一樣的，在MATLAB中輸入它們時(shí)，將二者設(shè)定為相等數(shù)目；E′為彈性模量；δ為彈性趨近量。

需要注意的是，在計(jì)算時(shí)要滿足pj≥0，因?yàn)楫?dāng)pj＜0時(shí)，代表此單元格已經(jīng)超出實(shí)際接觸區(qū)域，沒有受到外載荷的作用，在下一次程序迭代時(shí)就不再進(jìn)行計(jì)算。

整個(gè)關(guān)鍵參數(shù)在MATLAB程序中運(yùn)算時(shí)，按以下流程進(jìn)行：首先輸入外載荷Q以及物理幾何參數(shù)，計(jì)算向量Zi，并且將柔度系數(shù)Dij計(jì)算出來后轉(zhuǎn)化成矩陣形式；程序中δ雖然作為未知量，但是它是通過基本方程和誤差修正最終確定的，因此需要賦初值；其次通過解線接觸的基本方程求出接觸應(yīng)力pj，并判斷pj是否小于0，若pj＞0則進(jìn)行下一步，如果pj＜0，則將pj的值賦為0，循環(huán)到解式（8）；在解方程時(shí)，pj=0表示在接下來的計(jì)算時(shí)，將pj＜0的單元格跳過，直接進(jìn)行下一單元的計(jì)算，以此類推。最后一步的判斷條件為取ε為0.001，若滿足條件，則結(jié)束程序；若不滿足，循環(huán)跳到求解式（8），重復(fù)進(jìn)行循環(huán)，直至將所有的接觸區(qū)域解算完畢。

3 Hertz接觸問題的求解方法

利用Hertz理論同樣也可以求解線接觸問題，式（7）和式（8）是接觸問題積分方程數(shù)值解的最簡單形式之一，它已被許多學(xué)者以不同的形式成功地應(yīng)用于分析Hertz問題和非Hertz問題[2-5]。然而基于對線接觸問題的理解，還可以對圓柱滾子軸承有限長線接觸問題的數(shù)值解進(jìn)行進(jìn)一步的簡化，這也就是對圓柱滾子有限長線接觸問題進(jìn)行處理的一維處理方法。

它的求解思想也是將所有可能的接觸區(qū)域進(jìn)行單元格劃分，不同的是在單元j內(nèi)假定接觸應(yīng)力沿滾子素線方向（y軸）為均勻分布，沿橫向（x軸）按Hertz分布[1]，條形單元?jiǎng)澐秩鐖D2所示。

圖2 條形單元?jiǎng)澐?/p>

圖2中p0j為單元中心處的最大接觸應(yīng)力，單元j上的應(yīng)力在其他任意單元格i的中心處產(chǎn)生的位移與非Hertz接觸一樣。

式（9）中：hj為單元格寬度的1/2；aj為半長。

此處柔度系數(shù)與Hertz不同，主要是在x軸方向上，將應(yīng)力分布看成是按Hertz分布的，那么就可以先對y′積分，然后對x′積分，就可以求解出柔度系數(shù)，再根據(jù)單元格劃分，將Dij在MATLAB軟件中輸入，并且同樣轉(zhuǎn)化為矩陣的形式。

Hertz線接觸的基本方程為：

它也是將點(diǎn)接觸的方程轉(zhuǎn)化為求和的形式，而且與非Hertz方程比較相似。

此套計(jì)算程序比非Hertz計(jì)算稍微復(fù)雜，主要流程為：首先輸入外載荷Q、向量Zi、網(wǎng)格劃分個(gè)數(shù)n及幾何材料參數(shù)，并且對δ賦初值，求解向量Si=[δ－Zi]，這一步是為了簡化式（11）的計(jì)算，并且對aj賦初值；其次計(jì)算柔度矩陣，求解式（11），并判斷解出的p0j是否大于0，大于0則進(jìn)行下一步，小于0時(shí)回到計(jì)算柔度矩陣這一步重新計(jì)算；如果滿足條件求解式（10），若不滿足則同樣跳回到計(jì)算柔度系數(shù)矩陣這一步，重新循環(huán)，目的是為了更精準(zhǔn)地自動(dòng)跳過非接觸區(qū)域的計(jì)算；最后，與非Hertz問題一樣，判斷根據(jù)需要調(diào)整δ的值，從程的值，從程序第一步重新計(jì)算。

4 實(shí)例計(jì)算

某二維轉(zhuǎn)臺(tái)的圓柱滾子軸承，已知滾子曲率半徑R11=10 mm，R11=R21=R22=∞，外載荷Q=2 000 N，單元格劃分個(gè)數(shù)n=18，彈性模量E=2.06×105N/mm2，滾子素線長L=30 mm。計(jì)算該圓柱滾子在外載荷作用下的接觸應(yīng)力及彈性趨近量。

根據(jù)編制好的2種計(jì)算程序，將所有初始條件輸入到MATLAB軟件中，通過對2種計(jì)算方法的仔細(xì)分析，把對結(jié)果有直接影響的參數(shù)一一列舉，并整合出曲線圖。網(wǎng)格劃分個(gè)數(shù)與單元格中心最大接觸應(yīng)力之間的關(guān)系如圖3所示。

圖3 網(wǎng)格劃分與最大接觸應(yīng)力之間的關(guān)系

由圖3可知，利用2種計(jì)算方法進(jìn)行單元格劃分時(shí)，不同單元格中心點(diǎn)處的最大接觸應(yīng)力是不同的，而且它們都隨著網(wǎng)格數(shù)的增加整體呈逐漸遞增的趨勢。兩者的數(shù)值比較接近，由于在非Hertz計(jì)算程序中，它的網(wǎng)格劃分?jǐn)?shù)不能超過27個(gè)，當(dāng)網(wǎng)格個(gè)數(shù)超過27個(gè)時(shí)，就會(huì)導(dǎo)致柔度系數(shù)矩陣的行列式為0，程序就無法進(jìn)行計(jì)算，因此具有一定的局限性，而Hertz計(jì)算程序中網(wǎng)格劃分?jǐn)?shù)可以取到更多。

圓柱滾子的主曲率半徑與彈性趨近量之間關(guān)系的曲線圖如圖4所示。由圖4可知，隨著圓柱滾子的主曲率半徑逐漸增大時(shí)，利用2種計(jì)算方法得出的彈性趨近量都呈逐漸遞減的趨勢，但是在非Hertz計(jì)算方法中曲線下降比較明顯，而Hertz計(jì)算中曲線下滑比較緩慢，彈性趨近量的隨主曲率半徑增大波動(dòng)很小，2種計(jì)算方法雖然有一定的誤差，但是誤差較小。一般采用Hertz計(jì)算方法進(jìn)行求解。

圖4 主曲率半徑與彈性趨近量之間的關(guān)系

網(wǎng)格劃分個(gè)數(shù)與彈性趨近量之間關(guān)系的曲線圖如圖5所示。由圖5可知，網(wǎng)格劃個(gè)分?jǐn)?shù)不同時(shí)，利用2種計(jì)算方法中得出的彈性趨近量也是有差異的，隨著矩形單元格個(gè)數(shù)增加，彈性趨近量都呈遞增趨勢，而且2種計(jì)算方法的結(jié)果比較接近，這里也需要注意在非Hertz計(jì)算方法中單元格的個(gè)數(shù)不能超過27個(gè)，否則會(huì)影響計(jì)算結(jié)果。

圖5 網(wǎng)格劃分與彈性趨近量之間的關(guān)系

主曲率半徑與單元格中心最大接觸應(yīng)力之間的關(guān)系如圖6所示。由圖6可知，主曲率半徑改變時(shí)，2種計(jì)算方法中的單元格中心處的最大接觸應(yīng)力也隨著改變，隨著主曲率半徑逐漸增大而呈明顯上升趨勢，當(dāng)主曲率半徑為70～80 mm之間時(shí)，最大接觸應(yīng)力有下降趨勢，在80 mm之后又逐漸增大。

圖6 主曲率半徑與單元格中心最大接觸應(yīng)力之間的關(guān)系

5 結(jié)論

彈性趨近量和接觸應(yīng)力是后期求解軸承剛度的基礎(chǔ)，在2種計(jì)算方法中，對結(jié)果有影響的主要參數(shù)包括主曲率半徑和網(wǎng)格劃分個(gè)數(shù)，隨著網(wǎng)格劃分個(gè)數(shù)逐漸增多，對應(yīng)的單元格中心處的最大接觸應(yīng)力逐漸增大，雖然有個(gè)別突然減小，但整體上都呈上升趨勢；對應(yīng)的彈性趨近量也隨著單元格個(gè)數(shù)增多而逐漸增大。當(dāng)主曲率半徑逐漸增大時(shí)，兩者的彈性趨近量都逐漸減小，其中非Hertz計(jì)算方法減小比較明顯，而Hertz計(jì)算方法中減小幅度較小，2種曲線雖然看上去差別較大，但數(shù)值相差較小。2種計(jì)算方法雖然存在一定的誤差，但誤差較小，可以滿足實(shí)際要求。主曲率半徑對結(jié)果的影響在后期軸承的選擇上有一定的借鑒意義，根據(jù)需要對網(wǎng)格劃分個(gè)數(shù)進(jìn)行合理的劃分和求解。