郝燁江，李強，鄭靜

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

符號說明

B——軸承寬度

C0——軸承額定靜載荷

d——軸承內徑

D——軸承外徑

Dpw——滾子組節(jié)圓直徑

Fa——軸向載荷

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

Fc——滾子離心力

FL——轉臂效應載荷

Fr——徑向載荷

Ft——牽引制動載荷

g——重力加速度

Gr——徑向游隙

Hz——軸承承載區(qū)寬度

I——節(jié)點內力

k1，k2——載荷系數(shù)

K——載荷-位移系數(shù)

Lt——滾子有效承載長度

Lz——轉臂長度

ma——軸重

ms——簧下質量

M——節(jié)點質量矩陣

n——軸承轉速

ni——內圈轉速

nc——保持架轉速

P——系統(tǒng)外力

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

ü——節(jié)點加速度矩陣

Z——滾子數(shù)

α——接觸角

γ——計算參數(shù)

Δt——時間增量

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

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

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

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

1 軸箱軸承徑向載荷

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

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

(1)

圖1 軸箱軸承受力示意圖

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

FL=FaLz/Hz。

(2)

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

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

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

(3)

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

2 工作狀態(tài)下軸承內部載荷分布

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

圖3 軸承內部承載區(qū)示意圖

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

γ=Dcosα/Dpw=0.144 385，

(4)

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

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

(5)

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

(6)

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

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

(7)

根據(jù)載荷-位移系數(shù)換算得到

(8)

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

(9)

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

圖4 軸承內部載荷分布

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

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

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

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

3.1 ABAQUS-explicit求解器算法描述

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

Mü=P-I，

(10)

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

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

(11)

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

(12)

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

(13)

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

3.2 有限元模型的建立

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

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

圖5 軸承有限元離散結構

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

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

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

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

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

3.3 計算結果分析

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

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

圖6 滾子接觸應力分布

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

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

圖8 第6 027號節(jié)點位置

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

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

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

圖10 軸承接觸應力分布

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

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

4 結束語

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

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

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