基于Nastran的球軸承接觸問題有限元分析

段磊，陳曉陽，張濤

(上海大學 機電工程與自動化學院，上海 200072)

滾動軸承是機械傳動中重要的零部件[1]，軸承接觸應力影響著軸承的接觸疲勞和磨損，決定著軸承的可靠性壽命。

滾動軸承內、外圈和滾動體之間的接觸面為橢圓。目前Hertz理論法相對成熟[2]，但要求作用力與接觸面垂直是其局限性，不能計算考慮摩擦問題的軸承接觸問題。不少學者利用有限元軟件對軸承進行接觸分析[3-6]，當2個物體接觸時，其最大接觸應力、接觸變形不同，且最大接觸應力也不在接觸點上，不符合Hertz接觸理論。針對這一問題，文中首先利用Patran & Nastran軟件,建立球-板有限元接觸模型,探討有限元網格的劃分方法以及最大接觸應力和接觸區域形狀不等的原因，然后以某型深溝球軸承為例，說明此有限元網格劃分方法的可行性。

1 球-板模型有限元分析

1.1 建模和網格劃分

對62306-2RS軸承進行建模分析，球半徑為5 mm，板長寬高均為0.5 mm，作用力為200 N，材料均為GCr15，彈性模量為207 GPa，泊松比為0.3。由于球與板接觸具有對稱性，故只對1/4球-板模型進行分析，其模型如圖1所示。

圖1 球-板接觸模型

為兼顧計算精度和速度，采用六面體劃分網格，非接觸區域網格劃分較粗，接觸區域網格劃分較細。采用Patran軟件對球體劃分六面體網格的要求較嚴格[7]，如果直接在一個面上劃分面網格，然后掃掠為體網格，則掃掠軸線位置為五面體網格；用自動生成網格功能劃分則有四、六面體存在。因此，對1/4球體不能直接劃為完全六面體網格。首先將球分割為2部分，中間部分用meshing中sweep>loft劃分，另一部分可先在面上劃分六邊形網格，再旋轉得到六面體網格。板為規則形狀，可直接劃分六面體網格。此外，采用布種子節點的方法來控制接觸區域網格尺寸，如圖2所示。通過調整網格數量，獲得30組網格疏密不同的接觸有限元模型，將求解結果與Hertz理論結果對比發現，當球的單個網格尺寸等于a/3.9時(a為橢圓接觸區域的長半軸)，接觸形狀、最大接觸應力和彈性趨近量最接近Hertz理論解。

圖2 球-板接觸模型網格劃分

1.2 約束與載荷

將垂直于z軸的圓平面設為剛性面，然后在該平面上沿z軸施加50 N的載荷(1/4模型)。分別約束模型在x，y方向的位移以及板底面的位移。

1.3 接觸和邊界條件

球和板均設置為變形體；接觸問題為邊界非線性問題，未考慮摩擦因素，不設置摩擦類型和摩擦參數；接觸區域應變很小，非線性幾何影響因素設置為小位移小應變類型；設置輸出節點和單元的力、位移、應力、法向接觸力、法向接觸應力。

1.4 結果分析

接觸區域球的單個網格尺寸為a/3.9時，球-板模型接觸區域應力云圖如圖3、圖4所示。

圖3 球應力云圖

圖4 板應力云圖

由圖可知，球的最大接觸應力為 2 650 MPa，不在初始接觸點(圖2中黑點處)。板的最大接觸應力為2 840 MPa，在初始接觸點上。這與Hertz理論分析結果不同。其原因如下：有限元法是先把物體分為有限個單元體，計算出總體剛度矩陣K和所有節點的等效節點力矩陣F，然后根據F=Kδ得出所有節點的位移矩陣δ，最后利用應力矩陣σ=Sδε(S為常數矩陣，δε表示某個單元的節點位移矩陣，指數ε與單元有關)算出單元的應力矩陣σ，但計算結果要求得出單元節點的應力。根據有限元近似解性質[8]，應力和應變近似解一定是在精確解上下浮動，單元的應力在高斯積分點上，然后利用單元形函數將高斯節點的應力插值到單元節點上。這樣相鄰單元僅在相鄰的節點與邊界上變形連續，而非相鄰的邊界上并不一定相等，導致兩相鄰單元在相鄰節點上計算的應力不連續，一般采用繞節點平均或應力磨平等方法進行處理。

通過Patran后處理提取球和板在接觸面上接觸區域的接觸點編號及應力如圖5、圖6所示。圖中，括號內為節點的法向接觸應力，括號上面為對應的節點編號。

圖5 球在接觸區域的網格

圖6 板在接觸區域的網格

由圖可知，球和板的最大接觸應力分別為2 645 MPa和2 843 MPa，分別在節點233和6 607上。通過有限元理論可知，節點233和6 607上的接觸應力均是由對應網格上高斯節點的應力等效得到。單元上高斯點應力插值和繞節點平均或應力磨平等方法導致其最大接觸應力不在初始接觸點232上。因此造成球和板上最大接觸應力數值不等和位置(初始接觸點)不同。

球和板法向接觸應力隨邊界變化如圖7、圖8所示。由圖可知，球和板接觸區域長半軸a分別為0.193 2 mm和0.211 4 mm。通過Patran后處理提取球和板接觸面上接觸應力可得，球和板所有節點的接觸應力之和均為49.9 N，約等于施加的外力50 N，可以說明在接觸面上球和板的等效節點力F基本相同，又由于球和板在接觸區域的網格不同，則其剛度矩陣K不同，由F=Kδ可知，接觸區域半徑會有微小差別，因此造成球和板接觸區域長半軸不等。

圖7 球法向接觸應力隨接觸邊界變化圖

圖8 板法向接觸應力隨接觸邊界變化圖

接觸區域球的單個網格尺寸為a/3.9時，球-板模型接觸區域位移云圖如圖9、圖10所示。

圖9 球位移云圖

圖10 板位移云圖

由圖可知，板的絕對最大位移為2.93 μm，在初始接觸點6 007上，球的絕對最大位移為6.54 μm，在球心點上。則板的中心(無窮遠處)相對于初始接觸點的位移為2.93 μm，球的中心相對于初始接觸點的變形量為3.61 μm。由Hertz接觸理論可知，2個接觸物體的彈性趨近量為初始接觸點處各自所對應的圓弧中心在接觸過程中所移動的距離之和，球-板接觸模型的彈性趨近量為6.54 μm。

由于單元網格大小和形函數不同，最大接觸應力值不同。球-板模型最大接觸應力及接觸半徑取球和板的平均值，其與Hertz理論結果對比見表1。

表1 球-板有限元解與解析解

由表1可知，有限元解與 Hertz 理論解誤差在工程允許的范圍內，故利用有限元軟件來解決接觸問題是可行的。

2 實例計算

2.1 建模

以某型深溝球軸承模型為例，其參數為：外徑D=9 mm，內徑d=4 mm，寬度B=2.5 mm，球數n=8，球徑Dw=1.3 mm，內外圈和球的材料為GCr15，彈性模量為207 GPa，泊松比為0.3，額定靜載荷189 N，仿真中徑向載荷取40 N。

由于軸承具有對稱性，且僅有3個球受力，為減少網格數量，對其3/8模型進行有限元分析。利用上述方法進行網格劃分得出其模型如圖11所示。

圖11 3/8軸承有限元模型

2.2 設置邊界條件

為模擬軸對內圈的作用力，將內圈內表面設置為剛性，然后在軸承中心點上沿x軸(圖11)方向施加40 N載荷。

外圈外表面所有節點約束x方向位移；內外圈端面所有節點約束z方向的位移；在柱坐標系下約束球與內外圈接觸點連線上所有節點的周向和軸向自由度；約束球在xOy平面上所有節點，z方向的位移，如圖 12所示。

圖12 軸承有限元模型

2.3 結果分析

通過Nastran軟件仿真計算后得到球與內圈接觸的應力和位移云圖如圖13～圖16所示。

圖13 內圈上接觸應力云圖

圖14 內圈上接觸位移云圖

圖15 球上接觸應力云圖

圖16 球上接觸位移云圖

由圖可知，球與內圈的接觸區域為橢圓，接觸應力和變形由接觸中心向外逐漸減小，符合實際情況。球與外圈接觸情況與此相同。

軸承中受力最大球的最大接觸應力及變形有限元解與解析解的對比見表 2。有限元解與解析解誤差在工程允許范圍內，故利用有限元軟件來解決軸承接觸問題能滿足工程實際應用，而且可以圖形的形式直觀地顯示接觸應力、接觸尺寸、彈性趨近量等。

表2 受力最大球有限元解與解析解比較

軸承中受力最大球相鄰2個球的最大接觸應力、變形有限元解與軸承解析解見表 3。

表3 受力最大球相鄰球有限元解與解析解

通過對軸承中與受力最大球相鄰的2個球的有限元解和解析解對比可知，其誤差較大。這是由于接觸區域網格是以最大球與內圈接觸為基準劃分的，相對相鄰球接觸區域網格劃分更粗，誤差相對較大。

3 結束語

通過有限元法與軸承解析法對比可知，利用有限元軟件來解決軸承接觸問題是可行的。

1)采用有限元軟件對多組不同網格尺寸的球-板接觸模型分析可知，當球的單個網格尺寸等于a/3.9時，得到最佳結果。

2)在有限元軟件分析中，2個接觸物體同一接觸點最大接觸應力不同和最大接觸應力不在初始接觸點是由于網格差異和軟件插分算法引起的。

3)網格劃分對接觸區域、接觸應力影響很大，為保證最大接觸應力準確，對于軸承劃分網格，應以受力最大球的接觸區域為基準。

由于球-板模型和軸承模型彈性趨近量誤差相對較大，需進一步分析原因。