基于導航定位原理的火箭渦輪泵軸承故障診斷

戴屹梅，張和生，方 柯

(北京交通大學電氣工程學院，北京 100044)

0 引 言

滾動軸承是火箭渦輪泵的關鍵部件，與普通民用軸承相比，火箭發動機渦輪泵的軸承采用固體潤滑，當軸承瞬間受到外力時易導致結構破壞，或瞬間干摩擦引起燒蝕。其故障模式體現在振動信號上的特點是：1) 振動數據中同時含有低頻成分和幾萬赫茲量級的高頻成分；2) 軸承故障特征頻率振動量級小，容易被復雜的背景信號淹沒；3) 軸承表面一旦出現故障，其劣化速度很快；4)振動信號中故障特征頻率具有很強的突發性。這些特點使傳統的振動信號分析方法不能完全適用于火箭發動機渦輪泵軸承故障特征提取。因此，如何準確地實現液體火箭發動機低溫高速軸承的故障特征提取成為相關領域亟待解決的問題。

軸承DN(軸承內徑與軸轉速的乘積)值是軸承技術水平的主要指標，DN值越大，研制難度越高，從設計技術、材料技術、加工技術、試驗技術等各方面都要面臨很高的要求。高DN值軸承技術水平的限制，將會極大地制約發動機性能的進一步提升，高DN值軸承的研制需要充分的試驗驗證能力，因此設計具有超高速采樣頻率的高速軸承試驗臺是實現軸承研制目標的重要基礎，同時需要提高軸承的早期故障診斷能力。

高能振動可能會對其他部件造成損害,并由此產生新的振動頻率信號。這意味著多個振動頻率之間可能存在因果關系。因此,也有必要判斷不同振動頻率發生的時間順序,以找出它們之間可能的因果關系,用于振動信號的故障診斷。

提高這些能力需要探討新的理論途徑,進行跨專業的技術創新，為新型的高速低溫軸承試驗臺設計尋求理論支撐。

常用的渦輪泵軸承故障診斷方法是利用振動信號、聲學信號、溫度信號和介質參數變化量等對故障進行預警[1-3]。振動參數作為液體火箭發動機的重要參數應用于故障診斷工作中，源于航天飛機主發動機(SSME)的研制初期[4]。文獻[5]是NASA/馬歇爾空間飛行中心研究人員提出的一種CPLE分析方法，CPLE技術可以做出一種能夠同時保留振幅和相位的二維頻譜，這在傳統的頻譜分析中是不可能實現的。

安塔瑞斯(Antares)130 首飛失敗，一級火箭的兩臺AJ-26發動機點火15 s 渦輪泵發生故障，E15發動機氧泵轉子與靜子發生碰磨，產生火花點燃液氧導致泵爆炸。專家認為，如果能對飛行前的振動試驗數據發現的異常頻率信號進行振源位置定位，完全有可能避免這次重大失敗。

為了對故障進行診斷，近些年來出現了大量的分析動態非平穩信號的方法，比如短時傅里葉變換[6]，Wigner-Ville分布[7]，小波變換[8]，EMD分解法[9]，隨機共振[10]小波閾值降噪[11]等。形態濾波器(MF)是一種時域算法，能夠直接提取沖擊特性的幾何結構[12]。文獻[13]介紹了一種改進的MF，該方法能夠從低信噪比信號中提取故障特性。文獻[14-18]介紹了幾種先進的時頻分析技術并將其應用到軸承故障診斷上。

由于軸承的過度疲勞使用,軸承部件中往往發生剝落,如保持架、外圈、內圈、滾動體等,因此滾動軸承中遇到的故障類型為保持架故障、外圈故障、內圈故障、滾動體故障[19]。故障診斷的關鍵步驟是從振動信號中提取故障特征頻率，傳統的方法是對實測的頻率分量和理論計算值進行比較。由于滑動和滾動接觸,實測值與理論計算結果不完全相等,有約 5%～10% 的差異。另一方面,由于振動傳感器能夠準確獲得機械設備的工作信號,這些信號不僅包含一些顯示機器健康狀況的信息,而且還含有一些無用的噪音和干擾信號,比如來自齒輪箱和環境等的振動信號。在復雜的機械系統中,這些方法具有一定的內在制約，軸承元件的特征頻率往往與變速箱或其它部件的特征頻率非常接近,這就很難完全依靠特征頻率對故障進行準確的判斷。因此,能夠準確找到振源位置是故障診斷的終極方法。

振源定位法思路清晰，方法簡便，具有通用性。但是定位精度低，要想提高精度，必須提高振動信號頻率，提高時延估算精度，因此出現了各種超聲波定位儀器，包括醫學上使用的超聲波檢查設備。火箭發動軸承的轉速一般為數萬轉，最高可達8萬轉，其故障特征頻率小于6 kHz。這樣的頻率，用傳統的時延估計算法[20]，其精度無法達到故障定位要求。

根據文獻[21] 推導出的時幅曲線表達式：

因此，從理論上講，可以通過不同位置的振動傳感器測得的同一個突發振動信號的相位差來計算出振源的位置。基于這一點，本文提出一種新的突變信號軸承故障診斷方法，將信號處理技術和衛星導航定位原理相結合，發揮二者的優點。此方法類似于定位算法，即將四個振動傳感器和被測軸承設計在一個與衛星導航定位WGS-84類似的直角坐標系中，將系統采樣頻率提高到128 MHz以上，利用時幅曲線分析法精確捕獲故障頻率在各傳感器的出現時刻，將此時間作為輸入條件，計算振動波源位置，進而判斷出軸承故障。

為了表述嚴謹和簡潔，本文提出一個理想模型，在工程應用中，可以利用該文給出的方法，根據精度要求設計不同的實施方案。在實際工程中，多個振動之間往往存在因果關系。因而，利用振動信號進行故障診斷，同樣需要判斷不同振動出現的先后順序，尋找它們之間的因果關系及引起故障的原因。

1 系統構成

任何一次振動均包含四個要素：振源位置、振動強度(幅值)、振動發生的時間、振動頻率。了解這四個要素是徹底了解一次振動的必要條件。從振動數據中提取四個要素中的某一個或幾個，稱之為振動信號特征提取，用提取的特征信號進行故障診斷還需要其他相關知識的配合。目前，對于軸承振動數據的故障特征提取基本集中在特征頻率，振動幅值和振動發生的時間方面，對于確定振源位置幾乎沒有涉及，這是因為火箭發動機渦輪泵軸承的最高轉速為8萬轉/分量級，振動信號頻率遠遠達不到超聲波量級，如果沒有相位信息，幾乎無法達到定位所需求的毫米級精度。

對于氫氧火箭發動機渦輪泵高速軸承試驗臺，在試驗時將軸承置于液氮冷卻箱中，試驗臺的這一結構特點使得在軸承周圍分布多個振動測點成為可能。將試驗系統在結構上按圖1設計，利用振動定位算法與高速采樣同步測控技術以及時幅拐點數值分析方法相結合對軸承故障點進行診斷定位。

圖1 系統坐標系及軸承傳感器結構布局示意圖

以軸承的幾何中心為坐標原點，傳動軸的軸線為X坐標軸，垂直向上為Y軸，垂直向外為Z軸，建立立體直角坐標系。振動傳感器1，2，3，4分布在以R為半徑、以軸承的幾何中心為圓心的球面上。振動傳感器1，2，3，4的坐標分別為1(R,0,0)，2(0,-R,0)，3(0,R,0)，4(0,0,R)。

2 空間坐標定位原理

2.1 空間坐標定位數學模型

設4個振動傳感器測得同一振動的時刻分別為t1，t2，t3，t4，振動波在冷卻液中傳播速度為V，振源的坐標為(X,Y,Z)，4個傳感器的坐標為(Xi,Yi,Zi)，i=1,2,3,4，則各傳感器與振源之間的距離為

設各個傳感器測得同一振動信號的時刻為ti(i=1,2,3,4,)，t0為距振源最近的傳感器測得振動信號的時刻，則t0=min{t1,t2,t3,t4}；各個傳感器測得同一振動信號的時刻與t0的時間差為Ti0=ti-t0(i=1,2,3,4,)。設T為振動從振源傳到距振源最近的振動傳感器所用的時間，那么

dsi=(T+Ti0)V，i=1,2,3,4，即

(1)

式中：T，X，Y，Z是需要求解的未知數。

整理原方程組(1),并將常數記為kij，則：

(2)

將前三個方程聯立：

(3)

(4)

(5)

(6)

(7)

(8)

將式(8)代入式(1)整理得：

T2+k80T+k81=0

(9)

解出T，X，Y，Z，坐標(X,Y,Z)即為振動源所在位置坐標。

2.2 判定故障部位數學模型

根據振源坐標(X,Y,Z)判定故障部位方法。

圖2是軸承在坐標系中的位置圖和其尺寸說明。設軸承內圈內徑為2r1，內圈外徑為2r2，外圈內徑為2r3，外圈外徑為2r4，軸承座外徑為2r5，軸承寬度為h。

圖2 軸承各部件在坐標系中占據的空間位置圖

不難看出，通過坐標(X,Y,Z)的值可以判定故障所在部位。

1)軸承座所占據的坐標空間同時滿足以下數學表達式：

(10)

2)軸承內圈所占據的坐標空間同時滿足以下數學表達式：

(11)

3)軸承外圈所占據的坐標空間同時滿足以下數學表達式：

(12)

4)軸承滾動體和保持架所占據的坐標空間同時滿足以下數學表達式：

(13)

式(1) 和式(10)～(13) 可用于推導故障識別算法,算法流程如圖3所示。

在精度較低的系統中,如果計算出的故障位置坐標位于軸承兩個部件的交界面區域，則可將計算出的振源位置坐標和軸承元件的理論特征頻率進行對比，對故障部位做出準確的判斷。因為軸承不同部件的特征頻率數值差別較大。

圖3 計算振動源坐標和判斷故障部位的流程圖

2.3 模數轉換精度和多通道同步問題

在工程上實現本文提出的方法，需要高速多通道同步采樣板，采樣頻率大于100 MHz，通道數大于4，各通道的同步誤差小于0.1 ns。以下是SPECTRUM INSTRUMENTATION公司的一款產品，M4i.44xx-x8-14/16性能完全滿足要求。

四通道采樣頻率為500 MS/，所有通道同步采樣，通道之間同步采樣時間誤差小于60 ps。

2.4 軸承故障信號時域波形

軸承故障信號時域波形一般較為復雜，為振蕩衰減的周期信號，文獻[21]做了詳細的介紹。振蕩衰減周期信號經過傅里葉變換，其能量大部分集中在基頻正弦信號上。所以在頻域提取的故障信號都是正弦信號，故障仿真信號用3589 Hz正弦信號。

2.5 坐標精度問題

本文給出的是理論直角坐標系的計算公式和故障定位判定方法，為了簡化方程,便于理解,本文選取了4個特殊坐標。在實際工程實現中，四個振動傳感器和軸承所占據的空間位置與理論位置存在誤差。用三坐標機測量四個傳感器和軸承幾何中心的坐標,并根據這些測量結果再次建立較為復雜的定位方程如式(1)和四故障判斷不等式如式(10)～(13)。目前,市場上三坐標機的精度可達(2.5±L/300) μm,本文取L=R=300, 所以坐標測量精度達到 3.5 μm,完全滿足精度要求。

2.6 噪聲的影響

真實的試驗數據中包含各種噪聲，為了考核噪聲環境下該方法的可行性，本文采用了工程仿真的方法，構造仿真數據時域信號時，將故障信號疊加到真實的試驗數據中，可最大限度模擬真實環境下的故障信號，然后再用時幅曲線拐點分析法對故障信號發生的時刻進行提取。

3 仿真校驗

3.1 構造傳感器仿真時域信號

(14)

設ts為振源發出振動信號的時刻，各個傳感器測得信號的時刻為：

(15)

各傳感器捕獲到同一信號的時間差分別為

由此得出仿真程序的輸入輸出參數。

輸入參數：(X′，Y′，Z′)，R，V，K

將特定頻率的振動信號依據時間差加到真實的試驗數據中，即可獲得仿真計算所需要的時域信號試驗數據。

設軸承結構參數為r1=17.5 mm，r2=21.5 mm，r3=27 mm，r4=30 mm，r5=34 mm，h=14 mm。振源坐標為(5,0,28)，位于軸承外圈區域，R=300 mm，振動波在液氮中的傳播速度為V=1168000 mm/s，采樣頻率為k=128000000 Hz。

將以上參數輸入到求解仿真參數程序，計算結果:

2.580011×10-4, 2.580011×10-4,

2.329161×10-4}=2.329161×10-4

(2.329161×10-4)=0.207875×10-4

(2.329161×10-4)=0.250850×10-4

(2.329161×10-4)=0.250850×10-4

(2.329161×10-4)=0.0×10-4

3.2 合成振動傳感器故障時域信號

圖4 真實試驗數據時域圖

圖5 傳感器1合成后試驗數據時域信號

用同樣的方法合成其余三個傳感器的時域數據。

3.3 仿真計算和故障診斷

利用第3.2節合成的傳感器仿真時域信號計算振源位置并進行故障診斷。

3.3.1時幅曲線拐點分析算法簡介

根據參考文獻[21]所介紹的時幅曲線拐點分析算法,以Δt為時間步長，信號采樣頻率Fs。

在(T1+nΔt，T2)時間段內逐次取信號序列fj的傅里葉變換中頻率k的幅值Akn，則有

(16)

其中L為分析點數，L=N-nΔtFs；N為原序列的總采樣點數，N=(T2-T1)Fs，n≤N。

只要采樣頻率Fs足夠高,計算步長Δt足夠小,就可以得到所需精度的信號出現和消失時刻。

為了直觀理解該算法和時幅曲線，構造函數f0(t)

(17)

對函數f0(t)進行12800 Hz采樣，取分析時間從0～15 s，其時幅曲線如圖6所示。

圖6 f0(t)在200 Hz的時幅曲線(0～15 s)

分析時間從0～9 s，此時時幅曲線如圖7所示。

圖7 f0(t)在200 Hz的時幅曲線(0～9 s)

從圖6和圖7可以看出,時幅分析法能夠準確計算出信號的出現時刻和消失時刻。不同的分析時段對應不同的曲線形狀，但信號拐點出現的時刻是不變的。

3.3.2用時幅曲線分析算法求出各傳感器合成數據中3589 Hz信號出現的時刻

為了節省計算時間,先用Δt=0.001 s計算0～0.03 s時段范圍內傳感器1的3589 Hz時幅曲線，得到3589 Hz信號出現時刻大致在0.01 s，如圖8所示。再縮小分析時段范圍，在0.00975～0.01015 s之間，用Δt=1×10-8s計算3589 Hz的時幅曲線如圖9～圖12所示，選取幅值達到峰值時的時刻為ti(i=1,2,3,4)。

圖8 0～0.03 s 3589 Hz時幅曲線

1)傳感器1數據分析

從圖9所示傳感器1的時幅曲線可以看出，3589 Hz信號出現的時間為t1=1.000567×10-2。

圖9 傳感器1數據在3589 Hz的時幅曲線

2)傳感器2數據分析

從圖10所示傳感器2的時幅曲線可以看出，3589 Hz信號出現的時間為t2=1.000825×10-2。

圖10 傳感器2數據在3589 Hz的時幅曲線

3)傳感器3數據分析

從圖11所示傳感器3的時幅曲線可以看出，3589 Hz信號出現的時間為t3=1.000979×10-2。

圖11 傳感器3數據在3589 Hz的時幅曲線

4)傳感器4數據分析

從圖12所示傳感器4的時幅曲線可以看出，3589 Hz信號出現的時間為t4=9.984540×10-3。

圖12 傳感器4數據在3589 Hz的時幅曲線

t0= min{t1，t2，t3，t4}=min{1.000567×10-2，

1.000825×10-2，1.000979×10-2，

9.984540×10-3}=9.984540×10-3

T10=t1-t0=(1.000567×10-2)-

(9.984540×10-3)=0.2113×10-4

T20=t2-t0=(1.000825×10-2)-

(9.984540×10-3)=0.2371×10-4

T30=t3-t0=(1.000825×10-2)-

(9.984540×10-3)=0.2371×10-3

T40=t2-t0=(9.984540×10-3)-

(9.984540×10-3)=0.0×10-4

3.4 利用振動定位算法對軸承故障進行定位診斷

將以上時幅曲線分析法找出的各個傳感器測得3589 Hz信號出現的時刻，代入公式計算出振源坐標，判斷是否是軸承故障及故障部位。圖3為定位故障算法流程圖。

將時間t1,t2,t3,t4輸入定位算法故障診斷程序，經定位算法故障診斷程序計算,得到振源位置坐標是X=3.905743，Y=-0.902871，Z=27.395615，與軸承結構參數r1，r2，r3，r4，r5比對后判定，該振源位置位于軸承外圈區域。第3.1.2節用于構造模擬數據的振源坐標為(5，0，28),與仿真計算出的振源坐標(3.905743，-0.902871，27.395615)的最大誤差為 1.1 mm (5-3.905743=1.094257)。

4 結 論

由于時幅曲線分析法能夠依據采樣頻率精度準確找出信號發生突變的時刻，這為振源空間定位分析提供了可能。其理論定位精度只與傳感器的坐標精度、采樣頻率和多路采樣器的同步精度有關；與傳感器位置、振源的位置均無關。

通過以上仿真計算，在采樣頻率為128 MHz，傳感器的坐標誤差為零、多路采樣器的同步誤差為零、波的傳播速度為1168 m/s的條件下，最大誤差為±1.1 mm。

在實際應用中，只要判定振源位置在軸承結構包絡范圍內，即可將故障定位于軸承；通過增加定位精度可進一步判定軸承故障的具體部位，此方法亦可與故障特征頻率判別法進行相互驗證。

每一個振動信號都包含四個要素，即振源位置、振動發生的時間、振動的幅值(能量)、振動的頻率。掌握了這四個要素也就徹底了解了與這四個要素相關的振動。利用現代技術手段從振動信號中發掘出這四個要素是通過振動信號進行故障診斷的最根本的途徑。本文利用信號時幅曲線所包含的相位信息計算振源位置就是在這一方面做的一次嘗試。仿真計算結果表明該方法理論上可行，有實際工程價值。