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

段磊，陳曉陽，張濤

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

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

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

1 球-板模型有限元分析

1.1 建模和網(wǎng)格劃分

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

圖1 球-板接觸模型

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

圖2 球-板接觸模型網(wǎng)格劃分

1.2 約束與載荷

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

1.3 接觸和邊界條件

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

1.4 結(jié)果分析

接觸區(qū)域球的單個網(wǎng)格尺寸為a/3.9時，球-板模型接觸區(qū)域應(yīng)力云圖如圖3、圖4所示。

圖3 球應(yīng)力云圖

圖4 板應(yīng)力云圖

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

通過Patran后處理提取球和板在接觸面上接觸區(qū)域的接觸點編號及應(yīng)力如圖5、圖6所示。圖中，括號內(nèi)為節(jié)點的法向接觸應(yīng)力，括號上面為對應(yīng)的節(jié)點編號。

圖5 球在接觸區(qū)域的網(wǎng)格

圖6 板在接觸區(qū)域的網(wǎng)格

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

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

圖7 球法向接觸應(yīng)力隨接觸邊界變化圖

圖8 板法向接觸應(yīng)力隨接觸邊界變化圖

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

圖9 球位移云圖

圖10 板位移云圖

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

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

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

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

2 實例計算

2.1 建模

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

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

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

2.2 設(shè)置邊界條件

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

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

圖12 軸承有限元模型

2.3 結(jié)果分析

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

圖13 內(nèi)圈上接觸應(yīng)力云圖

圖14 內(nèi)圈上接觸位移云圖

圖15 球上接觸應(yīng)力云圖

圖16 球上接觸位移云圖

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

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

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

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

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

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

3 結(jié)束語

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

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

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

3)網(wǎng)格劃分對接觸區(qū)域、接觸應(yīng)力影響很大，為保證最大接觸應(yīng)力準(zhǔn)確，對于軸承劃分網(wǎng)格，應(yīng)以受力最大球的接觸區(qū)域為基準(zhǔn)。

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