滾動軸承復合故障振動建模及Lempel-Ziv 復雜度評價

朱永生，袁 幸，張優云，朱川峰，周 智

(西安交通大學 潤滑理論及軸承研究所，西安 710049)

滾動軸承被廣泛應用于機械工業，是基礎件也是關鍵核心件。由于其通常工作在高速、高溫、重載的惡劣環境下，滾動軸承疲勞導致的失效往往在整機故障中占大部分比率。因此故障診斷對于滾動軸承的安全、連續運行和預知更新至關重要［1］。故障建模研究故障狀態下動力學參數與響應征兆的內在聯系，籍此獲取服役正常與異常的證據信息［2］，其目的在于根據機械設備性能退化響應特點來識別潛在故障和構建故障庫。

滾動軸承動力學建模一直為學者所重視，Fukata［3］較早分析了滾動軸承的變柔性振動。研究結論顯示，滾動軸承支承下的不平衡轉子具有兩個激勵源，即來自不平衡的轉頻激勵和來自軸承總體剛度周期變化的內部參數激勵。Tiwari等［4］采用兩自由度模型研究軸承游隙造成的轉子系統的非線性行為和穩定性。然后Lioulios等［5］探討了速度波動對軸承非線性動力學的影響。Patil等［6］模擬內、外圈局部故障，研究單點缺陷所激勵的時、頻域特性。Liew等［7］進一步研究非線性Hertzian接觸、游隙、轉速等因素對動態行為的影響機制，并有效模擬了軸承元件的振動形態。Sopanen等［8－9］提出一種軸承模型用于分析局部及分布式故障，模型考慮了非線性Hertzian接觸，但是缺乏外圈－軸承座的耦合效應。澳大利亞UNSW大學Sawalhi等［10］采用附加自由度的方法建立5DOF模型，從而得到了相對精確的振動響應。一些新的建模手段也被引入，Mohsen等［11］用鍵合圖功率流構建滾動軸承振動模型，分析了軸承元件的動態接觸行為。機械設備出現故障時，故障的演變存在一個由輕微到嚴重、由單一到復合的過程［12］。對于長期在變載工況下服役的機床軸承、機車軸承尤為顯著，監測到的故障往往以多故障耦合的形態出現，研究滾動軸承復合故障建模有著重要工程意義［13］。然而，現有的動力學模型均以單點缺陷為研究對象，復合故障建模的研究鮮見報道。

為能夠有效預示多耦合故障誘發的振動特征，基于軸承運動學與Hertzian接觸理論，建立一種計及內圈－外圈－軸承座的6DOF振動模型。通過漸進過程描述缺陷區域，以深溝球軸承SKF6205－2RS為實例計算振動響應。與實驗數據對比表明，利用該方法可深入揭示軸承局部缺陷動力特性與宏觀動態性能之間的映射關系。進一步著重研究復合故障建模方法及其響應特點。在此基礎上引入非線性動力系統中復雜性刻畫方法——Lempel-Ziv復雜度指標用于評價不同故障狀態動力系統的變化。研究表明，這種特征量是表征不同故障信息的有效參數。

1 滾動軸承接觸建模

1.1 非線性彈性恢復力

圖1 滾動軸承運動示意圖Fig.1 Schematic diagram of rolling bearing motion

如圖1，第i個滾動體－套圈彈性趨近量δi是套圈的徑向位移(xs，ys)，滾動體位置角θi和游隙γ的函數，由運動關系可知:滾動體位置角ωct，其中:ωc為公轉速度，n為滾動體個數。設ω為相配軸的轉速，Db為滾動體直徑，Dp為軸承節圓直徑，ωc設軸承座運動坐標為(xp，yp)，則內圈

－軸承座的綜合彈性趨近量就為:需要說明的是，由于載荷分布區具有角度范圍，δi只能取正值。

1.2 載荷 －位移關系［14］

根據Hertzian理論，點接觸載荷－位移關系為:

式中:δ為彈性趨近量，K為Hertzian接觸剛度。為了描述兩個嚙合表面之間的接觸狀況，定義Kinner、Kouter、∑ρ分別為滾動體與內、外圈的接觸剛度與曲率和，則:

其中:E為彈性模量，ν為泊松比，μ*為變形系數。這樣，一旦軸承幾何尺寸確定，Kinner、Kouter就容易得到。進而求出軸承整體等效接觸剛度為:

滾動軸承運動副的等效約束力可由每個滾動體所受彈性接觸力求得，故總的非線性彈性恢復力可表示為:

1.3 接觸阻尼

滾動軸承接觸阻尼隨平均負荷變化，其值很小，Kramer通過動態測試方法給出了經驗計算公式［14］

2 振動系統模型

滾動軸承振動性能是其內部運動關系在外部載荷作用下與傳遞路徑的綜合作用。內圈－外圈－軸承座為非線性強耦合關系，圖2表達了軸承接觸界面的多自由度振動形態。其中，c0、k0、cp、kp分別為外圈、軸承座的阻尼和剛度。基于多自由度振動理論，構建多體系統環境下的軸承精細化模型如圖2所示。由于滾動軸承接觸問題的復雜性及各種非線性因素的影響，無疑會增加系統建模及求解的難度。鑒于此，本研究做以下基本假設。① 忽略滾動體滑移現象。② 忽略制造誤差，如表面波紋度的影響。③ 假設滾動體與內、外圈之間只存在局部彈性接觸變形。由第一節的分析結合圖2～3，建立振動方程如下:

圖2 滾動軸承振動示意圖Fig.2 Schematic diagram of rolling bearing vibration

式中:ms為內圈(相配軸)質量，mo、mp分別為外圈和軸承座質量，e為不平衡距離，fn為外載荷。

3 局部故障嵌入

滾動軸承的失效大部分表現為局部缺陷，例如疲勞剝落、凹痕、點蝕等。缺陷的存在會形成一個有限的區域，具有長度l，如圖3所示。當滾動體運行至缺陷位置時，彈性趨近量瞬間增加，設d為相應的實際增加值，即式(2)中的游隙γ增加d，這將導致沖擊產生和能量突變。在外圈固定內圈旋轉的情況下，外圈缺陷所對應的中心角α=2arcsin(l/2R)，R為外圈半徑。設缺陷區域相對初始滾動體的位置角為，則滾動體通過外圈缺陷區域的漸進過程可描述為:

圖3 表面缺陷區域Fig.3 The surface defects region

滾動體出現故障時，缺陷區域以轉速ωb自轉一周與內、外圈分別接觸一次，每一次接觸彈性趨近量均增加d，此時判別函數由兩部分組成。設滾動體缺陷區域與外圈初始相對角位置為，則缺陷區域通過外圈的漸進過程為:

同理，滾動體缺陷區域通過內圈的漸進過程為:

復合故障指同時發生兩個或多個互相關聯、交叉影響的故障。以內、外圈復合故障為例，動態建模表現為式(10)～(11)的級聯，即同時運行內、外圈的漸進過程。

4 計算分析及驗證

式(9)所代表的狀態方程組非線性強，難以解析，本文采用數值方法獲取系統響應。計算方法為定步長Runge–Kutta法，步長設為1/12 000 s，初始值為10－6m。算例采用美國凱斯西儲大學(CWRU)軸承研究中心的深溝球軸承(SKF6205－2RS)，轉速1 772 r/min，采樣頻率12 000 Hz。故障軸承采用電火花加工制成，內、外圈及滾動體的故障頻率依次為5.415 2、3.584 8、4.713 5乘以工頻，CWRU提供了詳細設計參數可供參考［15］。外圈發生故障時頻譜結構相對簡單，主頻為故障特征頻率，如圖4所示。可以看出，仿真信號包絡譜體現了不平衡激勵的轉頻fr=29.66 Hz，內圈故障頻率fo=106 Hz及其諧波2fo、3fo、4fo成分，這與實驗信號非常吻合。與外圈故障的區別在于內圈缺陷連同轉軸一起運行，故障頻率必然受到工頻的調制，因此包絡譜的主頻以 fi=159.9 Hz、fi+fr=189.5 Hz、fi－ fr=130.4 Hz及其組合諧波成分 2fi+fr、2fi－ fr、3fi+fr、3fr－ fr等為主。圖5的仿真信號包絡譜成功突出了主頻成分，并且頻率預測達到了很高的精度。仿真和實驗數據的對比表明，振動模型對顯著頻率成分實現了良好的預測。為驗證模型的魯棒性，考察轉速為1 797 r/min時滾動體故障的振動響應。從圖6中清晰看出滾動體缺陷激勵引起的振動在包絡譜中表現為保持架轉頻fc=11.9 Hz對故障頻率fb=141.2 Hz的調制。這是由于滾動體上缺陷位置隨保持架旋轉不斷改變，缺陷引起的沖擊響應幅值必然受到保持架旋轉頻率調制。以上分析充分體現了本文建模方法的可靠性。

圖4 外圈故障仿真與實驗信號包絡譜Fig.4 Simulation and experimental envelope spectrum of outer race defect

圖5 內圈故障仿真與實驗信號包絡譜Fig.5 Simulation and experimental envelope spectrum of inner race defect

設置相同的缺陷長度l=0.5 mm和深度d=0.1 mm，轉速1 772 r/min時內、外圈復合故障振動響應如圖7所示。滾動軸承復合故障的信號特征比單一故障相對復雜。從剛度損失的角度可以解釋為復合故障導致軸承接觸剛度發生更大變化，在不同缺陷位置引入局部柔度，使得動力學性能改變。當內、外圈復合故障出現時，局部柔度同時增加，激勵系統出現兩種故障頻率fo和fi。包絡譜低頻部分fo及其諧波比較明顯，fi的幅值相對較小。可以解釋為復合故障由單點故障呈矢量疊加方式構成，以內、外圈復合故障為例，頻域總幅值為:

圖6 滾動體故障仿真與實驗信號包絡譜Fig.6 Simulation and experimental envelope spectrum of rolling element defect

圖7 內、外圈復合故障信號及包絡譜Fig.7 Vibration response of inner and outer race compound defects

圖8 外圈、滾動體復合故障信號及包絡譜Fig.8 Vibration response of outer race and rolling element compound defects

圖9 內圈、滾動體復合故障信號及包絡譜Fig.9 Vibration response of inner race and rolling element compound defects

圖10為內圈、外圈、滾動體三種故障復合時的振動響應，包絡譜清晰顯示了三種特征頻率的2倍頻，fb=139.2 Hz相對微弱，并且各頻率的幅值均不同。計算結果充分表明滾動軸承復合故障具有相互耦合、相互制約和強度不平衡性，也說明這些故障的產生具有內在確定性機制。同樣顯示，包絡譜所提取的故障頻率倍頻分量可作為識別復合故障的有效依據。

5 復合故障Lempel-Ziv復雜度評價

滾動軸承復合故障往往使其頻率組成及各成分幅、頻特性復雜化，從而使時域波形發生不同程度的畸變，這種畸變也是軸承故障在時域信號中的直接反映。事實上，信號的波形、頻譜結構發生變化時，信號的復雜程度也發生了變化。人們開始探索采用一些指標定量地描述軸承振動信號的復雜性［2］。姜建東等［2］用Lempel-Ziv復雜性測度定量評估大機組運行狀態。Hong等［16］研究了Lempel-Ziv復雜度指標用于評估故障軸承損傷狀態。實踐表明，這一指標是衡量有限時間序列復雜程度的高效工具。

圖10 內圈、外圈、滾動體復合故障信號及包絡譜Fig.10 Vibration response of inner race，out race and rolling element compound defects

序列 SN={s1，s2，…，sN}的 Lempel-Ziv 復雜度可按CN(r)(r≤N)經過 N 次循環得到，計算過程如下［2、16］:

(1)初始化 Sv，0={}，Q0={}，CN(0)=0，r=1。令 Qr={Qr－1Sr}，由于 Qr不屬于 Sv，r－1，則 CN(r)=CN(r－1)+1，Qr={}，r=r+1。

(2)令 Qr={Qr－1Sr}，判斷 Qr是否屬于 Sv，r－1=Sv，r－1SrSr－1;若屬于，則 CN(r)=CN(r－1)，r=r+1，重復過程(2)。

(3)若不屬于，則 CN(r)=CN(r－1)+1，Qr={}，r=r+1，重復過程(2)。

例如序列{0000…}的分段形式為{0·000…}，CN=2。符號序列{010101…}的分段形式為{0·1·0101…}，CN就為3。出于對數據嚴格的可比性考慮，Lempel和Ziv進一步提出了一種歸一化基準公式，將CN值界定在［0 1］之間，詳細算法為:

其中:k為SN中元素的個數，對二進制序列k=2。根據已有的經驗，一般N大于3 600可保證上式成立［16］。根據上述Lempel-Ziv復雜度算法和文獻［2，16］可知，一個序列的復雜度越大，說明它的添加操作越多，序列的內涵分量也越多，描述給定符號序列所需最少的、互不相同的分段形式越多，給定的序列周期性越弱，出現新內涵模式的速度也就越快。Lempel-Ziv復雜度的物理意義在于它反映了一個時間序列隨著長度的增加出現新模式的速率。復雜度越大，說明數據在窗口長度時期內隨時間出現的新變化越多，發生新變化的速度越快，表明這一時期的數據變化是無序而復雜的。反之復雜度越小，則說明發生新變化的速度越慢，數據變化是規則的，周期性越強。因此，振動信號的Lempel-Ziv復雜度指標能夠客觀反映出系統狀態發生變化的情況。

為進一步完善軸承故障評估方法，本文借助上述Lempel-Ziv復雜度指標的建立思想，并將之推廣到復合故障診斷領域。同樣設置缺陷長度l=0.5 mm和深度d=0.1 mm的微弱故障進行數值計算。圖11從左往右依次對應正常、外圈、內圈、滾動體、內－外圈復合、外圈－滾動體復合、內圈－滾動體復合、內外圈－滾動體復合故障的Lempel-Ziv復雜度變化趨勢，表1為其具體數值。

圖11 不同狀態的Lempel-Ziv數值變化趨勢Fig.11 Lempel– Ziv complexity vs.different fault

表1 不同狀態的Lempel-Ziv復雜度Tab.1 Numerical value of Lempel-Ziv complexity

從結果中可以看出，正常狀態的振動最為簡單。滾動體故障雖然激勵沖擊能量不大，但是滾動軸承中滾動體的運動比較復雜，既有公轉又有自轉。在其運動過程中，故障區域既要與內圈接觸又要與外圈接觸，漸進過程由兩部分組成(式(12)～(13))。不斷變化的接觸行為會引起強烈的調制作用，反映在頻譜中會出現越來越多的譜峰及其諧波，波形中的成分越來越雜亂，振動的周期性減弱，非線性增強，相對于內、外圈更加傾向于較高的Lempel-Ziv值。復合故障破壞了單一故障的周期沖擊振動結構，故障沖擊加劇，直觀上可以看出，信號的復雜度增加了。不同形式的復合故障其響應的復雜度也不同，這是因為不同的復合故障形成機制不同，導致信號的沖擊振動結構有所差異。內外圈－滾動體三種缺陷復合故障的振動響應與兩種缺陷復合故障的Lempel-Ziv指標差別不太大，體現了故障之間既相互干擾又相互制約。對復合故障的復雜度評價與我們的工程經驗相吻合，再次證明了本文故障建模方法的可靠性。由上述分析可知，振動信號的復雜度較為敏感地反映了軸承狀態演變過程中動力系統的變化歷程，是一種有效的特征量。實際應用中整個信號的Lempel-Ziv指標單獨作為診斷特征量稍微少了一些。基于波形匹配思想，兼備優良性質的多小波擁有多個尺度函數和小波函數，可以很好地匹配復合故障的多種特征波形［13］。利用多小波變換消除噪聲和分解信號，從而構造不同尺度下的Lempel-Ziv復雜度向量，可以同時完成單點及復合故障的監測診斷。此方法的實驗研究加之工程應用，亦是下一步的研究課題，將另文給出。

6 結論

(1)建立一種計及內圈－外圈－軸承座的滾動軸承6DOF振動模型。從軸承幾何學、運動學的角度推導了局部故障的漸進過程，通過與實驗數據的比較表明，該方法能有效預測振動頻譜的主要成分，為后續復合故障的研究奠定了基礎。

(2)進一步研究了復合故障的建模方法和振動形態。研究結果顯示，復合故障具有相互耦合性、相互制約性和強度不平衡性。包絡譜或共振解調譜可同時解耦復合故障特征頻率及其倍頻分量。一定條件下，包絡譜的倍頻分量可以作為診斷復合故障的依據。

(3)滾動軸承局部故障的出現相當于局部柔度的引入，進而在不同漸進過程瞬間改變系統動力學特性。不同故障信號的Lempel-Ziv復雜度指標反映了不同故障的動力學形成機制。軸承振動信號的復雜度能有效地評價軸承運行狀態，靈敏地反映其狀態的變遷過程。

［1］RandallR B， Antoni J. Rolling element bearing diagnostics—A tutorial［J］.Mechanical Systems and Signal Processing，2011，25(2):485－520.

［2］機械監測診斷中的理論與方法:屈梁生論文集［M］.西安:西安交通大學出版社，2009.

［3］Fukata S.On the radial vibration of ball bearings(computer simulation)［J］.Bulletin of the JSME，1985，28:899－904.

［4］Tiwari M，Gupta K.Dynamic response of an unbalanced rotor supported on ball bearings［J］.Journal of Sound and Vibration，2000，238(5):757－779.

［5］Lioulios A N，Antoniadis I A.Effect of rotational speed fluctuations on the dynamic behaviour of rolling element bearings with radial clearances［J］.International Journal of Mechanical Sciences，International Journal of Mechanical Sciences，2006，48(8):809 －829.

［6］Patil M S，Mathew J.A theoretical model to predict the effect of localized defect on vibrations associated with ball bearing［J］.International Journal of Mechanical Sciences，2010，52(9):1193－1201.

［7］Liew A，Feng N，Hahn E.Transient rotor dynamic modelling of rolling element bearing systems［J］.Transactions of the ASME，Journal of Engineering for Gas Turbines and Power，2002，124(4):984 －991.

［8］Sopanen J，Mikkola A.Dynamic model of a deep-groove ball bearing including localized and distributed defects.Part 1:theory，Proceedings of the Institution of Mechanical Engineers［J］.Part K:Journal of Multi-body Dynamics，2003(217):201－211.

［9］Sopanen J，Mikkola A.Dynamic model of a deep-groove ball bearing including localized and distributed defects.Part 2:implementation and results，Proceedings of the Institution of Mechanical Engineers［J］.Part K:Journal of Multi-body Dynamics，2003(217):213 －223.

［10］Sawalhi N，Randall R B.Simulating gear and bearing interactions in the presence of faults:Part I.The combined gear bearing dynamic model and the simulation of localized bearing faults ［J］. MechanicalSystems and Signal Processing，2008，22(8):1924－1951.

［11］Mohsen N，Michael D.Dynamic modeling of rolling element bearings with surface contact defects using bond graphs［J］.ASME，Journal of Tribology，2011，133(1):1－12.

［12］何正嘉，訾艷陽，陳雪峰.機械故障預示中的若干科學問題［J］.振動與沖擊，2010，29(S):144 －147.HE Zheng-jiia， Z1 Yan-yang，CHEN Xue-feng. Some scientific issues in mechanical fault prognosis［J］.Journal of Vibration and Shock，2010，29(S):144 －147.

［13］王曉冬.基于多小波變換的機械設備早期故障與復合故障診斷研究［D］.西安:西安交通大學，2009.

［14］Harris T A.Rolling bearing analysis［M］.4th ed.New York:John Wiley and Sons，2001.

［15］The Case Western Reserve University Bearing Data Center Website.Bearing data center seeded fault test data［EB/OL］.http:http://www.eecs/cwru/edu/laboratory/bearing.

［16］Hong H B，Liang M.Fault severity assessment for rolling element bearingsusing the Lempel-Ziv complexity and continuous wavelet transform ［J］.Journal of Sound and Vibration，2009，320(1 －2):452 468.