陳 偉, 顧曉華, 張永強
(中國航發商用航空發動機有限責任公司,上海 200241)
振動模態參數對于結構動態特性優化有著重要且不可替代的作用。自20世紀70年代發展至今,試驗模態分析成為實驗室條件下靜態獲取結構模態參數的可靠手段。通過力錘或激振器實現人工激勵,同時測量結構的響應信息,基于系統的輸入和輸出構造頻響函數,最后使用多樣的參數辨識算法進行模態參數的辨識。然而,在役結構工作環境下的邊界條件有別于實驗室環境,載荷條件也更加復雜。
為了解決在役結構模態參數辨識問題,工況模態分析(或稱工作模態分析、環境激勵下模態參數辨識)在20世紀90年代應運而生[1-2],僅基于被測結構工況環境下的響應便可辨識其模態參數。工況下結構的模態參數辨識可真實反映結構在真實邊界條件下的動力學特性,且認為在結構響應中占優的頻域成分即為對結構影響較大的模態。對處于復雜環境或載荷條件下的試驗件,如轉子葉片的“硬化”特性、流固耦合、熱變形等,工況模態參數辨識算法依然適用,且能夠捕捉結構在復雜環境下的真實動力學特性。此外,該測試方法依賴于環境激勵或結構本身作為載荷,極大地減少了人工激勵的成本,避免了對試驗件的損傷[3]。
大型機械裝備往往含有轉子部件,如航空發動機、燃氣輪機、直升機、汽車等等。轉子部件在運行過程中必然帶來周期性諧波激勵[4]。工況模態分析的理論框架建立在寬頻隨機激勵(或白噪聲)的基礎上,因此認為測量數據所識別的模態皆與被測系統的極點有關。然而,由旋轉部件帶來的諧波導致載荷“非白”時,依賴原有方式進行模態分析就違背了白噪聲假設,使得模態參數辨識變得更加困難,導致的困難可總結為3點:首先需要從辨識的結果中區分出結構模態和“諧波模態”;其次,若輸入諧波的頻率接近甚至與系統的固有頻率重疊時,辨識方法將無法準確識別模態參數;最后,由于強諧波載荷的存在,結構測量響應中諧波響應往往占據主導地位,導致模態參數辨識需要更高階次的模型進行曲線擬合,辨識算法可能失去魯棒性。
針對諧波激勵下的模態參數辨識問題,學者們進行了廣泛的研究,如于亮亮等[5]針對激勵中含有諧波載荷時的工況模態分析,論證了人為添加諧波極點的理論可行性,并對特征系統實現(eigensystem realization algorithm,ERA)模態參數辨識算法進行了修改和驗證。張義民等[6]基于隨機響應與諧波信號的概率密度函數不同,區分了諧波頻率,并將該方法用于國產轎車的振動模態參數辨識中。夏遵平等[7]基于諧波信號和結構模態譜峭度的不同,提出了基于譜峭度的諧波識別方法,并應用于工況下模態參數辨識。Pintelon等[8]提出了非參數模型的諧波去除方法,并將該方法應用于直升機測試數據中。類似在工程結構中的應用已有許多,如風機、汽車等,這類方法大致可歸類為諧波識別類算法、諧波去除算法[9],但鮮有算法可同步實現諧波的識別與去除。
本文所提出的基于倒頻譜的工況模態參數辨識方法,旨在識別含轉子結構在運行狀態下所受諧波激勵造成的響應,以及去除諧波載荷對結構響應的影響,優化工況下振動模態參數辨識。本文首先對倒頻譜的定義進行概述,并敘述了倒頻譜編輯法的流程。然后,通過三自由度系統算例驗證方法的可行性。接著,引入實驗室條件下自由邊界實心梁的試驗進一步驗證方法的可靠性。最后,應用于國產某型航空發動機試車數據,驗證該方法的適用性。
倒頻譜在齒輪和軸承的故障診斷中有著廣泛的應用。倒頻譜早在1963年便被引入,定義為“對數功率譜的功率譜”[10],被用于在信號中檢測回聲信號,該種定義的缺陷在于倒頻譜不可逆。在快速傅里葉變換發展之后,倒頻譜更為廣泛的定義為幅值譜取對數運算后的逆傅里葉變換[11]。本文中倒頻譜定義為傅里葉譜取對數運算后的逆傅里葉變換,即
Cc(τ)=F-1{In[X(f)]}=F-1{In[A(f)]+jφ(f)}(1)
式中:τ對應倒頻譜的時間;X(f)為時域信號x(t)的頻域形式,即
X(f)=F{x(t)}=A(f)ejφ(f)
(2)
式中:A(f)為頻域幅值函數;φ(f)為頻域相位函數。
式(1)所定義的倒頻譜為復數形式,有幅值和相位信息。實倒頻譜為幅值譜取對數運算后的逆傅里葉變換,定義為
Cr(τ)=F-1{In[A(f)]}
(3)
周期譜成分,如諧波和其倍頻,在倒頻譜中的表現為間距相同的峰值,因此可通過陷波器將其在倒頻譜中去除。通過對實倒頻譜的編輯可實現對信號幅值的調節,分離出確定性周期信號和隨機信號,再結合原有頻譜的相位實現不改變原信號的相位信息。
使用陷波器對頻譜進行處理時,面對多組諧波易出現處理不徹底,或誤編輯的情況,從而更改了原有模態成分。將時域信號變換到倒頻譜域,會使結構模態信息在信號頭部集中。隨著式(3)中τ的增大而使模態信號成分降低,因此本文使用前端矩形窗結合后端指數窗的形式對倒頻譜進行編輯處理。修改后的倒頻譜再變換到時域即實現了對諧波成分的過濾,極大地保留了信號中的模態響應成分,諧波響應對應的頻譜峰值會極大降低。再通過對比處理前后的響應譜,從而實現諧波峰值的識別。因此,倒頻譜編輯法可同時實現諧波響應的識別和去除,其具體步驟如下:
步驟1將時域信號進行快速傅里葉變換,得到響應的傅里葉譜。
步驟2對傅里葉頻譜的幅值和相位進行分離,即分離出幅值譜和相位譜。
步驟3對幅值譜取對數,并進行逆傅里葉變換,得到實倒頻譜。
步驟4對實倒頻譜進行編輯,即加窗(矩形-指數窗),實現倒頻譜中諧波成分的去除。
步驟5對編輯后的倒頻譜進行快速傅里葉變換,得到對數幅值譜。
步驟6將處理后的幅值譜進行指數運算,并結合相位譜得到去除諧波的傅里葉譜。
步驟7對傅里葉譜進行逆傅里葉變換,獲得去除諧波的時域信號,該信號可用于常規工況模態參數辨識。
具體的流程圖如圖1所示

圖1 倒頻譜編輯方法流程圖Fig.1 Flow chart of the cepstrum editing method
關于實倒頻譜的編輯,主要采用矩形窗結合指數窗進行。矩形窗保留了原信號中的模態成分,指數窗則消減了信號中的諧波成分,通過兩者結合實現諧波的濾除。窗函數的表達式如下
(4)
式中:m為正整數,表示矩形窗的長度,因實倒頻譜的不同而具體確定,文中m對應響應包絡衰減90%時的信號點數;λ為正實數,表示衰減系數。
衰減系數可表達為
(5)
式中:ξ為阻尼比系數,根據諧波的強弱進行調節;fs為采樣頻率;fh為諧波基頻的頻率值。
為了驗證所提出倒頻譜方法的有效性,建立一個三自由度系統,其固有頻率和阻尼比如表1所示,其質量矩陣M、剛度矩陣K和阻尼矩陣C分別為

表1 構造系統的前3階固有頻率及阻尼比
在3個質量塊作用了標準差為10 N滿足正態分布的寬頻隨機噪聲,在第一個自由度還作用了兩組諧波:第一組基頻為7.1 Hz的1倍~8倍頻諧波激勵,第二組的基頻為16 Hz的1倍~4倍頻諧波,兩組諧波的幅值為1~50 N的隨機整數。數據記錄的采樣頻率為1 024 Hz,采樣時間為300 s。作為例子,圖2展示了第一個自由度m1測點響應的傅里葉譜。

圖2 系統在諧波和隨機載荷共同作用下響應譜Fig.2 Response spectrum of the system under the combination of the harmonic excitations and the random excitations

將倒頻譜編輯法應用于算例中的響應信號。首先,將時域響應通過快速傅里葉變換轉換到頻域中,再將幅值譜取對數,并通過逆傅里葉變換得到實倒頻譜,變換過程參考流程圖(圖 1),得到的實倒頻譜見圖3(a)。將長度為500的矩形窗結合指數窗的窗函數(見圖3(b))施加于實倒頻譜,獲得修改實倒頻譜,見圖3(c)。最后,對修改實倒頻譜進行快速傅里葉變換并進行指數運算,獲得修改幅值譜,修改前后幅值譜的對比見圖 3(d)。

圖3 倒頻譜編輯法的處理過程Fig.3 Results using cepstrum editing method
由圖 3(d)可知處理后的頻譜中不含有諧波成分,僅含有模態峰值。因此在未知頻率諧波激勵下,通過該方法可識別并去除頻譜中由諧波激勵引起的峰值,進而確定結構固有模態導致的峰值。此外,通過諧波補償法處理后的信號更加光滑,更易于識別頻譜中的結構模態參數。
下面對處理后的響應進行工況模態辨識,使用了相關函數法[12-13](或稱自然激勵技術),即基于相關函數識別結構的模態參數。對去除諧波響應后的響應進行工況模態參數辨識,首先提取了響應的相關函數,通過ERA辨識算法處理相關函數,提取系統的極點,辨識的穩態圖如圖4所示。穩態圖極點清晰,系統的前3階模態處聚集著穩定的極點。識別的固有頻率和阻尼比結果如表2所示。為了說明諧波處理算法的有效性,表 2同時展示了系統真實固有頻率和阻尼比,以及直接辨識的結果,即直接對相應的相關函數使用ERA辨識算法的結果。

表2 識別結果的對比Tab.2 Comparison of the identified results

圖4 模態參數辨識穩態圖Fig.4 Stabilization diagram of the modal parameter identification
由表2結果可知,在不去除諧波影響下辨識,第一階模態的阻尼比為0.01%,遠小于結構的模態阻尼比,易造成較大誤差。對于諧波與結構固有頻率重疊的模態,其阻尼比識別誤差較大,當使用所提方法去除諧波影響后提取的阻尼比接近真實值。
下面進行概率密度函數法(probability density function, PDF)、譜峭度法(spectrum kurtosis,SK)、修改ERA(modified ERA, MERA)法、非參數模型法[8](non-parameter,NP)與本文所提方法倒譜編輯算法(cepstrum editing,CE)在識別諧波、去除諧波和處理諧波頻率靠近模態頻率情況下的對比分析。
PDF方法通過選擇頻譜中的峰值,以峰值頻率為中心進行窄帶濾波,并對濾波后的響應以幅值為橫坐標以數量為縱坐標繪制直方圖代替概率密度曲線。圖 5(a)展示了使用PDF方法的處理結果,由于第一階模態和第三階模態附近存在諧波峰,PDF處理結果展示為諧波的特性,易造成結構模態的誤判。
譜峭度是基于響應譜的平均計算獲得,圖 5(b)的虛線展示了譜峭度,圖5中每個諧波峰值處對應的譜峭度值接近-1,可實現諧波的檢測。但諧波靠近固有模態頻率時,易發生將結構模態判斷為諧波響應的情況。

圖5 不同方法處理結果圖Fig.5 The results of different methods

圖6 測試裝置布置圖Fig.6 The test setup
圖 5(c)展示了基于MERA方法的穩態圖,相較于本文所提出處理算法的穩態圖,MERA方法需要更高階次的擬合模型獲取穩定的極點,雖然諧波的模態被單獨考慮,但由于諧波頻率接近第一和第三階模態且諧波幅值顯著高于模態響應幅值,結構模態依然很難識別。此外,MERA方法需要在已知諧波頻率的情況下進行,在諧波數量多且頻率變化的情況下很難適用。
圖 5(d)展示了非參數模型法的處理結果,該方法直接對響應譜進行編輯,手動去除諧波峰并進行插值。該方法在處理諧波頻率靠近結構模態頻率時易產生較大誤差,如圖 5(d)中的放大圖所示,第一階模態峰值直接被手動去除并使用線性插值補全,模態峰值處平坦,人工修改痕跡明顯,相較之下本文所提的倒譜編輯法在去除諧波影響情況下保留了結構模態響應,如放大圖中虛線所示。
表3總結了上述PDF、SK、MERA、NP方法和本文所提出的CE方法在檢驗諧波頻率、去除諧波影響和處理諧波頻率靠近模態頻率情況的效果。根據表3,本文所提方法是對比方法中同時可以實現檢測諧波和去除諧波影響的方法,且具備處理諧波頻率靠近模態頻率的能力,識別精度優于MERA方法。

表3 不同方法的對比Tab.3 Comparison of different methods
試驗件為實心鋼梁,其尺寸為720.0 mm×25.4 mm×25.4 mm。如圖 6所示,梁的兩端通過彈性繩懸掛,模擬自由邊界條件。試驗通過激振器模擬轉子激勵,激振器通過頂桿和力傳感器連接在自由梁的3號位置,7個100 mV/g的微型加速度傳感器均布在梁的長度方向,并通過數采系統記錄梁的響應信號。
信號采樣頻率為6 400 Hz,每次測試記錄120 s。在引入諧波載荷之前,對試驗件進行傳統寬頻隨機模態測試,獲取梁的前4階模態參數,如表4所示。隨機模態測試結果將作為工況模態辨識結果的參考值。

表4 系統前4階固有頻率和阻尼比Tab.4 The first 4 natural frequencies and damping ratios

表5 穩態圖中對應穩定極點開始階次Tab.5 The start order for stable poles
在控制器端加入235 Hz, 470 Hz,620 Hz的諧波信號及信號有效值為諧波信號1/10的寬頻噪聲信號模擬多轉子系統激勵,設計的信號通過激振器作動于試驗件之上。由于激振器與試驗件的相互作用,結構響應譜中會出現激勵諧波倍數相關的復雜諧波響應,測點1的響應時程和功率譜如圖7所示。

圖7 測量響應信號Fig.7 The measured response from the test
下面將倒頻譜編輯法應用于獲取的響應信號。倒頻譜編輯使用了點數為500的矩形窗結合衰減系數λ為0.2的指數窗。使用倒頻譜編輯算法處理前后信號的響應譜對比見圖 8,相較于原始信號(圖8中實線)處理后的信號(虛線)中諧波成分有了顯著下降,且模態峰值保留較好。

圖8 倒頻譜編輯法使用前后頻譜對比圖Fig.8 The comparison of the power spectrum density functions of using and without using the proposed cepstrum editing method
下面對去除諧波和未去除諧波的響應分別進行工況模態分析,使用了相關函數法進行前處理,并結合ERA辨識算法形成穩態圖進行參數選取。表 5對比了去除諧波和未去除諧波獲得穩定極點的ERA擬合階次,第一階模態進行諧波處理后,擬合階次從39降低至6。數據顯示去除諧波后獲得穩定極點的擬合階次顯著下降。
以傳統模態測試獲取的模態參數為參考,計算諧波處理前后辨識的前4階模態對應的固有頻率和阻尼比與參考值之間的誤差,結果如圖9所示。未進行諧波處理時,識別的第3階模態的頻率誤差為0.41%,為辨識的最大頻率誤差。頻率誤差整體相對較小,尤其體現在諧波處理后,頻率辨識誤差均在0.1%以內。圖 9(b)中,進行諧波處理后辨識阻尼比的誤差有了顯著下降,第1階模態阻尼比的相對誤差從204.0%下降至12.2%,去除諧波后阻尼比的辨識相對誤差均小于25%。

圖9 去除諧波前后辨識固有頻率和阻尼比相對誤差對比Fig.9 The comparison of relative errors of identified natural frequencies and damping ratios from with processing harmonic data and without

圖10 測試系統布置說明Fig.10 The explanation for the experimental setup
因此,通過文中所提的倒頻譜編輯算法可以有效去除諧波影響,使辨識算法的擬合階次更低,辨識模態參數的精度更高。
發動機的核心機進行大氣進氣試驗,測試布局如圖 10所示,壓電式高溫加速度傳感器布置在相應測試截面用于采集振動響應信號,并通過電荷放大器將電荷信號轉換為電壓信號傳入數采調理儀,試驗中設定采樣率為25 kHz,數采調理儀實現模擬信號轉換為數字信號并傳輸至服務器端,服務器進行數據的存儲和計算,根據需要將實時處理的結果通過監視電腦展示。
分析測點位于試驗件后承力機匣順航向順時針9點鐘位置,關注頻帶為0~2 000 Hz,因此對獲取的信號進行了重采樣,處理后采樣率為5 000 Hz,在試驗件額定轉速8 250 r/min,10 150 r/min,11 500 r/min記錄了3 min數據。實測發動機試驗件受到熱載荷、氣動載荷、轉子激勵等復雜激勵,轉速會存在一定波動,以8 250 r/min轉速下測點響應的功率譜為例,具體如圖11所示。由圖11可知,結構所受載荷復雜,除了轉子基頻及分量諧波,還有多樣周期性動載荷激勵。

圖11 響應的功率譜密度函數Fig.11 Power spectrum density function of the recorded response
將所提倒頻譜編輯算法應用于測量的響應信號,處理后的實倒頻譜如圖12(a)所示,對實倒頻譜加窗編輯,其中矩形窗的長度為10,指數窗的衰減系數λ為0.18。將編輯后的倒頻譜通過圖1的流程變換回時域,計算修改后信號的功率譜密度函數,修改前和修改后的功率譜對比如圖12(d)所示。經過倒頻譜編輯后結構模態成分被保留,而周期諧波信號被濾除,可依據處理前后功率譜中峰值的變換判斷對應峰值是否由周期諧波激勵導致。圖12(d)中235 Hz和328 Hz對應的頻譜峰未見顯著變化,因此該處頻譜峰為結構模態導致,進一步可基于去除諧波影響的信號進行模態參數辨識。需說明,工況模態分析僅提取試驗件被激發的主要模態,235 Hz和328 Hz為發動機運行過程所激發起來靜子機匣的主要模態特征,被激發的模態與機匣本身結構和運行環境有關。

圖12 諧波去除過程圖Fig.12 Process of harmonic response removal
對去除諧波后的時域響應進行基于相關函數的工況模態分析,即對響應信號做相關運算后,再使用ERA算法結合穩態圖進行模態參數辨識。辨識穩態圖如圖13所示,辨識的模態參數結果如表6所示。圖 13中235 Hz和328 Hz峰值處系統辨識結果清晰穩定,去除諧波影響后模態參數辨識更為簡單。

表6 核心機在不同轉速下機匣模態參數辨識結果

圖13 模態參數辨識穩態圖Fig.13 Stabilization diagram of the modal parameter identification
通過對比去除諧波成分前后響應譜,可判斷頻譜中對應結構模態的峰值,獲得了清晰的穩態圖,從工況試車數據中識別了結構的振動固有頻率及阻尼比參數。由于發動機結構動力學特性受到其運行狀態影響(溫度、壓力等變化),在不同轉速下其固有頻率和阻尼比也發生不同程度變化,在10 150 r/min和11 500 r/min轉速下識別的模態參數見表6,隨著轉速的升高,結構固有頻率有著較為顯著的下降趨勢。綜上,所提方法可在發動機運行工況下識別其固有頻率和阻尼比,可在實際環境下識別發動機機匣的動力學特性,有利于追蹤其動力學特性隨試驗條件的變化。
針對含轉子結構在運行工況下受諧波載荷影響的模態參數辨識問題,本文提出了基于倒頻譜編輯的工況模態參數辨識方法,可實現對恒定諧波激勵下諧波響應的識別和去除,通過三自由度系統算例、自由梁試驗和發動機試車數據驗證了方法的適用性和可靠性。具體結論如下:
(1)本文所提倒頻譜編輯法可識別諧波頻率和去除諧波的影響,算法操作簡單,且在多組諧波共同作用的情況下,該方法依然適用。
(2)對于諧波激勵頻率與結構模態頻率重疊的情況,該方法依然可以有效去除諧波的影響,提取結構固有模態參數。
(3)所提方法對轉子頻率波動不敏感,當實際含轉子結構轉速波動時依然適用,可識別諧波和去除諧波影響并辨識結構的固有模態。