韓東穎，李 庚，時培明

(1.燕山大學 車輛與能源學院，秦皇島 066004;2.燕山大學 河北省測試計量技術及儀器重點實驗室，秦皇島 066004)

旋轉(zhuǎn)機械轉(zhuǎn)子系統(tǒng)中轉(zhuǎn)軸裂紋、動靜碰摩和支撐部件松動為常見故障，是導致機械系統(tǒng)失效甚至造成嚴重事故的主要原因。目前國內(nèi)外學者［1－3］對裂紋、碰摩、松動故障單獨存在的系統(tǒng)研究較多，提出了許多有效的診斷方法。

但在實際轉(zhuǎn)子中，常常出現(xiàn)兩種甚至三種故障同時存在的情況。這種耦合故障轉(zhuǎn)子的動力學行為較單一故障轉(zhuǎn)子更加復雜，而且相互影響，不容易診斷。文獻［4］采用有限元方法分析了裂紋碰摩耦合故障轉(zhuǎn)子系統(tǒng)的響應，在嚴重碰摩故障信號中診斷出了裂紋故障。文獻［5］建立了帶有裂紋－支承松動耦合故障的雙跨彈性轉(zhuǎn)子系統(tǒng)動力學模型，利用數(shù)值仿真對故障非線性響應進行研究。文獻［6］通過分岔圖、Poincaré映射圖和幅值譜圖，分析了碰摩－松動故障機組軸系隨轉(zhuǎn)頻比和轉(zhuǎn)子質(zhì)量偏心變化的非線性動態(tài)響應。

EMD是近年來發(fā)展起來的處理非平穩(wěn)、非線性信號的時頻分析方法［7］。該方法克服了傳統(tǒng)時頻分析方法中的不足，具有很強的自適應性，并在機械故障診斷領域得到了廣泛應用［8－10］。分形的本意為不規(guī)則的、分數(shù)的、支離破碎的物體，可以理解為局部和整體在某個方面有相似性［11］。分形維數(shù)是度量分形的重要指標，有許多種分形維數(shù)可用來刻畫信號的復雜度［12］，如Hausdorff維、盒子維、容量維、信息維、關聯(lián)維數(shù)等，其中盒維數(shù)由于其計算簡單、高效而應用最為廣泛［13－15］。

針對耦合故障信號復雜，具有強非線性的特點，本文提出一種基于EMD和分形盒維數(shù)的耦合故障診斷方法。該方法先利用EMD將故障信號分解，然后對有效IMF進行盒維數(shù)計算，得到耦合故障信號的故障特征向量，判斷故障狀態(tài)，達到故障診斷的目的。

1 系統(tǒng)動力學模型和運動微分方程

如圖1所示，建立含有裂紋－碰摩－松動耦合故障的剛性支承轉(zhuǎn)子－軸承系統(tǒng)模型，轉(zhuǎn)子圓盤與軸承之間為無質(zhì)量的彈性軸。模型左端發(fā)生松動，軸承座與基礎之間的松動最大間隙為ξ1。轉(zhuǎn)子圓盤左側有一弓形橫向裂紋，其深度為a，并且圓盤受到碰摩故障影響，其碰摩函數(shù)為P。圖中O1為軸瓦幾何中心，O2為轉(zhuǎn)子幾何中心，O3為轉(zhuǎn)子質(zhì)心，kc為定子剛度，k為彈性軸剛度，m1為兩端軸承處的轉(zhuǎn)子集中質(zhì)量，m2為轉(zhuǎn)子圓盤的等效集中質(zhì)量，m3為軸承支座的等效集中質(zhì)量。模型還考慮了左端滑動軸承作用在轉(zhuǎn)軸上的非線性油膜力，為別為 Fx、Fy。

圖1 含裂紋碰摩松動故障的轉(zhuǎn)子－軸承系統(tǒng)模型Fig.1 Rotor-bearing system consider in crack rub-impact and pedestal looseness

設轉(zhuǎn)子右端的徑向位移為x1，y1;轉(zhuǎn)子圓盤的徑向位移為x2，y2;松動端軸心位移為x3，y3;軸承支座在豎直方向位移為y4，則具有裂紋碰摩松動耦合故障的轉(zhuǎn)子－軸承系統(tǒng)運動微分方程為:

圖2 裂紋截面模型Fig.2 Crack section model

式中:u為轉(zhuǎn)子的偏心量，c1為轉(zhuǎn)子在軸承處的阻尼系數(shù)，c2為轉(zhuǎn)子圓盤的阻尼系數(shù)，cs為支座松動阻尼系數(shù)，ks為支承剛度。ω為轉(zhuǎn)子轉(zhuǎn)速，g為重力加速度，ε、δ為僅與裂紋深度a有關的相對剛度參數(shù)，由式(3)、(4)所得。如圖2所示，F(xiàn)(Ψ)為裂紋開閉函數(shù)，本文采用余弦波模型來表示裂紋開閉過程，粗略地考慮裂紋半開半閉的過渡過程，忽略了裂紋的全閉和全開是一個持續(xù)過程。

余弦波模型的數(shù)學表達式為:

式中:φ0為初相位，β為裂紋方向與偏心之間的夾角，x，y為轉(zhuǎn)子初始位置。

式(1)中的轉(zhuǎn)定子碰摩函數(shù)P為:

式中:Pn為法向碰摩力，Pτ為切向摩擦力，ξ為靜止時轉(zhuǎn)子與定子之間的間隙，e為轉(zhuǎn)子軸心位移，由圖 3所示，當位移e大于等于間隙 ξ時系統(tǒng)發(fā)生碰摩。

碰撞發(fā)生時(e≥ξ)，其正向碰摩力與切向摩擦力分解在x－y坐標系可以表示為:

圖3 碰摩力模型Fig.3 Rub － impact model

式中:kc為定子徑向剛度，f為不考慮速度影響時的摩擦系數(shù)，b為速度影響系數(shù)為轉(zhuǎn)子靜件間的相對滑動速度為轉(zhuǎn)子的徑向位移，其中x和y分別對應式(1)中的x2和y2。

松動故障等效成剛度和阻尼的變化，支承間隙系統(tǒng)在位移條件下ks、cs為分段性，其表達式為:

式(1)中油膜力沿x和y兩個方向的分量為:

式中:μ為潤滑油粘度，c為軸承徑向間隙，R為軸承半徑，L為軸承長度。

2 經(jīng)驗模式分解

經(jīng)驗模式分解EMD是一種自適應分解方法，可以把復雜的信號分解為有限個IMF分量。IMF信號一般滿足兩個條件:① 從全局特性上看，極值點數(shù)必須和過零點數(shù)一致或者至多相差一個;② 在某個局部點，極大值包絡和極小值包絡在該點的算術平均值是零，即兩條包絡線關于時間軸對稱。

可以把任何信號x(t)按下面步驟分解:

(1)用三次樣條線將所有的局部極大值點連接起來形成上包絡線。

(2)用三次樣條線將所有的局部極小值點連接起來形成下包絡線。

(3)上下包絡線的平均值記為m1，求出:

理想情況，如果h1是一個IMF，則h1就是x(t)的第一分量。

(4)如果h1不滿足IMF的條件，把h1作為原始據(jù)，重復(1)、(2)、(3)，得到上下包絡線的平均值m1再判斷h11=h1－m11是否滿足IMF的條件，如不滿足，重復循環(huán) k次，得到 h1k=h1(k－1)－m1k，使得 h1k滿足IMF條件。記c1=h1k，則c1為信號x(t)的第一個滿足IMF條件的分量。

(5)將c1從x(t)中分離出來，得到:

將r1作為原始數(shù)據(jù)重復以上過程，得到x(t)的第二個滿足IMF條件的分量c2，重復循環(huán)n次，得到n個滿足IMF條件的分量。這樣就有:

當rn成為一個單調(diào)函數(shù)不能再從中提取滿足IMF條件的分量時，循環(huán)結束。這樣由式(10)和(11)得到:

3 分形盒維數(shù)及其計算

3.1 分形盒維數(shù)

設X是Rn的非空有界子集，記N表示最大直徑為ε且能覆蓋X集合的最少個數(shù)，則X的盒維數(shù)定義為:

3.2 盒維數(shù)計算方法

設離散信號y(i)?Y，Y是n維歐氏空間Rn上的閉集。用盡可能細的ε網(wǎng)格劃分Rn，Nε是集合Y的網(wǎng)格計數(shù)。由于式(13)中的極限無法按定義求出，所以在計算時需要采用近似的方法。以ε網(wǎng)格作為基準，逐步放大到kε網(wǎng)絡，其中 k∈Z+。這樣，令 Nkε為離散空間上的集合Y的網(wǎng)格計數(shù)，則由式(14)和式(15)可以計算得到。

式中:j=1，2，…，N/k;N 為采樣點數(shù)。k=1，2，…，M，M＜N。

網(wǎng)絡計數(shù)Nkε為:

式中:Nkε＞1。

在lgkε－lgNkε圖中確定線性較好的一段為無標度區(qū)，設無標度區(qū)的起點和終點分別為k1、k2，則:

最后，用最小二乘法確定該直線的斜率:

盒維數(shù)dB為:

4 數(shù)值仿真和故障診斷

由方程(1)可以看出，含有裂紋碰摩松動耦合故障的轉(zhuǎn)子系統(tǒng)是一個有復雜外激勵的非線性系統(tǒng)。目前分析這種系統(tǒng)最有效的方法就是數(shù)值仿真，本文采用變步長四階龍格－庫塔法對方程(1)進行數(shù)值求解，系統(tǒng)參數(shù)選取如下:m1=4 kg，m2=32.5 kg，m3=50 kg，R=0.025 m，L=0.012 m，c=0.11 mm，a=0.015 m，μ=0.018 Pa·s，c1=1 050 N·s/m，c2=2 100 N·s/m，cs1=350 N·s/m，cs2=500 N·s/m ，f=0.1，b=0.01，k=2.5 ×107N/m，u=0.05 mm，kc=4.5 ×107N/m，ks1=2.5 ×107N/m，ks2=2.5 ×109N/m，ω =789.3 rad/s，ξ=0.5 mm，ξ1=1 mm，φ0=0，β =0。

考慮到三種耦合故障模型結構特點，轉(zhuǎn)子左端的Y方向位移最能體現(xiàn)耦合故障特性。圖4為數(shù)值解得到的轉(zhuǎn)子左端徑向位移y3的時域圖。由圖4可以看出由于機械故障的影響，y3的振動有很強的非線性。

圖4 y3時域圖Fig.4 Time domain chart of y3

EMD方法用于處理非線性、非平穩(wěn)信號有良好的效果。如圖5所示，把由龍格－庫塔法解出的y3振動信號經(jīng)過EMD方法分解，得到含有故障特征的IMF。由于EMD方法本身原因產(chǎn)生虛假模態(tài)，故只給出IMF1～IMF4。

圖5 y3的EMD處理結果Fig.5 Empirical mode decomposition(EMD)results of y3

由分解得到的IMF可以看出信號每隔一段時間就會有頻率和幅值上的明顯變化，且這種變化具有周期性。根據(jù)已知的單一故障特征可以判斷這是由于松動或者碰摩故障引起的。要對故障進一步診斷，傳統(tǒng)方法是對每個IMF進行希爾伯特變換，求出邊界譜，通過邊界譜的倍頻關系進行故障診斷。

圖6是轉(zhuǎn)子左端的Y方向位移y3的邊界譜圖。由圖可以看出振動主要是由1/2倍頻組成的，其次是1倍頻，2倍頻等高頻成分則不是主要振動頻率。為了更好地進行故障診斷，需要把裂紋、松動、碰摩兩兩耦合情況下振動信號的邊界譜與三種故障同時存在的情況作比較。

把三種耦合故障模型進行簡化并分析，可以得到兩種耦合故障的邊界譜(如圖7～9)，由于篇幅有限不做詳細說明。圖7為碰摩松動故障信號的邊界譜，由圖可知信號振動主要由1/2倍頻構成，和三種耦合故障情況相似，但1倍頻成分不如圖6明顯，說明1/2倍頻振動主要是由于碰摩或松動引起的。圖8為裂紋松動故障信號的邊界譜。通過與圖6對比可以明顯看出裂紋松動故障信號不僅含有1/2倍頻、1倍頻成分，也含有3/2倍頻、2倍頻等高頻成分，由此可以診斷出裂紋故障是引起高頻振動的主要原因。圖9為裂紋碰摩故障信號的邊界譜，從圖中可以清晰地看出振動信號同時具有高頻和低頻振動。通過上面的分析可知高頻是有裂紋故障引起的，低頻是有碰摩故障引起的。

通過對四幅邊界譜圖的分析，可以得出裂紋松動碰摩三種故障同時存在時與其兩兩存在時故障特征的區(qū)別。但由于圖像抽象、頻率繁雜，不能直觀的把這四種情況區(qū)分開，所以本文在EMD處理的基礎上結合分形盒維數(shù)的方法，求出每種情況下IMF1～IMF4的盒維數(shù)，準確、清楚地反映出四種耦合故障特征。

圖6 y3邊界譜Fig.6 Boundary spectrum of y3

圖7 碰摩松動故障信號的邊界譜Fig.7 Boundary spectrum ofrub-impact and pedestal looseness signal

圖8 裂紋松動故障信號的邊界譜Fig.8 Boundary spectrum of crack and pedestal looseness signal

圖9 裂紋碰摩故障信號的邊界譜Fig.9 Boundary spectrum of crack and rub－impact signal

由表1可知，EMD處理后的故障信號分形盒維數(shù)有明顯的差異，其中IMF1到IMF4的盒維數(shù)相差最大的為裂紋碰摩故障，最小的為裂紋松動故障。三種故障耦合的信號盒維數(shù)值與碰摩松動故障和裂紋碰摩故障差異較大，雖與裂紋松動故障差異相對較小，但也能夠清晰分辨，所以分形盒維數(shù)可以作為故障的特征向量，對耦合故障進行診斷。

表1 四種故障信號經(jīng)EMD處理后的盒維數(shù)Tab.1 Processed box dimension by EMD of four kinds of fault signal

5 結論

(1)提出一種基于EMD和分形盒維數(shù)的旋轉(zhuǎn)機械耦合故障診斷方法。該方法先用EMD對故障信號進行分解，得到IMF，再根據(jù)分形計算方法求出IMF的盒維數(shù)，通過對盒維數(shù)的分析比較進行故障診斷。

(2)建立含有裂紋－碰摩－松動故障耦合的轉(zhuǎn)子－軸承系統(tǒng)動力學模型。在考慮單個故障特性的基礎上得出模型運動微分方程，并用龍格－庫塔法解出含有耦合故障特征的振動信號。

(3)把新方法應用到所建故障模型的診斷中，得到裂紋、碰摩、松動三種故障耦合故障特征向量。通過與傳統(tǒng)診斷方法比較，證明該方法能夠有效、準確地對耦合故障進行診斷。

［1］Hu N Q，Wen X S.The application of Duffing oscillator in characteristic signal detection of early fault［J］.Journal of Sound and Vibration，2003，268:917 －931.

［2］Sekhar A S，Mohanty A R，Prabhakar R.Vibrations of cracked rotor system:transverse crack versus slant crack［J］.Journal of Sound and Vibration，2005，279:1203－1217.

［3］馬 輝，孫 偉，任朝暉，等.多盤懸臂轉(zhuǎn)子系統(tǒng)支座松動故障研究［J］.航空動力學報，2009，24(7):1512－1517.MA Hui，SUN Wei， REN Zhao-hui et al. Research on pedestal looseness fault of multi-disc over-hung rotor system［J］.Journal of Aerospace Power，2009，24(7):1512－1517.

［4］姚紅良，李 鶴，李小彭，等.裂紋碰摩耦合故障轉(zhuǎn)子系統(tǒng)診斷分析［J］.振動工程學報，2006，19(3):307－212.YAO Hong-liang，LI He，LI Xiao-peng et al.Diagnoses of coupling fault of crack and rub-impact in rotor systems［J］.Journal of Vibration Engineering.2006，19(3):307 －212.

［5］羅躍綱，聞邦椿.雙跨轉(zhuǎn)子系統(tǒng)裂紋－松動耦合故障的非線性響應［J］.航空動力學報，2007，22(6):996－1001.LUO Yue-gang，WEN Bang-chun.Nonlinear responses of two-span rotor-bearing system with coupling faults of crack and pedestal looseness［J］.Journal of Aerospace Power，2007，22(6):996－1001.

［6］黃志偉，周建中，張孝遠，等.水輪發(fā)電機組軸系松動－碰摩耦合故障的動態(tài)響應［J］.西南交通大學學報，2011，46(1):121－126.HUANG Zhi-wei，ZHOU Jian-zhong，ZHANG Xiao-yuan，et al.Dynamic responses of shaft system of hydraulic generator set with coupling faults of bearing looseness and rotor rubimpact［J］.Journal of Southwest Jiaotong University，2011，46(1):121－126.

［7］Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis［J］.Proceedings of the Royal Society London A，1998，454(1971):903－995.

［8］Yang Y，He Y G，Chen J S，et al.A gear fault diagnosis using Hilbert spectrum based on MODWPT and a comparison with EMD approach［J］.Measurement，2009，42:542－511.

［9］唐宏賓，吳運新，滑廣軍，等.基于EMD包絡譜分析的液壓泵故障診斷方法［J］.振動與沖擊，2012，31(9):44－48.TANG Hong-bin，WU Yun-xin，HUA Guang-jun，et al.Fault diagnosis of pump using EMD and envelope spectrum analysis［J］.Journal of Vibration and Shock，2012，31(9):44－48.

［10］邵忍平，曹精明，李永龍.基于EMD小波閾值去噪和時頻分析的齒輪故障模式識別與診斷［J］，振動與沖擊，2012，31(8):96 －101，106.SHAO Ren-ping，CAO Jing-ming，LI Yong-long.Gear fault pattern identification and diagnosis using time-requency analysis and wavelet threshold de-noising based on EMD［J］.Journal of Vibration and Shock，2012，31(8):96 －101，106.

［11］Mandelbrot B B.The Fractal Geometry of Nature［M］.San Francisco:Freeman，1982.

［12］陳 奇，趙 韓，黃 康.齒輪結合面切向接觸剛度分形計算模型研究［J］.農(nóng)業(yè)機械學報，2011，42(2):203－206.CHEN Qi，ZHAO Han，HUANG Kang.Fractal model of tangential contact stiffness between Gears’s Joint Surfaces［J］.Transactions of the Chinese Society for Agricultural Machinery，2011，42(2):203－206.

［13］鐘明壽，龍 源，謝全民，等.基于分形盒維數(shù)和多重分形的爆破地震波信號分析［J］.振動與沖擊，2010，29(1):7－11.ZHONG Ming-shou，LONG Yuan，XIE Quan-min，et al.Signal analysis for blasting seismic wave based on fractal boxdimension and fractal［J］.Journal of Vibration and Shock，2010，29(1):7－11.

［14］郝 研，王太勇，萬 劍，等.分形盒維數(shù)抗噪研究及其在故障診斷中的應用［J］.儀器儀表學報，2011，32(3):540－545.HAO Yan，WANG Tai-yong，WAN Jian，et al.Research on fractalbox dimension anti-noise performance and its application in faultdiagnosis［J］. Chinese Journalof Scientific Instrument，2011，32(3):540 －545.

［15］趙春暉，馬 爽，楊偉超.基于分形盒維數(shù)的頻譜感知技術研究［J］.電子與信息學報，2011，33(2):475－478.ZHAO Chun-hui，MA Shuang，YANG Wei-chao.Spectrum sensing in cognitive radios based on fractal box dimension［J］.Journal of Electronics＆ Information Technology，2011，33(2):475－478.