局部剝落故障對(duì)滾動(dòng)軸承接觸與振動(dòng)特性的影響

劉靜，唐昌柯，師志峰，王林峰，徐子旦，邵毅敏

(1.重慶大學(xué) 機(jī)械傳動(dòng)國家重點(diǎn)實(shí)驗(yàn)室，重慶，400044；2.重慶大學(xué)機(jī)械工程學(xué)院，重慶，400044)

滾動(dòng)軸承是重要基礎(chǔ)部件之一，廣泛應(yīng)用于工業(yè)、能源和航天等行業(yè)的關(guān)鍵機(jī)械裝備。滾動(dòng)軸承的運(yùn)行狀態(tài)和壽命直接影響機(jī)械裝備的運(yùn)轉(zhuǎn)狀態(tài)，30%的旋轉(zhuǎn)機(jī)械故障是由軸承故障引起的[1-3]。滾動(dòng)軸承的故障形式包括銹蝕、斷裂、磨損、疲勞剝落和膠合等，許多研究者對(duì)滾動(dòng)軸承滾道表面的局部剝落故障進(jìn)行研究。當(dāng)滾動(dòng)軸承出現(xiàn)局部剝落故障時(shí)，會(huì)導(dǎo)致軸承剛度發(fā)生變化，使得軸承的振動(dòng)水平增大[4-5]，甚至危害整套機(jī)械裝備。因此，需要深入研究局部剝落故障誘發(fā)的滾動(dòng)軸承振動(dòng)特征。對(duì)滾子與局部剝落故障周圍區(qū)域接觸特性的深入認(rèn)識(shí)將有助于理解軸承局部剝落故障的振動(dòng)特征，為軸承的故障診斷和振動(dòng)特征識(shí)別提供有益的參考。針對(duì)滾動(dòng)軸承部件局部剝落故障與滾子的接觸特性問題，許多學(xué)者進(jìn)行了大量研究工作。BRANCH等[6]利用試驗(yàn)和有限元的方法，研究了球軸承滾動(dòng)體沖擊剝落故障邊緣時(shí)引起的臨界應(yīng)力和塑性應(yīng)變。PATIL 等[7]建立了帶有局部剝落故障的球軸承動(dòng)力學(xué)模型，研究了局部剝落故障尺寸對(duì)軸承振動(dòng)特征的影響規(guī)律。牛藺楷等[8]以剛體動(dòng)力學(xué)碰撞理論為基礎(chǔ)，在分析滾動(dòng)軸承運(yùn)動(dòng)學(xué)特性的基礎(chǔ)上，對(duì)滾動(dòng)體通過局部剝落故障時(shí)產(chǎn)生的沖擊力進(jìn)行了分析和計(jì)算。樊莉等[9]運(yùn)用ANSYS/LS-DYNA軟件，建立了含局部剝落故障的滾動(dòng)軸承有限元模型，研究了局部剝落故障對(duì)軸承滾動(dòng)體與滾道之間接觸應(yīng)力的影響規(guī)律。張福星等[10]利用有限元方法，分析了深溝球軸承滾動(dòng)體和內(nèi)、外圈不同部位接觸應(yīng)力和變形的分布。張雪等[11]通過計(jì)算軸承接觸表面及次表面應(yīng)力等證明了有限元分析的正確性。伍濟(jì)鋼等[12]將軸承的表面損傷故障考慮為沖擊函數(shù)分別加載到軸承的內(nèi)外圈及滾動(dòng)體上，計(jì)算了軸承振動(dòng)和頻譜分析。剡昌鋒等[13]計(jì)算了考慮彈流潤滑及滑動(dòng)作用下滾動(dòng)軸承系統(tǒng)局部缺陷位移激勵(lì)的動(dòng)力學(xué)響應(yīng)。盡管上述研究采用集中質(zhì)量法和有限元分析方法討論了局部剝落故障過渡區(qū)對(duì)滾子與滾道之間接觸應(yīng)力及其誘發(fā)的振動(dòng)特征，但是缺乏滾子通過軸承局部缺陷故障時(shí)滾子與滾道之間接觸特性的相關(guān)研究。針對(duì)上述問題，本文作者采用有限元分析方法，建立圓柱滾子軸承滾子與滾道接觸二維有限元等效模型，分析單個(gè)滾子進(jìn)入和退出滾道局部剝落故障區(qū)域的過程中滾子與滾道之間的接觸特性。研究局部剝落故障寬度對(duì)滾子與滾道之間的接觸區(qū)寬度、接觸變形、接觸應(yīng)力分布狀態(tài)與接觸剛度的影響規(guī)律，討論局部剝落故障過渡區(qū)對(duì)軸承振動(dòng)特性的影響規(guī)律，為深入認(rèn)識(shí)局部剝落故障誘發(fā)的滾動(dòng)軸承的振動(dòng)特征提供參考。

1 模型描述

1.1 圓柱滾子與滾道等效接觸模型

圖1所示為局部剝落故障出現(xiàn)時(shí)圓柱滾子通過局部剝落故障區(qū)域的過程。圖1中：Q為作用于滾子上的徑向載荷；v為滾子運(yùn)動(dòng)方向；O1和O2分別為滾子在通過局部剝落故障過程中的不同位置；δ1和δ2分別為不同位置的滾子與滾道之間的接觸變形量；w1和w2分別為不同接觸位置處滾子與滾道之間的接觸寬度。本文假設(shè)滾道表面的局部剝落故障采用矩形輪廓進(jìn)行表征[4]。滾子進(jìn)入局部剝落故障區(qū)間時(shí)，滾子與故障邊緣位置發(fā)生接觸。由于故障邊緣的曲率半徑不同于正常滾道，故滾子與故障滾道的接觸變形和滾子與正常滾道的接觸變形量存在差異，滾子與故障滾道的接觸寬度和滾子與正常滾道的接觸寬度也存在差異。為了準(zhǔn)確地描述滾子通過局部剝落故障過程中滾子與滾道的接觸狀態(tài)，忽略滾子在運(yùn)動(dòng)方向的速度沖擊對(duì)接觸的影響，建立了圓柱滾子軸承滾子與滾道接觸的等效模型，如圖2所示。

圖1 圓柱滾子與滾道之間的位置關(guān)系Fig.1 Geometric relationship between roller and defect

圖2 圓柱滾子軸承局部剝落故障等效模型Fig.2 Equivalent model of cylindrical roller bearing with localized spalling defect

圓柱滾子軸承滾子與滾道接觸的等效模型[14]，如圖2所示。從圖2可知：滾道表面分為局部剝落故障區(qū)、過渡區(qū)與無故障區(qū)，矩形缺陷表示軸承局部剝落故障，XP為等效模型中軸承滾子與滾道接觸的中心位置與軸承局部剝落故障中心位置的距離；Pmax為等效滾子中心位置的最大赫茲接觸應(yīng)力。等效模型的長度為10b，寬度為7b，其中，b為用赫茲接觸模型計(jì)算得到的軸承接觸半寬。接觸半寬b表示為[15]

式中：F為作用在軸承上的徑向載荷；R1為滾子半徑；R2為內(nèi)圈滾道半徑；l為滾子等效接觸長度；E′為等效彈性模量。E′的表達(dá)式為

式中：E1和E2分別為滾子與滾道的彈性模量；ν1和ν2分別為滾子與滾道的泊松比。與滾道對(duì)應(yīng)的滾子等效半徑用R表示，其值取為

1.2 建立有限元模型

根據(jù)圓柱滾子軸承滾子與滾道接觸等效模型(見圖3)，以圓柱滾子軸承N306為例，建立含局部剝落故障的滾子與滾道接觸有限元二維模型。該模型的寬度和高度與圖2所示等效模型一致。其中，滾子半徑為R1，局部剝落故障的深度為d，在接觸位置對(duì)網(wǎng)格進(jìn)行加密處理，單元網(wǎng)格邊長應(yīng)不大于接觸橢圓的半軸長度[16]，單元長度取為1.5 μm，遠(yuǎn)離接觸位置的單元長度逐漸增大。約束模型底部邊緣，將載荷豎直施加在滾子中心，滾子與滾道之間接觸的目標(biāo)單元與接觸單元分別為TARGE169和CONTA172[17]。

圖3 滾子與滾道接觸有限元模型Fig.3 Finite element model for contact between roller and raceway of cylindrical roller bearing

1.3 有限元計(jì)算理論

滾子與滾道接觸有限元模型的求解過程中，有限元模型的剛度矩陣、位移矩陣和載荷矩陣的方程表達(dá)式為[18]

式中：K，u和F分別為滾子與滾道之間的接觸剛度、位移和載荷矩陣。在獲取有限元模型的節(jié)點(diǎn)位移之后，根據(jù)幾何方程和材料的應(yīng)力-應(yīng)變關(guān)系，可計(jì)算獲得有限元模型各結(jié)點(diǎn)的應(yīng)變和應(yīng)力。

為了準(zhǔn)確求解滾子與滾道接觸有限元模型中的接觸變形，必須選擇合適的接觸算法求解。有限元方法常用的接觸算法包括：罰函數(shù)法、拉格朗日乘子法以及增強(qiáng)拉格朗日算法[19]。罰函數(shù)法是將接觸非線性問題轉(zhuǎn)化為材料非線性問題，在接觸面之間利用接觸彈簧單元建立連接關(guān)系，該彈簧剛度被定義為接觸剛度。接觸體的接觸壓力與接觸變形尺寸之間的關(guān)系可表示為[17]

式中：p為接觸壓力；Kn為接觸剛度；un為接觸變形尺寸。為了找到精確的拉格朗日乘子，增強(qiáng)拉格朗日算法對(duì)罰函數(shù)進(jìn)行的一系列修正迭代，其接觸體的接觸壓力與接觸變形尺寸之間的關(guān)系可表示為[17]

式中：

式中：?為輸入容差；λi為第i次迭代的拉格朗日乘子分量。然而，拉格朗日乘子法則不需要定義接觸剛度，僅需定義顫振控制參數(shù)。理論上，拉格朗日乘子法在接觸閉合時(shí)強(qiáng)制零穿透，在發(fā)生黏著接觸時(shí)強(qiáng)制定義為“零滑移”。但是，拉格朗日乘子法增加了模型的附加自由度，需要額外的迭代來穩(wěn)定接觸條件，會(huì)增加計(jì)算成本。

1.4 等效模型驗(yàn)證

滾子與滾道之間的接觸算法采用增強(qiáng)的拉格朗日接觸算法，將有限元等效模型與數(shù)值計(jì)算方法進(jìn)行對(duì)比驗(yàn)證。對(duì)于正常鋼制圓柱滾子軸承，常用計(jì)算滾子與滾道之間的接觸變形的方法為Palmgren經(jīng)驗(yàn)公式[20]，即：

式中：δ為滾子與滾道之間的接觸變形。利用式(8)計(jì)算正常滾子與滾道之間的接觸變形，將計(jì)算得到的結(jié)果與有限元方法計(jì)算得到的結(jié)果進(jìn)行對(duì)比，結(jié)果如表1所示。

表1 有限元模型和Palmgren法計(jì)算得到的接觸變形對(duì)比Table 1 Comparison of deformations between roller and raceway obtained from FEM and Palmgren’s method

從表1可知：當(dāng)軸承徑向載荷分布范圍為100～534 N時(shí)，有限元計(jì)算得到的接觸變形與Palmgren經(jīng)驗(yàn)公式計(jì)算得到的接觸變形之間的差異小于5%，故該有限元等效模型適用于研究圓柱滾子軸承滾子與滾道之間的接觸特性。

2 局部剝落故障對(duì)軸承接觸特性的影響

文獻(xiàn)[15]的研究表明：對(duì)于圓柱滾子軸承，滾子與滾道之間接觸寬度、接觸變形和接觸應(yīng)力取決于滾子的徑向載荷；接觸寬度、接觸變形和接觸應(yīng)力還與滾子半徑、滾道半徑和等效接觸長度相關(guān)。然而，滾道表面出現(xiàn)局部剝落故障時(shí)，滾子與滾道之間的等效接觸長度將出現(xiàn)變化，從而影響滾子與滾道之間接觸寬度、接觸變形和接觸應(yīng)力。

以圓柱滾子N306 研究對(duì)象，其幾何參數(shù)如表2所示。該軸承的滾子直徑為11 mm，內(nèi)圈滾道半徑為20.25 mm，滾子等效接觸長度為11.4 mm，XP取為-1.5b～1.5b。本次計(jì)算選取25個(gè)等效接觸位置進(jìn)行研究。

表2 圓柱滾子軸承N306的幾何尺寸參數(shù)Table 2 Geometric parameters of a CRB N306

2.1 局部剝落故障對(duì)接觸應(yīng)力和接觸變形的影響

利用圖3所示的有限元模型求解滾子與局部剝落故障周圍區(qū)域之間的接觸變形、接觸寬度和接觸應(yīng)力的分布狀態(tài)。選取滾子的外載荷為19.5 N，方向定義為豎直向下。選?。?.5b≤Xp≤1.5b，b為式(1)中F取19.5 N 時(shí)對(duì)應(yīng)的接觸寬度，滾子的外載荷保持不變。圖4所示為局部剝落故障寬度分別為0.5b，0.75b和b時(shí)滾子通過局部剝落故障區(qū)域滾道的過程中滾子與滾道之間接觸變形量、接觸應(yīng)力和接觸寬度。從圖4(a)可知：當(dāng)局部剝落故障寬度為b時(shí)，隨著滾子與局部剝落故障邊緣之間距離的減小，滾子與滾道之間的接觸變形從0.563 μm 增大為0.812 μm，后又減小到0.682 μm，呈先增大后減小的趨勢。從圖4(b)可知：當(dāng)局部剝落故障寬度為b時(shí)，隨著滾子與局部剝落故障邊緣之間距離的減小，滾子與滾道之間的接觸應(yīng)力從386 MPa 增大到1 644 MPa，呈增大的趨勢。從圖4(c)可知：當(dāng)局部剝落故障寬度為b時(shí)，隨著滾子與局部剝落故障邊緣之間距離的減小，滾子與滾道之間的接觸寬度從73 μm 減小為37 μm，呈減小趨勢。仿真結(jié)果顯示，隨著局部剝落故障寬度的增加，滾子與滾道之間的接觸變形量與接觸應(yīng)力也隨之增加，而滾子與滾道之間的接觸寬度則隨之減小。

圖4 滾道接觸位置的變形、應(yīng)力和接觸寬度Fig.4 Deformation,stress and contact width of contact

圖5所示為滾子通過局部剝落故障的過程中滾子與滾道上應(yīng)力分布云圖。從圖5可知：滾子與滾道之間的最大接觸應(yīng)力所在的位置發(fā)生了變化，當(dāng)接觸點(diǎn)距離局部剝落故障邊緣的距離較大時(shí)，最大接觸應(yīng)力位置位于接觸點(diǎn)附近，在滾子趨近局部剝落故障邊緣的過程中，最大接觸應(yīng)力位置逐漸移動(dòng)到局部剝落故障邊緣根部A，當(dāng)滾子進(jìn)入局部剝落故障區(qū)間時(shí)，最大接觸應(yīng)力位置到達(dá)局部剝落故障邊緣頂部B。

2.2 局部剝落故障對(duì)接觸剛度的影響

徑向載荷與滾子-滾道之間接觸變形的關(guān)系表達(dá)式為[21]

式中：n為載荷-變形系數(shù)。為了簡化擬合過程，對(duì)式(9)兩邊分別取對(duì)數(shù)，可將指數(shù)型公式轉(zhuǎn)化為線性擬合問題，從而獲得不同接觸位置時(shí)滾子與滾道之間的接觸剛度，其表達(dá)式為

取XP=-1.5b到XP=1.5b過程中，不同載荷作用下滾子與滾道的接觸變形，利用數(shù)值擬合的方式計(jì)算滾子與滾道之間的接觸剛度。當(dāng)滾子徑向載荷從10 N 增大到25 N 時(shí)，XP從-1.5b增大到1.5b過程中不同局部剝落故障寬度下對(duì)應(yīng)的滾子與內(nèi)圈滾道之間的接觸剛度，如圖6所示。從圖6可知：當(dāng)局部剝落故障寬度為b時(shí)，滾子與滾道之間的最小接觸剛度為1.18×108N/m1/0.7231，最大接觸剛度為4.98×108N/m1/0.8957。隨著滾子與局部剝落故障中心位置距離的減小，滾子與滾道之間的接觸剛度先減小后增大，且滾子位于局部剝落故障區(qū)間時(shí)，滾子與滾道之間的接觸剛度小于滾子與正常滾道之間的接觸剛度。

3 局部剝落故障過渡區(qū)對(duì)軸承振動(dòng)特性的影響

3.1 圓柱滾子軸承的動(dòng)力學(xué)方程

圖5 滾子與滾道接觸處應(yīng)力分布Fig.5 Stress distribution between roller and raceway

圖6 接觸剛度變化規(guī)律Fig.6 Contact stiffness between roller and raceway

為研究局部剝落故障對(duì)軸承振動(dòng)特性的影響，將軸承滾子考慮為非線性接觸彈簧，引入文獻(xiàn)[22]中的簡化圓柱滾子軸承2 自由度系統(tǒng)的動(dòng)力學(xué)方程，其表達(dá)式為

式中：min為內(nèi)圈與轉(zhuǎn)軸總質(zhì)量；c為系統(tǒng)阻尼；K(θ)為滾子與內(nèi)外圈滾道的綜合接觸剛度，其數(shù)值與圖6中的一致；Fx和Fy分別為沿x和y方向的徑向載荷分量；δi為第i個(gè)滾子與內(nèi)外圈的接觸變形。

任意角度θi下的滾子與內(nèi)外圈滾道的綜合接觸變形為

式中：Hd為滾子通過軸承滾道局部剝落故障位置的時(shí)變位移激勵(lì)。時(shí)變位移激勵(lì)Hd通常用半正弦函數(shù)表示[22]，考慮到軸承局部剝落故障過渡區(qū)對(duì)時(shí)變位移激勵(lì)的影響，將其表示為

式中：Δθ=wt/rd；wt為過渡區(qū)的寬度；rd為局部剝落故障所在的套圈的滾道半徑。

3.2 局部剝落故障對(duì)軸承振動(dòng)響應(yīng)的影響

運(yùn)用定步長的四階龍格庫塔方法對(duì)方程(11)進(jìn)行求解。時(shí)間步長Δt=5 μs，方程中軸承內(nèi)圈與轉(zhuǎn)軸的總質(zhì)量min=1.2 kg，系統(tǒng)阻尼系數(shù)設(shè)定為350 N?s/m[22-23]，轉(zhuǎn)軸速度為4 000 r/min，施加在轉(zhuǎn)軸上的載荷分別為Fx=400 N 和Fy=0 N。系統(tǒng)初始位移為x=1 μm 和y=1 μm。系統(tǒng)初始速度為0 m/s。圓柱滾子軸承N306的參數(shù)如表2所示。

取外圈滾道上局部剝落故障的寬度為0.5b和b，分別計(jì)算圓柱滾子軸承的振動(dòng)響應(yīng)。其中，局部剝落故障的深度都為b，過渡區(qū)寬度取局部剝落故障寬度的1/2，缺陷相對(duì)于X軸正方向的初始角位置θd0=0°，計(jì)算結(jié)果如圖7所示。從圖7可知：局部剝落故障的寬度為0.5b時(shí)，包含局部故障過渡區(qū)的軸承X向最大加速度為35.5 m/s2，而不含局部故障過渡區(qū)的軸承X向最大加速度為20.64 m/s2；局部剝落故障的寬度為b時(shí)，包含局部故障過渡區(qū)的軸承X向最大加速度為78.43 m/s2，而不含局部故障過渡區(qū)的軸承X向最大加速度為36.35 m/s2。對(duì)比結(jié)果顯示，考慮軸承局部剝落故障過渡區(qū)時(shí)軸承振動(dòng)加速度大于不考慮軸承局部剝落故障過渡區(qū)的軸承振動(dòng)加速度。

圖8 不同剝落故障寬度對(duì)應(yīng)的軸承X方向位移的FFT變換對(duì)比Fig.8 Comparisons of FFT inX direction of bearing with different spalling fault widths

外圈滾道上局部剝落故障的寬度為0.5b和b時(shí)，軸承X方向位移的頻譜圖如圖8所示。從圖8可知：2種故障寬度下都能獲得相同的外圈通過頻率及其倍頻成分。2種故障寬度下獲得的外圈通過頻率均為314.7 Hz，與采用文獻(xiàn)[15]的理論計(jì)算方法獲得的結(jié)果一致，驗(yàn)證了仿真結(jié)果的正確性。

4 結(jié)論

1)當(dāng)滾道存在局部剝落故障時(shí)，隨著滾子與局部剝落故障中心位置距離的減小，滾子與滾道之間的接觸應(yīng)力增大，接觸寬度減小，接觸變形先增大后減小。

2)在滾子與局部剝落故障中心位置距離的減小過程中，滾子與滾道之間的最大接觸應(yīng)力位置由故障邊緣根部移動(dòng)到故障邊緣頂部。

3)隨著局部剝落故障寬度的增加，滾子與滾道之間的接觸寬度減小，接觸應(yīng)力與接觸變形量增大。

4) 滾子在通過局部剝落故障位置的過程中，滾子與局部剝落故障區(qū)域的接觸剛度小于正常接觸剛度。

5)在圓柱滾子軸承的動(dòng)力學(xué)系統(tǒng)中，考慮軸承局部剝落故障過渡區(qū)時(shí)對(duì)應(yīng)的軸承振動(dòng)加速度大于不考慮局部故障過渡區(qū)的軸承振動(dòng)加速度。