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

李猛帥，鄭甲紅

（陜西科技大學，陜西 西安 710021）

1 確定相關參數

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

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

采用文獻[1]中的式1.21—式1.22b得出：

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

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

這樣就確定了Z的表達式，在MATLAB軟件中輸入時，需要將其以向量的形式輸入，而且它的個數也要和2種數值求解方法中的網格數一一對應，因為這樣才能滿足矩陣之間的基本運算，從而進一步求解基本平衡方程。

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

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

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

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

圖1 劃分矩形單元

則單元j上的接觸應力在任意一個單元i的中心上產生的位移為：

式（6）中：Dij為柔度系數。

柔度系數Dij通過MATLAB軟件進行計算，而且它的下標在MATLAB中代表著任意一個矩陣，也就是單元矩陣劃分多少個，那么對應的柔度系數矩陣在MATLAB中就有多少個。那么線接觸的基本方程就可以表示為：

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

需要注意的是，在計算時要滿足pj≥0，因為當pj＜0時，代表此單元格已經超出實際接觸區域，沒有受到外載荷的作用，在下一次程序迭代時就不再進行計算。

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

3 Hertz接觸問題的求解方法

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

它的求解思想也是將所有可能的接觸區域進行單元格劃分，不同的是在單元j內假定接觸應力沿滾子素線方向（y軸）為均勻分布，沿橫向（x軸）按Hertz分布[1]，條形單元劃分如圖2所示。

圖2 條形單元劃分

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

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

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

Hertz線接觸的基本方程為：

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

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

4 實例計算

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

根據編制好的2種計算程序，將所有初始條件輸入到MATLAB軟件中，通過對2種計算方法的仔細分析，把對結果有直接影響的參數一一列舉，并整合出曲線圖。網格劃分個數與單元格中心最大接觸應力之間的關系如圖3所示。

圖3 網格劃分與最大接觸應力之間的關系

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

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

圖4 主曲率半徑與彈性趨近量之間的關系

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

圖5 網格劃分與彈性趨近量之間的關系

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

圖6 主曲率半徑與單元格中心最大接觸應力之間的關系

5 結論

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