汽輪機長葉片顫振預測方法

劉長春，關淳，郭魁俊，李宇峰，馬義良

汽輪機長葉片顫振預測方法

劉長春，關淳，郭魁俊，李宇峰，馬義良

（哈爾濱汽輪機廠有限責任公司，黑龍江省 哈爾濱市 150046）

基于葉片氣動和結構振動仿真計算平臺，建立了汽輪機長葉片流固耦合的三維計算模型。以國外某60Hz汽輪機機組次末級葉片為例，分析了葉片流場的氣動特性和整圈葉片結構場的振動特性，并進一步結合能量法準則，提出了汽輪機長葉片顫振預測和評估的方法。結果表明，機組流量和葉片節徑振動形式對葉片顫振特性影響很大。對于當前次末級葉片，40%流量工況時，流場內壓力脈動誘導葉片顫振發生的概率為80%～90%，因此在小流量工況下需慎重使用該類葉片。基于單向流固耦合的方法，能夠有效且相對快速地分析流場中脈動壓力作用下的周期功，獲得葉片模態振動條件下的氣動阻尼系數，從而預測葉片顫振特性。

汽輪機；流固耦合；氣動特性；振動特性；顫振預測

0 引言

大型汽輪機葉片的工作條件非常惡劣，葉片事故也時有發生，葉片事故是造成機組強迫停機的主要原因之一。據不完全統計，運行中葉片事故約占汽輪機事故的40%，其中60%～80%葉片損壞與振動疲勞失效導致的斷裂有關，而葉片的檢修費用約占汽輪機檢修費用的50.5%。隨著大功率汽輪機的發展，汽輪機低壓缸末幾級葉片的展弦比越來越大，剛度越來越小，其顫振問題也因此愈發突顯，尤其是對于那些需要頻繁在小負荷、高背壓工況下運行的大型汽輪機，葉片顫振問題逐漸引起關注[1-4]。

葉片顫振是一種由流體誘發的振動，是流體和結構相互作用產生的不穩定現象。對葉片的顫振特性的研究涉及到流體動力學和結構動力學2個方面。關于汽輪機長葉片顫振的分析，目前國外自由葉片的計算及試驗考核體系基本完善。但在國內，顫振分析仍主要集中在壓氣機和航空發動機葉片，關于電廠汽輪機葉片顫振分析方法的研究仍然存在很大不足，尤其是對于整圈圍帶連接的長葉片，尚缺少業內公認的顫振考核標準和評估體系。

目前，葉片顫振問題的研究方法主要有兩大類：一是基于半經驗的變形激盤法[4-6]，這種顫振預測的方法已被成功應用于航空領域；二是數值計算方法，主要包括能量法和基于流固耦合的時域分析法。其中，能量法是根據葉片模態分析的結果，針對特定的葉片振動頻率、振幅和振型，通過分析一個振動周期內流體對葉片做功的大小來分析預測葉片的顫振[7-8]。近幾年來，計算機技術的快速發展使得計算流體/固體力學耦合計算方法也得到了迅速的發展?；诹鞴恬詈系臅r域分析法能夠同時考慮流體的氣動阻尼和葉片的機械阻尼，更加真實地反映流場與結構場之間的相互作用，已被應用于葉片顫振研究方面[9-11]。時域分析的耦合方法分為全耦合和離散耦合2種。采用全耦合方法時，流體和固體方程在一個統一的矩陣中求解，這種方法需要的假設條件少，更接近于實際過程，但是求解矩陣的建立非常困難，并且計算量巨大，計算收斂困難，對于流動和結構的耦合只能應用于一些非常簡單的研究。離散耦合法將耦合系統分解成單獨的子系統，用傳統方法逐一求解各子系統，利用流固交界面在子系統之間傳遞壓力和位移等耦合信息，通過迭代使整個系統達到平衡[12]。目前，幾乎所有的流固耦合分析都采用離散耦合方法。

由此可知，得益于計算機技術的快速發展，流固耦合方法已逐漸成為葉片氣彈性耦合計算的主流方法，并有望成為葉片顫振特性研究的重要手段。但目前關于汽輪機葉片顫振特性流固耦合方法的研究依然很少，尤其是對于整圈存在圍帶連接的葉片組，其顫振特性流固耦合方法鮮見報道。為此，本文基于葉片氣動和結構振動仿真計算平臺，建立了汽輪機長葉片流固耦合的三維計算模型，分析了葉片流道內流場的氣動特性和整圈葉片結構場的振動特性，并進一步結合能量法準則，提出了汽輪機長葉片顫振預測和評估的方法，有利于在汽輪機長葉片前期的開發和設計中降低葉片顫振的風險，為汽輪機長葉片的安全使用提供依據。

1 計算模型和方法

1.1 計算模型

本文以國外某60Hz汽輪機次末級動葉為研究對象，表1為該葉片、圍帶及輪槽等部件的材料參數和機組在該級的設計參數。基于葉片氣動和結構振動仿真計算平臺，分別建立了氣動分析子系統和結構分析子系統，并以葉片壓力面和吸力面為流固耦合面，將氣動分析壓力場結果直接傳遞給結構分析子系統，實現了流固耦合分析系統的構建。為了確保2個子系統耦合面之間數據傳輸準確，需要保證流場模型和結構場模型節點坐標保持一致。為此，本文在計算模型構建過程中，首先將葉片型線文件導入網格劃分模塊中，獲得流場計算域模型，而后經網格處理模塊實現葉片實體和流場計算域之間的分離，進一步將葉片實體導入UG軟件中進行二次模型修正，最終得到葉片結構場的計算域模型。通過這樣的方式建立模型，能夠保證在整個模型建立過程中葉片壓力面和吸力面的坐標始終保持不變，從而確保了耦合面之間數據映射準確。流固耦合分析流程如圖1所示。

表1 轉子材料和設計參數

圖1 流固耦合分析流程圖

氣動分析子系統采用計算流體力學(computational fluid dynamics，CFD)模塊進行三維、可壓縮、單相、黏性流場分析，湍流模型選取模型，邊界層采用可擴展壁面函數進行處理，壁面光滑、無摩擦，工質選用濕蒸汽，工質參數基于IAPWS IF97數據庫獲得。本文選擇2個工況點作為研究工況，分別為100%流量和40%流量，2個工況點的邊界參數見表2。對于100%流量工況，入口邊界采用Inlet模式，條件限定為壓力和溫度，出口邊界采用Outlet模式；但對于40%流量工況，考慮到小流量可能會引起汽流回流，入口和出口模式均改為Opening模式，使計算具備捕捉汽流回流的功能。流道網格為六面體網格，采用雙流道模型作為計算域，圖2為計算域網格劃分，經網格無關性驗證，確定單流道網格單元數約為320000。

表2 計算域邊界條件

結構分析子系統采用有限元仿真計算模塊計算應力特性和模態特性，設置循環對稱邊界條件，圍帶接觸面采用摩擦接觸，摩擦因數為0.2，接觸面剛度為1N/mm3，滲透限制為0.1mm，輪槽底部,,三個方向位移約束，轉速設置和氣動分析子系統保持一致，為3600r/min，模態分析過程忽略結構阻尼的影響，分析內容為第一階的所有節徑振動。結構場計算域網格通過映射技術生成的六面體網格，以一階M0頻率為目標，經網格無關性驗證，單只葉片網格單元數確定為8000。

1.2 計算方法

由圖1可知，基于已有的流場和結構場模型，可進一步實現雙向流固耦合和單向流固耦合分析。雙向流固耦合分析通過在子系統之間傳遞網格位移和表面壓力，實現流體域和固體域的雙向耦合[10]，這種方法所需要的假設條件少，更接近于實際過程，能夠獲取流場和結構場相互作用的時域結果。但對于存在節徑振動的整圈圍帶連接的葉片組，雙向流固耦合方法需要建立更多的流道模型才能獲得實際的節徑振動，這就會大大增加計算量[8,13]。相比而言，采用單向流固耦合方法可以顯著減少計算量，縮短計算周期，并能夠采用循環對稱的邊界條件獲得整圈節徑振動結果。但需要注意的是，采用單向流固耦合方法存在假設條件，如假定葉片做恒幅周期振動，振動頻率等于某一節徑頻率，這使得氣彈性作用的計算結果與真實情況存在一定偏差，對顫振特性的預測將更加保守[14]。表3匯總了2種耦合方法的對比結果?？紤]到本文研究的模態振動為整圈帶有圍帶的葉片組振動，雙向流固耦合方法對計算機的配置要求將會非常高，并且所耗費的時間和內存也是難以接受的，因此，本文出于工程應用的考慮，選用單向流固耦合方法，基于流場壓力的計算結果，分析結構場模態振動特性，而后將結構場分析的各個節徑振動下的振動位移導回流場分析子系統中，借助CFD軟件中顫振計算模塊，分析流場中脈動壓力作用下的周期功，最終獲得當前工況和模態振動條件下的氣動阻尼系數，以此來評估葉片的顫振特性。

圖2 計算域網格劃分

表3 2種評估顫振的流固耦合方法對比

2 計算結果與分析

2.1 流場氣動特性

圖3為2個研究工況條件下流道典型截面的壓力包絡線，其反映出葉片表面壓力沿流動方向的變化。典型截面包括50%葉高處(span=0.5)和85%葉高處(span=0.85)。由圖3(a)、(b)可以看出，在100%流量工況下，壓力面和吸力面的壓力沿流動方向逐漸降低，入口附近吸力面的壓力降低更為明顯，壓力面壓力能夠保持恒定。葉片出口附近，壓力面壓力開始迅速下降，在出口圓弧處，2個面的壓力逐漸趨于相等，但存在一定的波動，這表明當前工況下汽流在出口處并不能很好地過渡，壓力的波動可能會導致汽流尾跡，并造成端部損失。對比不同特征截面壓力包絡線的結果可以看出，100%流量工況條件下，吸力面和壓力面所組成的壓力包絡線的面積沿葉高方向逐漸增大，這表明從葉根到葉尖，葉片的負荷逐漸增大，做功能力逐漸增強[10]。

圖3 流道典型截面壓力包絡線圖

此外，由圖3(c)、(d)可以看出，在40%流量工況下，壓力面壓力較低，而吸力面的壓力變化非常劇烈，這是因為小流量工況會在入口處造成較大的負攻角，汽流沖擊葉片吸力面劇烈，而壓力面會造成漩渦甚至脫流區。40%流量工況下，2個截面壓力包絡線的面積相差不明顯，均存在負功區域，并且面積明顯小于100%工況下的包絡線面積，這說明蒸汽整體做功的能力嚴重降低，級效率差。

圖4為流道典型截面流線圖，其可以反映出流場沿軸向和周向的變化。由圖4(a)、(b)可以看出，100%流量工況下，流道內的流線較為規整，不存在明顯的渦流，速度變化與圖3的壓力變化相對應。而由圖4(c)、(d)可以看出，40%流量工況下，入口處汽流存在較大的負攻角，汽流沖擊葉片吸力面，導致壓力面形成徑向渦流，吸力面形成高速區，這些現象與壓力包絡線的變化相符合。

圖4 流道典型截面流線圖

圖5為流道子午面流線，其可以反映出流場沿軸向和徑向的變化。由圖5可知，100%流量工況下，流場整體分布均勻，沿流動方向速度逐漸增大。而40%流量工況下，流道出口出現汽流倒流現象，流場在葉根附近也相應出現渦流，并且渦流的方向與典型截面所看到的方向有所不同。典型截面所看到的渦流中心線沿徑向(徑向渦流)，這是由進口負攻角所導致的，而子午面所看到的渦流中心線沿周向(周向渦流)，通常是由出口汽流倒流所致，在小流量、高背壓、鼓風等工況下容易出現。由于小流量工況下往往都導致較大的負攻角，因此流道內周向渦流的出現往往都會伴隨著徑向渦流。

圖5 流道子午面流線圖

2.2 結構場振動特性

葉片振動特性的研究是葉片顫振特性研究的前提，只有獲得葉片振動頻率和振型，才能夠進一步計算葉片表面脈動壓力對葉片所做的周期功，以及流固耦合分析系統在當前工況和模態振動條件下的氣動阻尼，從而實現葉片顫振分析。葉片表面壓力分布是葉片振動特性分析的一個重要邊界條件，本文基于單向流固耦合方法，將流場分析子系統中葉片表面的壓力分布結果通過耦合面直接導入到葉片結構分析子系統中。流場和結構場耦合面之間壓力映射的結果見表4?？梢钥闯?，壓力數據的傳遞效率達到了100%，證實了1.1節中模型建立方法的有效性和數據傳遞的準確性。

表4 流固耦合面壓力映射結果

對于圍帶連接的整圈葉片，葉片振型為輪盤的節徑振動。當前研究的動葉一共有68只，因此整圈葉片每一階將存在34個前行波節徑振動和34個后行波節徑振動。其中，前行波由正節徑表示(ND>0)，后行波由負節徑(ND<0)表示?？紤]到葉片顫振主要發生在低頻區域，通常低階模態振動的影響較大[10,15]，為此，本文的模態分析過程包括2個工況下整圈葉片第一階的所有節徑振動頻率和振型，并隨機提取了一階5個節徑振動下的振動頻率和振型，用于進一步計算葉片氣動阻尼。表5列舉了所提取的5個節徑振動的靜頻和在離心力/汽流力作用下的動頻，表6為這5個節徑振動對應的振型。通過對比100%流量和40%流量工況的振型圖可以看出，2種工況下各個節徑振動的振型一致，但振動幅值存在差別。100%流量工況下的振動幅值大于40%流量工況下的幅值，這主要是因為當前所做的模態分析是基于穩態汽流力的定常分析結果，小流量工況下流道內的汽流發分布不均，局部渦流會導致葉片表面壓力降低，因而葉片振動能量較低，振動幅值較小。但這并不能說明小流量工況下葉片振動較弱，小流量工況下的渦流會導致葉片表面存在脈動壓力，當脈動壓力的頻率與葉片振動固有頻率一致時，反而會造成葉片共振，并且小流量工況下流場與葉片之間的氣彈性耦合更加復雜，往往還會強化葉片的自激振動，能量逐漸積累，振幅越來越大，從而發生顫振。為此，本將采用非定常分析方法進一步對比研究40%流量工況下葉片的顫振特性。

表5 葉片一階模態振動頻率

表6 葉片一階模態振動振型

2.3 葉片顫振特性

在不同流量工況下模態振動的節徑、頻率、振型以及工況邊界條件，可以通過CFD顫振分析模塊計算葉片在當前流量工況下的周期功和響應的氣動阻尼系數。圖6為40%流量工況各個節徑振動下的氣動阻尼系數、阻尼系數擬合線和部分計算點的周期功。周期功的計算值為一個振動周期內非定常平均積累功?；跉鈩幼枘岬恼?，本文將葉片某階振動下顫振發生的概率定義為氣動阻尼為正的節徑數占該階振動總節徑數的百分比。當前葉片一階振動最小氣動阻尼系數為ND=?17，即葉間相位角為?90°，這與文獻[9]基于時域分析方法得到的顫振預測結果一致。該工況下氣動阻尼系數擬合線表達式為

式中daero為系統氣動阻尼系數。

擬合線的擬合系數為2=0.999。由當前計算結果可以看出，若不考慮結構阻尼，在40%流量工況下，約90%節徑振動的氣動阻尼系數為負數，表明葉片在這些節徑振動下會從汽流中獲取能量，即汽流的非定常氣動力對葉片做正功，如此下去，葉片的振動會持續加劇，進而激發顫振。這類顫振主要是由動葉入口較大的負攻角而引發的失速顫振。在小流量工況下，汽流在入口處以較大負攻角流向動葉，在動葉壓力面形成徑向渦流和周向渦流，構成失速脫流區，引起葉片表面氣動力的改變，造成葉片彈性變形，進而又再次影響汽流攻角和繞流狀態，形成不穩定流場與葉片彈性變形之間的氣動耦合，當耦合作用使流場的能量不斷地輸入葉片系統中時，就會引發葉片顫振，導致葉片振幅顯著增大[16]。對于汽輪機末幾級長葉片，結構阻尼系數通常為0.01～0.02。因此，若考慮結構阻尼，系統總阻尼系數會有所增加，但即便如此，對于本文研究的次末級葉片，當處于40%流量工況下時，一階振動誘導顫振發生的概率依然高達80%，因而在非設計工況下，該類葉片需要慎重使用。

3 結論

基于葉片氣動和結構振動仿真計算平臺，建立了汽輪機長葉片流固耦合的三維計算模型，分析了葉片流道內流場的氣動特性和整圈葉片結構場的振動特性，并進一步結合能量法準則，提出了汽輪機長葉片顫振預測和評估的方法。主要結論如下：

1）流量變化會對汽輪機流道的氣動特性產生顯著的影響。小流量工況下會同時存在徑向渦流和周向渦流，其形成分別與入口負攻角和出口回流有關。

2）流量變化會對汽輪機長葉片模態振動特性產生影響。不同工況下葉片各個節徑振動的振型一致，但振動頻率和幅值存在差別。小流量下流道內存在局部渦流，導致葉片表面壓力降低，葉片振動能量較低，振動幅值較小。

3）流量工況和葉片節徑振動形式對葉片顫振特性影響很大。對于所研究的次末級葉片，當處于40%流量工況時，一階各節徑振動誘導顫振發生的概率為80%～90%，因此在小流量工況下需要慎重使用該類葉片。

4）基于單向流固耦合方法，將結構場模態振動特性和流場氣動特性相結合，能夠有效且相對快速地分析流場中脈動壓力作用下的周期功，獲得葉片模態振動條件下的氣動阻尼系數，實現葉片顫振特性的評估。

[1] 李宇峰，任大康，黃鋼，等．空冷汽輪機低壓末級變工況設計[J]．熱力透平，2004，33(1)：14-16．

LI Y F，REN D K，HUANG G，et al．The off-design for LP last stage blade of air-cooling units[J]．Thermal Turbine，2004，33(1)：14-16．

[2] 劉萬琨．汽輪機末級葉片顫振設計[J]．東方電氣評論，2007，21(4)：7-13．

LIU W K．Flutter design for steam turbine last blade[J]．Dongfang Electric Review，2007，21(4)：7-13．

[3] 陶德平，楊曉東，周盛．蒸汽輪機長葉片顫振預估方法研究[J]．航空動力學報，1991(2)：151-156．

TAO D P，YANG X D，ZHOU S．Prediction of long blade flutter in a steam turbine[J]．Journal of Aerospace Power，1991(2)：151-156．

[4] 張揚軍，陶德平．汽輪機葉片顫振研究的新進展[J]．科學通報，1996(23)：2204-2206．

ZHANG Y J，TAO D P．New progress in research on blade flutter of steam turbine[J]．Chinese Science Bulletin，1996(23)：2204-2206．

[5] 張揚軍，陶德平．變葉片間相角的蒸汽輪機葉片顫振預測方法[J]．水利電力機械，1993(4)：29-32．

ZHANG Y J，TAO D P．Prediction method of steam turbine blade flutter with variable blade phase angle

[J]．Hydraulic and Electric Machinery，1993(4)：29-32．

[6] 張揚軍，李克儉，陶德平．葉片間相角對蒸汽輪機葉片顫振的影響[J]．航空動力學報，1994(3)：277-280．

ZHANG Y J，LI K J，TAO D P．Effect of interblade phase angle on blade flutter of steam turbine[J]．Journal of Aerospace Power，1994(3)：277-280．

[7] 張興國，劉鋒，張陳安．基于CFD技術的葉片顫振分析[J]．航空計算技術，2009，39(4)：75-78．

ZHANG X G，LIU F，ZHANG C A．Flutter computation of turbo machinery blades based on CFD

[J]．Aeronautical Computing Technique，2009，39(4)：75-78．

[8] CARTA F O．Coupled blade-disk-shroud flutter instabilities in turbo-jet engine rotors[J]．Journal of Engineering for Gas Turbines and Power，2011，89(3)：419．

[9] 姜偉，謝誕梅，陳暢，等．基于時域分析法的汽輪機末級葉片顫振預測及分析[J]．振動與沖擊，2015，34(11)：194-199．

JIANG W，XIE D M，CHEN C，et al．Flutter prediction and analysis for a steam turbine last-stage blade based on time domain analysis method[J]．Journal of Vibration and Shock，2015，34(11)：194-199．

[10] 張帥，高麗敏，鄭天龍，等．基于雙向流固耦合方法的某風扇特性數值研究[J]．工程熱物理學報，2017，38(8)：1683-1691．

ZHANG S，GAO L M，ZHENG T L，et al．Numerical investigations on characteristic of fan based on two-way fluid-structure interaction approach[J]．Journal of Engineering Thermophysics，2017，38(8)：1683-1691．

[11] 鄭赟，楊慧．跨音速風扇全環葉片顫振特性的流固耦合分析[J]．北京航空航天大學學報，2013，39(5)：626-630．

ZHENG Y，YANG H．Full assembly fluid/structured flutter analysis of a transonic fan[J]．Journal of Beijing University of Aeronautics and Astronautics，2013，39(5)：626-630．

[12] SADEGHI M，LIU F．Coupled fluid-structure simulation for turbomachinary blade rows[EB/OL]．

[2021-01-01]．https://arc.aiaa.org/na101/home/literatum/publisher/aiaa/books/content/6.asm/2005/masm05/6.2005-18/staging/6.2005-18.fp.png_v01.

[13] 杜云祥，徐自力，焦玉雪，等．基于振動時滯法的非零葉間相位角葉片顫振[J]．航空動力學報，2020，35(8)：189-198．

DU Y X，XU Z L，JIAO Y X，et al．Blades flutter of non-zero inter-blade phase angle based on vibration time-delay method[J]．Journal of Aerospace Power，2020，35(8)：189-198．

[14] 周迪，呂彬彬，陸志良，等．能量法和時域法在葉片顫振計算中的比較研究[J]．航空計算技術，2019，49(5)：43-48．

ZHOU D，LV B B，LU Z L，et al．Comparative study on energy method and time-domain method for blade flutter predictions[J]．Aeronautical Computing Technique，2019，49(5)：43-48．

[15] 李迪，張曉杰，王延榮．壓氣機轉子葉片的抑顫設計[J]．推進技術，2020，41(9)：207-216．

LI D，ZHANG X J，WANG Y R．Design for flutter suppression of rotor blade in a compressor[J]．Journal of Propulsion Technology，2020，41(9)：207-216．

[16] 楊光海．汽輪機葉片的安全防護[M]．北京：機械工業出版社，1992．

YANG G H．Safety protection of steam turbine blade[M]．Beijing：China Machine Press，1992．

Flutter Prediction Method for Long Blade of Steam Turbine

LIU Changchun, GUAN Chun, GUO Kuijun, LI Yufeng, MA Yiliang

(Harbin Turbine Company Limited, Harbin 150046, Heilongjiang Province, China)

A three-dimensional fluid-solid coupling calculation model of long blade of steam turbine was established based on the aerodynamic and structural vibration simulation platform. Taking the second last stage blade of a foreign 60Hz steam turbine unit as an example, the aerodynamic characteristics of the flow field in the single flow channel and the vibration characteristics of the entire blade structure field were analyzed. Combined with the energy method criterion, a method for predicting and evaluating the flutter of the steam turbine long blade was further proposed. The results show that the unit flow rate and the vibration form of blade pitch diameter have a great influence on the blade flutter characteristics. For the current 2nd last stage blade, the probability of blade flutter induced by pressure fluctuation in the flow field is between 80% and 90% under 40% flow condition. Therefore, this kind of blade should be carefully used under the low flow condition. Based on the unidirectional liquid-solid coupling method, the periodic work under fluctuating pressure in flow field can be analyzed effectively and relatively quickly, and the aerodynamic damping coefficient under modal vibration of blade can be obtained. Thus, the flutter characteristics of blade can be predicted.

steam turbine; liquid-solid coupling; aerodynamic characteristics; vibration characteristics; flutter prediction

2021-01-19。

10.12096/j.2096-4528.pgt.21001

TK 05

(責任編輯 尚彩娟)