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

(1.重慶大學 機械傳動國家重點實驗室，重慶，400044；2.重慶大學機械工程學院，重慶，400044)

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

1 模型描述

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

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

圖1 圓柱滾子與滾道之間的位置關系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為等效滾子中心位置的最大赫茲接觸應力。等效模型的長度為10b，寬度為7b，其中，b為用赫茲接觸模型計算得到的軸承接觸半寬。接觸半寬b表示為[15]

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

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

1.2 建立有限元模型

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

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

1.3 有限元計算理論

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

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

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

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

式中：

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

1.4 等效模型驗證

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

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

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

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

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

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

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

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

2.1 局部剝落故障對接觸應力和接觸變形的影響

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

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

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

2.2 局部剝落故障對接觸剛度的影響

徑向載荷與滾子-滾道之間接觸變形的關系表達式為[21]

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

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

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

3.1 圓柱滾子軸承的動力學方程

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

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

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

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

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

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

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

3.2 局部剝落故障對軸承振動響應的影響

運用定步長的四階龍格庫塔方法對方程(11)進行求解。時間步長Δt=5 μs，方程中軸承內(nèi)圈與轉軸的總質(zhì)量min=1.2 kg，系統(tǒng)阻尼系數(shù)設定為350 N?s/m[22-23]，轉軸速度為4 000 r/min，施加在轉軸上的載荷分別為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，分別計算圓柱滾子軸承的振動響應。其中，局部剝落故障的深度都為b，過渡區(qū)寬度取局部剝落故障寬度的1/2，缺陷相對于X軸正方向的初始角位置θd0=0°，計算結果如圖7所示。從圖7可知：局部剝落故障的寬度為0.5b時，包含局部故障過渡區(qū)的軸承X向最大加速度為35.5 m/s2，而不含局部故障過渡區(qū)的軸承X向最大加速度為20.64 m/s2；局部剝落故障的寬度為b時，包含局部故障過渡區(qū)的軸承X向最大加速度為78.43 m/s2，而不含局部故障過渡區(qū)的軸承X向最大加速度為36.35 m/s2。對比結果顯示，考慮軸承局部剝落故障過渡區(qū)時軸承振動加速度大于不考慮軸承局部剝落故障過渡區(qū)的軸承振動加速度。

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

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

4 結論

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

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

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

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

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