基于ABAQUS的列車軸箱軸承動力學分析

郝燁江，李強，鄭靜

(北京交通大學 機械與電子控制工程學院，北京 100044)

符號說明

B——軸承寬度

C0——軸承額定靜載荷

d——軸承內徑

D——軸承外徑

Dpw——滾子組節圓直徑

Fa——軸向載荷

Fap——軸承的軸向承載能力

Fc——滾子離心力

FL——轉臂效應載荷

Fr——徑向載荷

Ft——牽引制動載荷

g——重力加速度

Gr——徑向游隙

Hz——軸承承載區寬度

I——節點內力

k1，k2——載荷系數

K——載荷-位移系數

Lt——滾子有效承載長度

Lz——轉臂長度

ma——軸重

ms——簧下質量

M——節點質量矩陣

n——軸承轉速

ni——內圈轉速

nc——保持架轉速

P——系統外力

Qi——第i個滾子的徑向載荷

ü——節點加速度矩陣

Z——滾子數

α——接觸角

γ——計算參數

Δt——時間增量

δi——第i個滾子的徑向位移

δr——沿徑向力方向產生的總位移

φi——第i個滾子的方位角

列車軸箱軸承是確保列車運行品質和安全性的關鍵技術部件，通常使用雙列圓錐滾子軸承或雙列圓柱滾子軸承。由于結構及使用工況的特點，軸承在運行周期內承受的載荷十分復雜，受多方面因素的影響。軸承失效形式多種多樣，但大多數集中在表面接觸疲勞損傷和次表面缺陷疲勞擴展方面。對于軸承，除保證良好工作條件并定期檢修外，針對性地研究其在不同工況下的動力學特性對于改進結構、提高壽命意義重大。文中以BC2-0375型列車軸箱軸承為對象建立分析模型，按照實際載荷工況加載邊界條件，利用有限元方法對其動力學特性進行求解分析[1-3]。

1 軸箱軸承徑向載荷

軸箱軸承承載狀況極為復雜，受到輪軸、轉臂等多方面因素的影響(圖1)，很難給出完全符合實際條件的載荷計算公式，通常采用的徑向力簡化計算式為

Fr=(mag-msg)/2+Ft+FL。

(1)

圖1 軸箱軸承受力示意圖

簧下質量根據列車型號基本固定不變；軸重取列車空載狀態下的質量；牽引制動載荷通過牽引功率和制動加速度求得；轉臂效應載荷由轉臂受到軸向載荷導致的軸承偏轉產生(圖2)，即

FL=FaLz/Hz。

(2)

圖2 轉臂效應載荷計算示意圖

參考軸承幾何尺寸并根據實際載荷狀況和潤滑條件的影響對參數加以修正，得到徑向力的最終計算式為

Fr={-Fap+104k1C0/[n(d+D)]}/k2。

(3)

BC2-0375型列車軸箱軸承為圓柱滾子軸承，其簡化模型的幾何參數為：d=130 mm,D=240 mm,B=80 mm; 按照n=900 r/min，C0=1 560 kN，代入相應的載荷系數，由(3)式得到穩定運行情況下等效的徑向載荷為25 kN，由此得到計算模型加載的載荷邊界條件[4-6]。

2 工作狀態下軸承內部載荷分布

在徑向力與轉速恒定的情況下，軸承會形成穩定的承載區間，如圖3所示。通常轉速不是很高的時候，滾子的離心力、摩擦力和力矩對載荷分布不會產生特別明顯的影響。數值計算的基本思路是建立滾子與滾道間的載荷-位移關系，通過迭代求解非線性方程組獲得位移分布，從而最終獲得軸承內部的載荷分布。

圖3 軸承內部承載區示意圖

假定軸箱軸承徑向游隙為零，根據幾何尺寸可得計算參數γ為

γ=Dcosα/Dpw=0.144 385，

(4)

保持架和內圈轉速間的運動關系為

nc=0.5ni(1-γ)=0.427 81ni，

(5)

繞軸承軸線以nm轉速運動的鋼制滾子的離心力為

(6)

在低速條件下，離心力對于滾子與滾道間的接觸載荷貢獻很小；但在高速時離心力的效應顯著。由于所計算軸承的轉速集中在1 200 r/min以內，產生的離心力效應基本可以忽略。

在軸承的徑向載荷作用方向，根據力的平衡條件可得

(7)

根據載荷-位移系數換算得到

(8)

軸承工作狀態下內、外滾道與滾子間的位移關系為

(9)

聯立 (8)～(9)式，建立由18個方程組成的非線性方程組，通過Newton-Rampson法對其進行迭代求解，得到滾子與外圈滾道的位移分布，再利用載荷位移系數進行修正就可得到軸承的內部載荷分布。將上述求解過程利用MATLAB程序語言進行描述，分別對420和900 r/min轉速下的軸承內部載荷分布進行數值計算，得到的結果如圖4所示。

圖4 軸承內部載荷分布

由圖4可知，最大載荷出現在方位角為0的位置，載荷隨著方位角的增大而減小，直至承載區之外滾子不受載荷作用，承載區域約為168°。中低轉速區間，轉速的提升對于軸承外圈所受載荷的大小及分布影響并不顯著，承載區域沒有發生大的變化。在忽略滾子滾動摩擦效應的情形下，位于承載區各位置滾子的接觸應力可以按Hertz接觸理論計算得出。方位角為0處的內圈的接觸應力為625.545 MPa，接觸區域半寬為0.126 8 mm；外圈接觸應力為542.174 2 MPa，接觸區域半寬為0.142 9 mm，這些結果是進行有限元離散求解的參考基礎[7-9]。

3 基于ABAQUS的顯式動力學計算

ABAQUS-explicit是一種高效的求解復雜非線性問題的數值計算工具。接觸問題是一個典型的復雜非線性問題，本構關系即使是處于線性區，接觸區域的局部效應也會使得接觸應力通常到達塑性區，尤其對于軸承來說，接觸應力達到1 380 MPa以上的情況并不少見；再加之其所處的邊界條件的非線性，無法事先知道應力的分布，接觸應力與接觸區域互相影響，呈現出高度的非線性特質。對于靜態基本幾何外形的非共形接觸，已經由Hertz理論得到了較為精確的解答。

軸承的動態接觸問題較之經典接觸對之間的接觸要復雜許多，軸承零件之間的接觸關系時刻發生變化，且受到內部載荷分布、潤滑效應以及材料等多方面因素的影響。使用隱式求解方法進行解答可以平穩地建立接觸關系，減小沖擊效應產生的應力波的影響，但對于接觸對建立條件的要求復雜很多，而顯式計算中的通用接觸模式(general contact)可以很方便地建立這種接觸關系，且不存在收斂問題，適合于求解內部結構關系多變且存在復雜相互接觸作用的瞬時動態響應問題。但是ABAQUS-explicit對于求解軸承接觸問題的缺點在于初始力邊界條件施加時會產生應力波干擾，這種沖擊的效應需要加大載荷步長度來進行削弱，計算量也隨之上升[10]。

3.1 ABAQUS-explicit求解器算法描述

ABAQUS-explicit應用中心差分算法對運動方程進行顯式的時間積分，由一個增量步的動力學條件計算下一個增量步的動力學條件。在增量步開始時，求解動力學平衡方程

Mü=P-I，

(10)

在當前增量步開始時，計算加速度

ü|(t)=M-1·(P-I)|(t)。

(11)

顯式算法并不需要同時求解聯立方程組，任何節點的加速度完全取決于節點質量和作用在節點上的合力。對加速度在時間上進行積分采用的是中心差分法，在計算速度變化時假定加速度為常數,將速度的變化值加上前一個增量步(increment)中間時刻的速度來確定當前增量步中間時刻的速度，即

(12)

再利用速度對時間的積分與增量步初始位移求和,即可得到增量步末的位移

(13)

得到節點位移后，即可通過應變速率和本構關系獲得節點的應力矩陣[11]。

3.2 有限元模型的建立

軸承運行過程中，滾子與內、外滾道及保持架之間存在著復雜的運動和接觸關系，且由于保持架和滾子之間的碰撞以及內、外滾道之間的阻尼效應，往往內圈和保持架轉動并不同步，這也就影響到了滾子承受循環載荷的頻率。為了盡量真實模擬軸承零件的動力學關系，將滾子，內、外滾道以及保持架均按照彈性體進行建模。

取單列軸承為對象，根據幾何參數利用ABAQUS集成的modeling模塊完成幾何模型的建立，忽略倒角、注油孔以及一些邊緣細節。徑向游隙為0，保持架與滾子間的間隙為0.5 mm。忽略微小區域塑性應變，材料選定在線彈性范圍內。單元類型選擇減縮積分單元C3D8R，掃掠形成網格。需要注意的是由于接觸區域附近的離散化質量直接決定了計算精度，為了得到準確的接觸應力，根據大量接觸問題的有限元計算結果，取接觸區域的單元尺寸為接觸區域半寬的一半，如果是橢圓形接觸區域，則需小于短半軸長度。離散后的有限元模型如圖5a所示。

圖5 軸承有限元離散結構

滾動體接觸區域附近網格的質量也會對應力分布產生很大的影響，故將滾子分為2個環形區域，外環形區域加密網格，保證均勻度;內環形區域疏化,減小計算規模，如圖5b所示，單元數目為138 939，節點數目為176 517。

內、外滾道以及滾子彈性模量為207 GPa，密度為7.8×10-6kg/mm3，泊松比為0.3。保持架彈性模量為120 GPa，密度為8.5×10-6kg/mm3，泊松比為0.34。

利用顯式計算的通用接觸建立軸承接觸模型，在運算過程中，ABAQUS將所有的零件生成一個實體模型，通過設置自接觸來判斷接觸區域，從而實現復雜接觸模型的建立。法向接觸屬性采用硬接觸，接觸剛度和罰函數系數將根據實際接觸區域的變形自動調整。切向接觸屬性中的摩擦模型采用庫侖摩擦，內、外圈滾道與滾子之間靜摩擦因數給定為0.1,動摩擦因數給定為0.05，衰減系數取0.01；滾子與保持架之間摩擦因數在中低轉速區間對結果影響不大，簡化后給定為0.002。

邊界條件：內圈節點耦合在圓心處,在初始分析步中僅放開沿y軸(垂向)的徑向自由度；保持架耦合在圓心處，只釋放繞軸線的旋轉自由度；外圈始終施加全約束；分析全程設定5個載荷步，step1至step5：step1(0.01 s)中對滾子施加重力加速度，使接觸可以穩定的建立；step2(0.02 s)中在內圈耦合點施加垂向載荷25 000 N；step3(0.1 s)放開內圈耦合點的繞軸線的旋轉自由度，施加200 r/min轉速;step4(0.1 s)在step3的基礎上繼續增加轉速達到420 r/min;step5(0.1 s)繼續將轉速增加到900 r/min。

所有的載荷施加形式均選定光滑曲線(smooth step)加載，盡量減小系統應力波的震蕩。

3.3 計算結果分析

顯式計算的速率與節點質量的開方值成反比，網格越精細，計算所需的時間增量步越多，計算時間越長，復雜的接觸模型計算消耗的時間較長，可以選擇相應的質量放大系數，在保持精度的條件下實現軸承接觸的動力學仿真。

(1)不同轉速等級下軸承內部的載荷分布情況與Newton-Rampson算法得到的載荷分布結果基本一致，且在加速過程中承載區也較為穩定，在中低速時滾子與滾道接觸產生的等效應力場水平并沒發生很大的變化，基本維持在1 000 MPa以內，且最大值出現在滾子邊緣處，這是由邊緣壓膜效應產生的應力集中現象造成的，如圖6所示。各個載荷步最終的應力分布如圖7所示。當轉速提升到900 r/min,平均應力水平增大，某些瞬時時刻的最大應力出現在承載區的滾子與保持架之間(由碰撞所產生)。這說明軸承持續加速過渡階段保持架起到了重要的作用，其強度決定了加速運轉的可靠性與安全性。

圖6 滾子接觸應力分布

輸出外圈內表面與對稱中面交線中點即6 027號節點(圖8)在不同轉速下的等效應力曲線，如圖9所示。

由圖9可知，隨著轉速升高，單位時間內的應力峰值數目增多，受載頻率變大，但峰值并未顯著增大，說明中低轉速產生的離心力作用并不顯著。

圖8 第6 027號節點位置

圖9 第6 027號節點在不同轉速下的等效應力曲線

step2中徑向力加載穩定后，外圈與方位角為0的滾子所構成的接觸對中，外圈接觸區域中心平均接觸應力為480.309 MPa，滾動體接觸區域中心平均接觸應力為645 MPa；內圈與方位角為0的滾子構成的接觸對中，內圈接觸區域中心平均接觸應力為630.702 MPa，滾動體接觸區域中心平均接觸應力為723.5 MPa(圖10)。理論上同一個接觸對，接觸區域的應力對于兩接觸體應相同，產生的差值由離散化計算導致，網格的密度與均勻程度決定了接觸應力計算的準確度，內、外圈接觸區域網格的大小為2 mm，約為滾子的10倍，導致滾子接觸應力偏大，滾道接觸應力偏小，將對應節點的接觸應力取平均值作為接觸應力數值計算的結果。則滾子與內、外圈間的接觸應力分別為677.101，562.65 MPa，與理論計算值相差分別為7%，3%。

(2)軸承即使是理想幾何形狀，也會由于各個滾子角位置的不斷變化而產生內、外圈的周期性相對運動，其振動頻率等于保持架運動速度與滾子數之積，即滾子通過外圈滾道的頻率，類似于彈簧的振動。軸承的振動頻率往往影響著轉向架系統尤其是軸向輪對及牽引系統的振動頻率，稱之為可變彈性柔度。采集step1至step5的內圈耦合點的垂向位移，如圖11所示。

圖10 軸承接觸應力分布

圖11 內圈耦合點的垂向位移時間歷程

由圖11可知，在5個載荷步的時間歷程中，step1重力加載時耦合點并未出現顯著的垂向位移；step2徑向力加載時出現較大的垂向位移，末端有一段時間的應力波激擾衰減過程；step3轉速施加后，位移呈現周期性變化，轉速為200 r/min時，軸承的可變彈性柔度頻率為120.48 Hz，與理論計算得到的頻率(145.758 Hz)基本吻合。且隨著轉速的提高，內圈耦合點位移變化周期逐漸減小，但是趨勢并不明顯，原因在于使用了平滑曲線加載，在載荷步末端才達到目標轉速，且計算時間較短，并沒有顯示出較為明顯的頻率變化。隨著轉速提高，由于離心力作用加大以及內圈的變形導致的垂向位移波動振幅變大，趨勢越發清晰[12-13]。這種高頻的振動需要在信號采集結果處理時濾掉，否則會影響分析精度。

4 結束語

(1)根據列車軸箱軸承結構特點和受力分析求解了單列軸承承受的徑向載荷。

(2)采用Newton-Rampson法，利用MATLAB求解了轉動狀態下軸承內部載荷的分布情況，并根據Hertz接觸理論求解了最大載荷作用接觸區域的接觸應力，將其作為有限元計算的參考。

(3)建立了有限元離散模型，加載計算得到的載荷和相應的轉速工況，研究軸承的動力學特性，獲得了不同轉速條件下軸承內部應力場與位移場，得到的等效應力分布與實際工況和理論解相符合，實現了軸箱軸承的顯式動力學模擬；并進行了相關動力學特性分析，以此為基礎，可以對軸箱軸承各種運行失效狀態進行模擬仿真，找出各種不同失效形式下振動和應力分布的特點，為故障診斷提供可靠的參考依據，并可借此提出相應的軸承在線檢測方法。