渦輪機葉片同步振動參數辨識方法研究(建模仿真)

王維民， 任三群， 陳立芳， 黃 山

(北京化工大學 高端機械裝備健康監控與治愈化北京市重點實驗室, 北京 100029)

渦輪機葉片同步振動參數辨識方法研究(建模仿真)

王維民， 任三群， 陳立芳， 黃 山

(北京化工大學 高端機械裝備健康監控與治愈化北京市重點實驗室, 北京 100029)

基于葉尖定時非接觸式葉片振動測量技術，提出一種需要較少傳感器又能提高測量精度的葉片同步振動參數分析方法。通過低轉速到高轉速運行過程，測量出葉片振動位移隨轉頻變化之間的關系，根據非線性最小二乘擬合算法獲得葉片諧共振點處的中心頻率、振動最大幅值、相位、振動恒偏量，但無法獲得葉片振動倍頻值。對諧共振動點附近處的頻率值進行恒速運轉，不同位置處的傳感器獲得不同的振動位移，通過GARIV算法獲得各葉片在對應諧共振頻率點處的倍頻值，結合掃頻擬合獲得的葉片振動參數可做出葉片坎貝爾圖。基于以上的識別方法，提出了葉尖計時傳感器的布置方法，并對實際測量中可能存在的干擾因素進行了詳細的仿真分析，驗證了方法的可行性。該辨識方法的實驗驗證將在《渦輪機葉片同步振動參數辨識方法研究(實驗研究)》一文中給出。

渦輪葉片； 同步振動； 參數辨識； 葉尖定時(BTT)； GARIV

葉尖計時測量法常被用于非接觸式應力測量系統中以評估葉片的工作應力[1]。也可以被用于監測葉片和輪盤的工作狀態，葉片振型預測以及阻尼比識別等[2]。該方法主要通過安裝在機匣上的葉頂傳感器監測葉片實際到達時間與理論到達時間的差值來獲得葉片的振動情況。在研究渦輪葉片振動特性的過程中，主要分為異步振動和同步振動兩種。同步振動較于異步振動而言，由于單支葉頂傳感器在監測過程中監測到的葉片振動值不變，這對葉片振動參數辨識帶來較大的難度，所以需要布置較多的傳感器才能辨識出葉片的同步振動參數。但在實際的工程測量中，可安裝的測點數量被嚴格的限制，進而言之使用較少的傳感器來辨識葉片振動特性是非常有必要的。

葉片同步振動參數識別方法有雙參數法[3-4]，但對噪聲較為敏感，抗干擾較差。正弦曲線擬合法[5-6]，要求布置較多的傳感器。單獨使用GARIV方法時[7]，可獲得葉片諧共振中心的頻率值，但葉片振動幅值等參數的求取誤差較大。文獻[8]在正弦曲線擬合法的基礎上提出了多傳感器掃頻擬合方法，雖然可以通過倍頻遍歷算法獲取葉片諧共振中心處的倍頻值，但需安裝傳感器的數量較多。

本文通過比較GARIV法和掃頻擬合法的優點，提出了新的同步振動參數辨識方法。該方法具體為先通過掃頻擬合確定同步振動中諧共振點處的中心頻率、幅值等參數；然后對存在的諧共振中心進行恒速運轉，通過GARIV方法獲得對應諧共振中心的倍頻值，進而獲得較為完整的葉片同步振動參數。在測試方案確定上，先將葉片對象模型化，通過對模型進行建模仿真來確定合適的傳感器布置，安裝角度的允許誤差，測量時合適的運轉轉速值，以及多倍頻存在時的辨識可行性等都做了驗證，并針對模型給出了具體的要求。對于掃頻擬合方法，通過仿真給出了合適的升速值，以及在數據處理過程中使用的擬合方法，以期提高擬合的準確性。最后通過一個仿真實例來綜合展示新方法的辨識準確性。

1 葉片振動倍頻值求解算法

1.1自回歸算法(AR法)求解振動倍頻值

Heath給出了自回歸方法對同步振動下求解倍頻值的具體推導過程。其中，假設葉片為單自由度無阻尼振動模型，且運轉在諧共振中心頻率處，通過推導可以獲得使用4支等間距分布的傳感器時，單葉片振動參數求解方程組為

(1)

式中:x1，x2，x3，x4分別為四支傳感器監測到的葉片振動位移。a1=2cos(ωnΔtp);ωn表示葉片諧共振中心處振動角頻率;Δtp表示葉片經過等間距傳感器的時間間隔;DB為葉片振動的恒偏量。通過式(1)可以求解出a1和DB的值。

葉片同步振動過程中諧共振中心對應的倍頻值可表示為

(2)

式中:n表示倍頻值;Ω表示轉頻值。將a1=2cos(ωnΔtp)代入式(2)中可得：

(3)

當傳感器等間距分布之間的夾角為θ時(與下文中的傳感器安裝夾角相同)，則Δtp表示為

(4)

代入式(3)得：

(5)

推導中假設葉片組運轉轉速為諧共振中心處轉速，但通過分析式(5)可知運轉轉速對振動倍頻的求解沒有影響。進而言之，只要a1值的求解準確，則在任何運轉轉速下都可以準確求解出振動倍頻值。但考慮到實際情況中的誤差、干擾及多倍頻等的存在，需要對可運行轉速做進行進一步的分析。

1.2GAR法及GARIV法

1.2.1 GAR方法

AR法是基于單葉片單圈數據進行分析的。考慮整機分析和測量分析的實時性，在AR的基礎上進行擴展，提出GAR方法。通過對式(1)擴展得到GAR法的計算式為

(6)

式(6)中xR,B,p表示葉片B在第R圈時經過傳感器p采集到的振動位移。上述公式僅僅在維度上進行了擴展，沒有在實際測量過程中可能存在的干擾噪聲等進行有效的處理。在GAR的基礎上為了減弱噪聲等干擾的影響，引入了輔助變量，即GARIV法。

1.2.2 GARIV方法

(7)

當測量過程中有誤差時，記誤差為εi，則有：

(8)

代入式(8)有：

(9)

上式可簡寫為如下形式：

Y=Xθ

(10)

引入輔助矩陣Z，該矩陣通過對矩陣X進行適當的平移后獲得的，具體如下：

(11)

其中平移k值的選取為k=NAReq×R，式中NAReq為單圈單葉片可成立等式個數，如4個傳感器時可成立的等式個數為2；5個傳感器時，可成立的等式個數為3；R值的選取不定，k值往往選取可移動總數的10%。在實測數據處理中發現，平移k值對倍頻值求取的影響較小。

對式(10)兩邊同時左乘ZT，則有：

ZTY=ZTXθ

(12)

由干擾的不相關性和假設測量誤差的均值為0，可通過式(12)將包含誤差εi的項消除，則測量誤差對采用GARIV方法的影響較小，進而起到提高測量精度的要求，提高抗干擾能力。詳細證明過程可參考文獻[7]。

2 葉片振動倍頻值測量方案

為保證辨識方法在實際測量使用時的可行性，需通過建模仿真對實際使用中可能出現的各種情況進行研究，最后給出有效可行的測量方案。

2.1仿真模型簡介

為模擬實際情況下葉片振動，對實驗臺中使用的葉片模型化，見圖1。建立質量-剛度-阻尼模型，不考慮葉片之間存在耦合，進而可對單個葉片做同步振動特性分析。

圖1 從實際煙機葉片到實驗臺使用的簡化葉片

2.2相關影響因素對AR法求取葉片振動倍頻值的影響

2.2.1 單倍頻時傳感器安裝夾角對倍頻值求解的影響

仿真中默認激振力僅激起單倍頻下的葉片同步振動，傳感器安裝夾角變化范圍為1°～120°，步長1°。仿真結果如圖2所示，圖(a)～(d)分別為1倍頻，10倍頻，20倍頻和30倍頻時，通過AR法求解的倍頻值與傳感器安裝夾角之間的關系。

通過圖2分析可知，隨著倍頻值的增大，為通過AR法獲得準確的倍頻值，則傳感器允許的安裝夾角在減小。例如分析1倍頻時傳感器允許安裝夾角不超過20°，分析30倍頻時傳感器安裝夾角不應超過6°。

(a)EQ=1(b)EQ=10

(c)EQ=20(d)EQ=30

圖2 不同倍頻下傳感器安裝夾角對倍頻值求解的影響

Fig.2 The influence of the angle between sensors on solving the EO under different EO

通過仿真計算發現傳感器1與鍵相傳感器的初始夾角不會對AR法分析的結果產生影響。當使用AR法(GAR、GARIV法)時，傳感器布置方面的要求為，根據葉片可能激起的倍頻值確定合適的傳感器布置夾角，這樣才能通過AR法較為準確的獲取葉片同步振動下的倍頻值。

2.2.2 多倍頻存在時AR法求解倍頻值的準確性

在實際測量時，更多的是激起多倍頻的情況，此時就需要研究AR法(GAR、GARIV法)能否很好的解決多倍頻存在時的倍頻值求取問題。

假設葉片模型中存在兩個倍頻值，分為兩種情況。第一種為兩個倍頻值是依次關系，如19、20，對應的諧共振中心頻率為78.95 Hz和75 Hz；第二種為兩個倍頻值相距較遠，如11、20，對應的諧共振中心頻率為136.36 Hz和75 Hz。激振力大小相同，同時傳感器安裝滿足2.2.1中的要求。對兩種情況下的兩個諧共振中心頻率附近進行恒速運轉及倍頻值分析。如19倍頻時恒速運轉轉頻范圍為78.5～79.5 Hz，步長為0.1 Hz(下文中轉頻步長同樣設為0.1 Hz)；20倍頻處進行的恒速運轉轉頻范圍為74～76 Hz。每一運轉轉頻通過AR法求得的倍頻值，如圖3所示。

11倍頻處進行的恒速運轉轉頻范圍為135～138 Hz，11倍頻和20倍頻共存時不同運轉轉頻下求取的倍頻結果,如圖4所示。

(1) 對上述兩種情況的分析可知，多倍頻存在時，通過AR法(適用AR的拓展法)可以分析出存在的葉片振動倍頻值；

(a)19倍頻(b)20倍頻

圖3 倍頻值為19和20時不同轉頻下求解的倍頻值變化情況

圖4 倍頻值為11和20時不同轉頻下求解的倍頻值變化情況

Fig.4 Obtained blade vibration EO at different operation speed for EO=11 and EO=20

(2) 在諧共振中心頻率附近進行運轉測量，對倍頻值的求取也是非常準確的，即不需要將運轉轉頻與諧共振中心頻率保持一致，在諧共振中心頻率附近也可以實現，這也與1.1節的結論對應。

2.2.3 測量過程中干擾(或噪聲)對倍頻值求取的影響

就實際測量系統而言，干擾或者噪聲的存在是不可避免的。通過使用GARIV法可減弱噪聲或干擾對測量結果的影響，提高測量精度。下面通過添加不同強度的高斯噪聲來比較AR、GAR、GARIV三種方法在不同運轉轉頻下的倍頻值分析結果。假設模型存在多倍頻19、20。添加噪聲情況為：無噪聲，信噪比為100，90，70四種情況。圖5為倍頻值20時不同噪聲下三種分析方法在不同運轉轉頻下的倍頻值分析結果。

分析可知無噪聲時，AR、GAR、GARIV三種方法獲得的結果一樣，如圖5(a)。隨著噪聲的逐漸增大，GAR法計算的結果出現較大的偏差，AR法和GARIV法與正確值接近，但AR法在計算過程中會由于噪聲的存在導致a1值超出允許范圍而失效，綜合比較可知GARIV法具有抗干擾能力，如圖5(b)和(c)。隨著噪聲強度的增大，即SNR減小到70時，三種方法的測量結果都出現較大的誤差，測量準確性降低，如圖5(d)。總結而言，GARIV法對一定噪聲下的測量分析結果較為準確，即具有較好的抗干擾能力，可用來分析實際測量中的振動倍頻值。通過比較相同噪聲，不同運轉轉頻下的結果可知，越接近諧共振中心頻率時，運轉測量結果的準確性越高，抗干擾能力越強。所以，可以通過將運轉轉頻接近諧共振中心頻率來提高倍頻值的測量精度，但也需要考慮到諧共振中心處的葉片振動局部最大，避免發生危險情況。

(a)無噪聲,不同運轉轉頻下倍頻值求解(b)SNR=100,不同運轉轉頻下倍頻值求解

(c)SNR=90,不同運轉轉頻下倍頻值求解(d)SNR=70,不同運轉轉頻下倍頻值求解

圖5 20倍頻，不同信噪比和運轉轉頻下倍頻值求解

Fig.5 Comparison of the results of blade vibration EO at different rotating speed and noise

2.2.4 傳感器安裝角度誤差對倍頻值求取的影響

針對安裝角度誤差，通過按一定范圍內生成隨機數，以達到產生隨機夾角誤差的目的。為避免特殊情況的發生，通過大量隨機產生的模擬誤差來研究安裝誤差對倍頻值求取的整體影響趨勢。

同時由于分析公式中對傳感器安裝夾角有要求，仿真中分別通過設計夾角和平均夾角兩種情況來對比哪種夾角下的分析結果更準確。平均夾角由各傳感器的夾角經實際測量并平均后獲得。圖6為倍頻值20，添加噪聲一定時，不同最大夾角誤差時的仿真結果對比。考慮到實際運轉轉頻不可能與諧共振中心頻率相同，故選取運轉轉頻為75.2 Hz。

(a) 夾角誤差最大0.5°，平均夾角下運行次數

(b) 夾角誤差最大0.5°，設計夾角下運行次數

(c) 夾角誤差最大1.5°，平均夾角下運行次數

(d) 夾角誤差最大1.5°，設計夾角下運行次數

通過圖6分析可知，相同噪聲下，采用平均夾角的計算結果要比采用設計夾角的計算結果更準確。同時隨著夾角誤差的增大，采用設計夾角的分析結果雖然較為平穩，但與正確值偏差增大。采用平均夾角的分析結果雖然與準確值接近，但也出現了較大的波動。

總結以上，在設計和安裝傳感器時，不能超出允許夾角的同時，可通過動態測量的方式獲得各傳感器之間的夾角，進而獲得平均安裝夾角。應保證各安裝夾角相等，且與設計值保持一致。

3 掃頻擬合分析與坎貝爾圖的獲取

3.1掃頻擬合分析

將旋轉葉片簡化成懸臂梁，不考慮葉片之間耦合時，葉片振動位移隨轉頻之間的響應函數可表示為

(13)

式中:A0為激振力幅值；Q為品質因數，Q=1/(2ξ);ξ表示阻尼比；A0Q表示振動幅值；ω為轉頻；ωn為中心轉速頻率；φ表示初始相位；DC為葉片振動恒偏量。η′表示如下：

(14)

當葉片運轉轉速經過諧共振中心轉頻時，葉片在激振力的作用下會產生較大的振動。根據振動位移與轉頻之間的響應函數，可通過非線性最小二乘擬合算法獲得諧共振中心頻率、振動幅值、初始相位、振動恒偏量等參數。非線性最小二乘擬合算法具體內容可參考LM法[9-10]。

在仿真實驗過程中發現，不同升速情況下，葉片振動相應曲線不同，如圖7所示。

圖7 不同升速情況下葉片振動位移響應曲線

響應曲線是擬合的數據源，對擬合結果影響很大，所以需通過仿真對升速值進行確定。可有效保證實際測量時葉片振動參數擬合的準確性。在上述仿真的基礎上，通過不同升速下的擬合結果來確定可允許的最大升速值。

圖8 不同傳感器在不同升速下的響應曲線擬合情況

通過圖8可知，在不同升速下葉片振動幅值的擬合結果不同，隨著角加速度的增大，擬合獲得的幅值逐漸減小，與實際值偏離逐漸增大，擬合的誤差逐漸增大，所以選擇合適的升速值可通過擬合獲得較為準確的葉片同步振動參數。

3.2綜合驗證

為綜合驗證掃頻擬合結合GARIV法來辨識葉片同步振動參數，進而獲得葉片振動坎貝爾圖。建立的仿真模型參數設置與2.1節一致，葉片運轉區間內被激起的倍頻值和振幅對應關系,如表1所示。

表1 模型參數設置

仿真時傳感器布置角度滿足前文中的要求，在振動位移的結果上添加10%的高斯白噪聲。

先通過掃頻運轉獲得葉片振動位移隨運轉轉速的變化情況，初步識別出存在的諧共振區間和粗糙的共振中心頻率，然后對提取出的諧共振區間進行非線性最小二乘擬合(LM法擬合)獲得葉片振動位移響應曲線的參數，包括諧共振中心頻率、振動幅值、相位、振動恒偏量等。對四支傳感器擬合的結果進行平均獲得葉片同步振動參數。根據掃頻結果針對存在的每一個諧共振區間選取合適的運行轉頻，然后進行恒速運轉仿真，四支葉頂傳感器監測到的振動位移經過GARIV方法分析可以獲得葉片各個諧共振中心對應的振動倍頻值。表2為通過掃頻擬合獲得的諧共振中心頻率和振動幅值，以及GARIV分析獲得的振動倍頻值(表中的頻率指諧共振中心頻率)。

表2 葉片振動參數辨識結果

通過上表可以發現LM法對葉片振動幅值和諧共振中心的擬合結果是非常準確的，GARIV方法也較準確的分析出各諧共振中心對應的振動倍頻值。通過上表可以獲得葉片振動坎貝爾圖，如圖9所示。

圖9 葉片振動坎貝爾圖

4 結 論

本文提出了掃頻擬合分析和GARIV法相結合來辨識葉片同步振動時的振動參數。主要完成以下：

(1) 通過建模仿真給出了不同倍頻下傳感器允許的安裝夾角，驗證了AR法及其擴展法可對單倍頻和多倍頻進行有效辨識，同時也對運轉轉速的選擇做了說明。

(2) 驗證了GARIV法對葉片振動倍頻識別的準確性，以及抗干擾能力；通過仿真對傳感器安裝夾角的誤差進行了規范，提出GARIV分析算法中計算采用傳感器平均夾角。

(3) 在使用掃頻擬合方法的過程中，通過比較不同升速值下的葉片振動參數擬合結果提出了合適的升速值。

通過結合兩種方法的優點可辨識出葉片同步振動參數，在準確辨識的基礎上將使用的傳感器數量減少到四個，可較好的滿足工程應用。

[1] PICKERING T M. Methods for validation of a turbomachinery rotor blade tip timing system[D]. Blacksburg: Blacksburg Virginia Tech University, 2014.

[2] ABDELRHMAN A M, HEE L M, LEONG M S, et al. Condition monitoring of blade in turbomachinery: A review[J]. Advances in Mechanical Engineering, 2014(2014)：210717.

[3] DIMITRIADIS G, CARRINGTON I B, WRIGHT J R, et al. Blade-tip timing measurement of synchronous vibrations of rotating bladed assemblies[J]. Mechanical Systems & Signal Processing, 2002, 16(4): 599-622.

[4] HEATH S. A new technique for identifying synchronous resonances using tip-timing[J]. Journal of Engineering for Gas Turbines and Power, 2000, 122: 219-225.

[5] 歐陽濤. 基于葉尖定時的旋轉葉片振動檢測及參數辨識技術[D]. 天津: 天津大學, 2011.

[6] ROBINSON W, WASHBURN R. A real time non-interference stress measurement system (NSMS) for determining aero engine blade stresses[C]//Proceedings of the 37th Symposium of the Instrumentation Society of America, San Diego, California. 1991: 91-103.

[7] CARRINGTON I B. Development of blade tip timing data analysis techniques[D]. Manchester: University of Manchester, 2002.

[8] 歐陽濤, 郭文力, 段發階, 等. 基于葉尖定時的旋轉葉片同步振動辨識新方法[J]. 振動與沖擊, 2011, 30(8): 249-252.

OUYANG Tao, GUO Wenli, DUAN Fajie, et al. New method for identifying rotating blades synchronous vibration based on tip-timing[J]. Journal of Vibration and Shock, 2011, 30(8): 249-252.

[9] FREUND O, MONTGOMERY M, MITTELBACH M, et al. Non-contact test set-up for aeroelasticity in a rotating turbomachine combining a novel acoustic excitation system with tip-timing[J]. Measurement Science and Technology, 2014, 25(3): 116-121.

[10] 張光澄. 非線性最優化計算方法[M]. 北京: 高等教育出版社， 2005： 159-169.

Synchronousvibrationparametricidentificationmethodforbladesofaturbine

WANG Weimin, REN Sanqun, CHEN Lifang, HUANG Shan

(Beijing Key Laboratory of Health Monitoring Control and Fault Self-recovery for High-end Machinery, Beijing University of Chemical Technology, Beijing 100029, China)

A new method used to identify a turbine’s blade synchronous vibration parameters based on the blade tip-timing (BTT) non-contact type blade vibration measurement technique and previous studies in this field was proposed here. According to the measured relationship between blade vibration displacement and rotating speed, the parameters, such as, the center frequency of blade’s harmonic resonance points, the maximum amplitude of blade vibration, and its initial phase were computed with the nonlinear least square fitting algorithm. If the blade assembly was operated at a constant rotating speed being close to the frequency of its resonance point, different blade vibration displacements were measured with sensors at different positions. The double-frequency vibration displacements of various blades at their corresponding harmonic resonance points were computed with the global autoregressive with instrumental variables (GARIV) method. A blade’s Campbell diagram was made with its vibration parameters obtained using GARIV method and the sweep frequency fitting method. Based on this method, the locating scheme of BTT sensors was proposed, interference factors existing in the actual measurements were analyzed in detail, the feasibility of the proposed method was verified. The test verification of this method will be published later.

turbine blades; synchronous vibration; parametric identification; blade tip-timing (BTT); GARIV method

國家自然科學基金(51275028; 51135001)

2016-11-01 修改稿收到日期：2016-12-16

王維民 男，博士，教授，1978年4月生

陳立芳 女，博士，講師，1973年3月生

TK14； V216.2

： A

10.13465/j.cnki.jvs.2017.17.019