姜　偉，謝誕梅，陳　暢，胡鵬飛，高　尚

(武漢大學(xué)動力與機械學(xué)院， 武漢　430072)

第一作者姜偉男，碩士生，1991年生

基于時域分析法的汽輪機末級葉片顫振預(yù)測及分析

姜偉，謝誕梅，陳暢，胡鵬飛，高尚

(武漢大學(xué)動力與機械學(xué)院， 武漢430072)

摘要：利用基于流固耦合的時域分析法對某汽輪機末級葉片的顫振問題進行了數(shù)值模擬研究。通過對不同工況、不同相位差下葉片振動響應(yīng)的計算，分析預(yù)測了顫振發(fā)生的臨界流量以及葉片的振型；通過對葉片受力與位移的分析，獲得了不同葉高處流體對葉片的做功，確定了葉片上氣動不穩(wěn)定區(qū)域。通過對氣動不穩(wěn)定區(qū)壓力分布的分析，解釋了葉片間相位差影響顫振發(fā)作的機理。此外，研究還發(fā)現(xiàn)葉片受力與位移之間的相位差是影響葉片氣動穩(wěn)定性的重要因素。

關(guān)鍵詞：時域分析法；末級葉片；顫振；流固耦合；氣動穩(wěn)定性

基金項目：國家自然科學(xué)基金項目(51376140)

收稿日期：2014-02-12修改稿收到日期：2014-05-23

通信作者謝誕梅女，博士，教授，1962年生

中圖分類號:Tp62

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2015.11.034

Abstract:A numerical simulation study was performed for a steam turbine last-stage blade with the time domain analysis method based on fluid-structure interaction. The critical flow to cause flutter and the vibration modal shapes of the blade were predicted by analyzing blade vibration responses under various operating conditions and interblade phase angles. Through analyzing load and displacement of the blade, the work done by fluid on the blade was obtained and the aerodynamic instability region on the blade was determined. The mechanism of interblade phase angle affecting flutter was illustrated through the analysis of pressure distribution in the aerodynamic instability region. The study also found that the phase difference between load and displacement of the blade is a key factor affecting its aerodynamic stability.

Flutter prediction and analysis for a steam turbine last-stage blade based on time domain analysis method

JIANGWei,XIEDan-mei,CHENChang,HUPeng-fei,GAOShang(School of Power and Mechanical Engineering, Wuhan University, WuHan 430072, China)

Key words:time domain analysis method; last-stage blade (LSB); flutter; fluid-structure interaction; aerodynamic stability

隨著大功率汽輪機的發(fā)展，汽輪機級末葉片的展弦比越來越大，其顫振問題變得越來越突出，特別對于那些不得不頻繁在低負荷工況下運行的大型汽輪機來說，葉片顫振問題越來越受到技術(shù)人員的關(guān)注[1-3]。

目前，研究葉輪機械葉片顫振問題的方法主要有兩大類，一是基于半經(jīng)驗的變形激盤法[4-5]，這種顫振預(yù)測的方法已被成功用于航空領(lǐng)域；二是數(shù)值方法，這類方法中又可細分為兩種，①能量法，其基本思路是根據(jù)葉片模態(tài)分析的結(jié)果，針對特定的葉片振動頻率、振幅和振型，通過考察一個振動周期內(nèi)流體對葉片做功的大小來分析預(yù)測葉片的顫振，文獻[6-8]采用了這種方法研究葉片顫振，但僅考慮氣動阻尼；②基于流固耦合的時域分析法，文獻[9-11]采用這種方法研究葉片顫振時，同時考慮了流體的氣動阻尼和葉片的機械阻尼。時域分析的耦合方法有全耦合(直接耦合)和離散耦合(間接耦合)兩種[12-13]。采用全耦合方法時，流體和固體方程在一個統(tǒng)一的矩陣中求解，這種方法需要的假設(shè)條件少，更接近于實際過程，但是求解矩陣的建立非常困難，并且計算量巨大。離散耦合法將耦合系統(tǒng)分解成單獨的子系統(tǒng)，用傳統(tǒng)方法逐一求解各子系統(tǒng)，在子系統(tǒng)之間傳遞壓力和位移等耦合信息，通過迭代使整個系統(tǒng)達到平衡。本文利用離散耦合的時域分析法，研究某600MW汽輪機末級葉片的顫振，預(yù)測其顫振邊界并分析顫振的機理及主要影響因素。

1數(shù)值模型及方法

(a)　流道　　　　(b)　葉片圖1　流道和葉片的網(wǎng)格模型Fig.1 Mesh model of flow pass and blade

本文的研究對象為某600MW汽輪機末級葉片，其葉高為1 016 mm，平均直徑為2 743 mm，葉片只數(shù)為94。利用ANSYS軟件分別構(gòu)建用于該葉片的流道氣動分析和葉片強度分析的子系統(tǒng)。通過在子系統(tǒng)之間傳遞網(wǎng)格位移和表面壓力實現(xiàn)流體域和固體域的雙向耦合。圖1(a)和圖1(b)分別給出了流道和葉片的網(wǎng)格模型。

圖2　數(shù)值方法的迭代層次Fig.2 Iteration levels of numerical method

求解程序的迭代關(guān)系見圖2，在最內(nèi)層的迭代中，子系統(tǒng)分別求解各自的計算域，當計算收斂后進入第二層迭代——“耦合迭代層”，在這一層迭代中通過反距離插值算法[14]由CFD(Computational Fluid Dynamics)求解器向CSD(Computational Structure Dynamics)求解器傳遞流體壓力，同時利用網(wǎng)格變形技術(shù)從CSD求解器向CFD求解器傳遞結(jié)構(gòu)變形數(shù)據(jù)，第三層迭代針對瞬態(tài)問題，可以獲得葉片受力和位移隨時間的變化，從而對顫振問題進行時域上的分析。

在氣動分析子系統(tǒng)中，采用ANSYS CFX模塊進行三維、可壓縮、單相、粘性流場分析，湍流模型選取k-ε模型，邊界層使用可擴展壁面函數(shù)處理。入口邊界為質(zhì)量流量，出口邊界為壓力出口，不同工況下的流量和背壓見表1。流道網(wǎng)格為ANSYS TurboGrid生成的六面體網(wǎng)格，單流道網(wǎng)格單元數(shù)約為40萬。

表1　不同工況下流量和背壓

在結(jié)構(gòu)分析子系統(tǒng)中，使用ANSYS Transient Structural模塊進行葉片動強度分析。葉片底邊剛性固定，表面受到流體壓力。采用四面體網(wǎng)格，葉底前后緣應(yīng)力集中處進行加密處理。

2葉片顫振預(yù)測及分析

圖3　流道子午面視圖Fig.3 Meridian plane view of flow pass

對于所研究的葉片，位移測點布置在葉片上50%～100%葉高的前緣和后緣處(見圖3)。測點命名為BxLEy和BxTEy(x：1、2，指定葉片；LE、TE：指定前緣、后緣；y：50、60、70、80、90、100，指定葉高)。

2.1顫振預(yù)測

由于該汽輪機末級葉片的數(shù)量較多，對全環(huán)葉柵建模需要10兆級別的網(wǎng)格數(shù)。為了減少網(wǎng)格數(shù)量，通常只會選取部分相鄰流道并設(shè)置周期性邊界條件。假設(shè)選取的流道數(shù)為n，全環(huán)葉片數(shù)為N，則該模型可以捕捉到的葉片間相位差α=2πk/n(k=0,1,2,…),其中要求N能被n整除。本文分別選取n=1、2，當n=1時，葉片間無相位差；當n=2時，可以捕捉到α=0、π兩種相位差。

圖4和圖5分別為測點B1LE100在n為1和2的情況下測得的位移數(shù)據(jù)。數(shù)據(jù)表明不考慮葉片間相位差影響時，在30%～100%THA工況(熱耗率驗收工況)下，葉片的振動響應(yīng)呈收斂趨勢，即可判別此時不會發(fā)生顫振；而當考慮葉片間相位角影響時，在45%～100%THA工況下，葉片不會顫振，在30%～45%THA工況下，葉片的振動響應(yīng)呈發(fā)散趨勢，會發(fā)生顫振，進而可以預(yù)測出該機組末級葉片顫振的臨界工況在43%～45%THA工況之間。

圖4　n=時不同工況下葉片的振動響應(yīng)Fig.4 Vibration response of blade under different work condition when n=1

圖5　n=2時不同工況下葉片的振動響應(yīng)Fig.5 Vibration response of blade under different work condition when n=2

該汽輪機末級葉片的模態(tài)分析結(jié)果見表2。由表2可知該葉片的一階彎振頻率比一階扭振頻率小得多，意味著其扭轉(zhuǎn)剛度比彎曲剛度大得多。圖6為測點B2LE100和B2TE100在43%THA工況下(如不加說明，后文所列數(shù)據(jù)均在該工況下測得)測得的部分振動響應(yīng)，該圖表明葉頂前緣和后緣測點間的振動相位差為零，即可判斷此時的葉片的振型為彎振。而對于航空領(lǐng)域的葉輪機械，比如壓氣機，其振型多以扭振為主[15-16]，相較于汽輪機末級葉片，這種振型上的差別主要是由展弦比的差異導(dǎo)致的。對該振動響應(yīng)進行傅里葉變換可以獲得葉片的振動頻率為40 Hz，接近于一階彎振頻率。圖7為測點B1LE100和B2LE100測得的振動響應(yīng)，該圖表明相鄰兩葉片相同位置測點的振動相位差為π，進一步表明葉間相位差是影響顫振的重要因素[17-18]。需要說明的是，本文并未考慮拉筋和圍帶的影響，理論上拉筋和圍帶可以在一定程度上控制葉間相位差，并增大葉片阻尼，抑制顫振的發(fā)生。

表2　葉片模態(tài)分析結(jié)果

圖6　葉頂前緣和后緣振動響應(yīng)Fig.6 Vibration response at leading point and trailing point of blade top

圖7　相鄰葉片葉頂前緣振動響應(yīng)Fig.7 Vibration response at leading point of adjacent blades

2.2壓力分布

計算中，在顫振振幅達到最大值后，本文選取一個振動周期提取壓力數(shù)據(jù)，周期用T表示。圖8為不同葉高處在一個周期內(nèi)的壓力分布，橫坐標中0到1為壓力面的壓力分布，1到2為吸力面的壓力分布。從圖中可以看出，自80%葉高處開始出現(xiàn)了比較明顯的壓力波動，且壓力的波動集中在壓力面0.5 -0.75處。前文的分析表明汽輪機末級葉片發(fā)生顫振時，葉片間存在相位差且葉片振型為彎振，這兩個條件導(dǎo)致葉片間的節(jié)距發(fā)生變化，從而導(dǎo)致葉片表面壓力的波動。

圖8　葉片表面壓力分布Fig.8 Pressure distribution at blade surface

圖9為80%葉高處某一周期內(nèi)t=0和t=T/2時刻的壓力分布。圖中實線表示t=0時的壓力場，虛線表示t=T/2時的壓力場。圖9顯示，在葉片的前緣處有密集的等壓線，意味著此處有非常大的壓力梯度，即存在激波，且該激波屬于脫體激波[19]。這道激波與相鄰葉片相交于0.5壓力面處，導(dǎo)致葉片壓力面壓力分布在0.5處驟增，圖8中80%葉高處的壓力分布以及圖9中3、4、5等壓線都表明了壓力驟增的過程。比較圖9中兩個時刻3、4、5等壓線可以發(fā)現(xiàn)，T/2時刻激波與葉片的相交位置前移了。激波作用于葉片位置的改變，導(dǎo)致80%葉高壓力面0.5～0.75處的壓力波動。激波與葉片相交位置之所以發(fā)生變化，一是因為葉片節(jié)距的改變，二是因為激波與葉片前緣距離的變化。葉片表面壓力的波動形成了葉片顫振發(fā)作的激振力。

圖9　不同時刻80%葉高處壓力分布圖Fig.9 Pressure distribution of 80% blade height at different time

2.3氣動不穩(wěn)定區(qū)

激振力的存在并非顫振發(fā)生的充分條件。下面通過對葉片受力與位移的分析，找出葉片上的氣動不穩(wěn)定區(qū)以及流體對葉片做正功的決定因素。

葉片壓力面和吸力面的壓力差可以通過公式(1)進行描述：

Δp(x,t)=-[(pu+ps)upper-(pu+ps)lower]

(1)

式中，下標u和s表示非穩(wěn)態(tài)量和穩(wěn)態(tài)量，upper表示吸力面，lower表示壓力面。

葉片表面某一截面上的受力可以通過式(2)進行計算：

(2)

式中:c為葉片弦長。由于葉片前緣和后緣處壓力變化大、噪聲多，本文只選取0.1～0.9弦長處壓力差進行積分。

不同葉高處的切向受力見圖10，結(jié)合相應(yīng)的位移數(shù)據(jù)，可以分析出不同葉高處的氣動穩(wěn)定性。

圖10　一個周期內(nèi)不同葉高處葉片切向力Fig.10 Tangential force of blade at different blade height in one period

圖11為一個周期內(nèi)不同葉高處葉片受力與位移之間的關(guān)系，箭頭方向表示時間。當沿箭頭方向逆時針旋轉(zhuǎn)時，流體在一個振動周期內(nèi)對葉片做正功(此時參考系隨葉片旋轉(zhuǎn)，后文討論的做功均在此參考系下獲得) ，氣動阻尼[20]為負，處于氣動不穩(wěn)定區(qū)；當沿箭頭方向順時針旋轉(zhuǎn)時，則相反。0.9、0.8、0.7葉高處均處于氣動不穩(wěn)定區(qū)，0.5和0.6葉高處曲線形成“8”字形，0.5葉高處的“8”字形分布對稱，此時處于氣動穩(wěn)定區(qū)，而0.6葉高處“8”字不對稱，處于穩(wěn)定與不穩(wěn)定之間的臨界狀態(tài)。同時可以根據(jù)曲線圍成面積的大小判斷流體對葉片做功的大小，在氣動不穩(wěn)定區(qū)(0.6葉高以上部分)，隨著葉高的增加，流體對葉片的做功也越大。

圖11　不同葉高處葉片位移與切向力Fig.11 Displacement and tangential force at different blade height

圖12　不同葉高處位移與切向力Fig.12 Displacement and tangential force at different blade height

圖12表明，不論是氣動穩(wěn)定區(qū)還是氣動不穩(wěn)定區(qū)都存在壓力波動，但是只有受力與位移滿足特定的相位差才能使得一個周期內(nèi)流體對葉片的做功為正。圖11比較了0.8葉高和0.5葉高處受力與位移之間的關(guān)系。在0.8葉高處，受力與位移相差0.4T的相位角，滿足0<0.4T<0.5T的關(guān)系；而在0.5葉高處，受力與位移的相位差在0.5T左右?？梢姎鈩臃€(wěn)定性不僅取決于激振力，而且還取決于受力與位移之間的相位差。形成不同相位差的原因有待進一步的研究。

3結(jié)論

本文以某國產(chǎn)600MW汽輪機末級葉片為研究對象，基于ANSYS的流固耦合和顫振的時域分析法，研究了汽輪機末級葉片顫振的判定、振型、影響因素和發(fā)生機理等問題。通過對葉片振動響應(yīng)、表面壓力分布和流體對葉片做功等數(shù)據(jù)的分析，得出以下結(jié)論：

(1)顫振的發(fā)作不僅需要周期性的壓力波動，還需要葉片受力與位移的相位差滿足特定的條件；

(2)葉片間相位差的存在改變了激波與葉片相互作用的位置，導(dǎo)致葉片表面壓力分布的波動；

(3)時域分析法可以相對更精確地預(yù)測顫振邊界，本文研究對象發(fā)生顫振的臨界工況在40%THA到45%THA之間；

(4)汽輪機末級葉片只有部分葉高處于氣動不穩(wěn)定區(qū)，對于本文所研究的葉片，60%葉高以上區(qū)域處于氣動不穩(wěn)定區(qū)。

參考文獻

[1]劉萬琨.汽輪機末級葉片顫振設(shè)計[J].東方電氣評論,2007,21(4):7-13.

LIU Wan-kun. Flutter design for steam turbine last blade[J]. Dongfang Electric Review,2007,21(4):7-13.

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

LI Yu-feng, REN Da-kang, HUANG Gang, et al. The off-design for LP last stage blade of air-cooling units[J]. Thermal Turbine,2004,33(1):14-163.

[3]王家楹.改善汽輪機低壓缸末級葉片調(diào)峰性能的探討[J].浙江電力,2001,(4):57-58.

WANG Jia-ying. Investigation on improving peak-load regulating performance of Final stage blade of steam turbine’s low pressure casing[J]. Zhejiang Electric Power,2001,(4):57-58.

[4]Whitehead D S. Vibration of cascade blades treated by actuator disc methods[J]. Proceedings of the Institution of Mechanical Engineers, 1959, 173(1): 555-574.

[5]楊曉東.預(yù)估葉片失速顫振的變形激盤法研究[D].北京航空航天大學(xué)，1984.

[6]陳佐一，劉紅，王繼宏,等. 汽輪機末級葉片失速顫振的全三維粘性流數(shù)值分析[J].中國電機工程學(xué)報,1999(3):18-20.

CHEN Zhuo-yi, LIU Hong, WANG Ji-hong. Three dimensional viscous flow numerical analysis of stall flutter in steam turbine blades[J]. Proceedings of the CSEE,1999(3):18-20.

[7]施永強.三維葉片顫振與葉片設(shè)計關(guān)聯(lián)性究[D].西安：西北工業(yè)大學(xué),2006.

[8]張興國,劉鋒,張陳安等.基于CFD技術(shù)的葉片顫振分析[J].航空計算技術(shù),2009,39(4):75-78.

ZHANG Xing-guo, LIU Feng, ZHANG Chenan. Flutter computation of turbomachinery blades based on CFD[J]. Aeronautical Computing Technique,2009,39(4):75-78.

[9]Bakhle M A, Reddy T S R, Keith T GK. Time domain flutter analysis of cascades using a full potential solver[J]. AIAA Paper 90-0984, 1990.

[10]弓三偉,任麗蕓,劉火星,等.汽輪機末級轉(zhuǎn)子葉片流固耦合計算[J].熱力透平,2007,36(3):153-157.

GONG San-wei, REN Li-yun, LIU Huo-xing, et al. Fluid-Solid interaction of rotor blade in last stage of steam turbine[J]. Thermal Turbine,2007,36(3):153-157.

[11]王征,吳虎,史亞鋒,等.基于CFD/CSD技術(shù)的壓氣機葉片流固耦合及顫振分析[J].航空動力學(xué)報,2011,26(5):1077-1084.

WANG Zheng, WU Hu, SHI Ya-feng, et al. Fluid-structureinteraction and flutter analysis of compressor blade based on CFD/CSD[J]. Journal of Aerospace Power, 2011,26(5):1077-1084.

[12]Massjung R. Discrete conservation and coupling strategies in nonlinear areoelasticity[J]. Computer Methods in Applied Mechanics and Engineering,2006,196(1/3):91-102.

[13]Sun D, Owen J S, Wright N G, et al. Fluid-structure interaction of prismatic line-like structures, using LES and block-iterative coupling[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2008, 96(6): 840-858.

[14]李磊,于明,李元生等.學(xué)科間信息傳遞的參數(shù)化空間散亂數(shù)據(jù)插值法[J].航空制造技術(shù),2011,(9):77-81.

LI Lei, YU Ming, LI Yuan-sheng, et al. Parameterized space scattered data interpolation method for inter- disciplinary information transferring[J]. Aeronautical Manufacturing Technology,2011,(9):77-81.

[15]Carta F O, Hilaire A O S. Experimentally determined stability parameters of a subsonic cascade oscillating near stall[J]. Journal of Engineering for Power, 1978,100:111-120.

[16]Carta F O. Unsteady aerodynamics and gapwise periodicity of oscillating cascaded airfoils[J]. Journal of engineering for power, 1983, 105(3): 565-574.

[17]Carta F O, Hilaire A O S. Effect of interblade phase angle and incidence angle on cascade pitching stability[J]. Journal of engineering for power, 1980, 102(2): 391-396.

[18]張小偉,王延榮.葉間相位角對葉片顫振的影響[J].航空動力學(xué)報,2010,25(2):412-416.

ZHANG Xiao-wei, WANG Yan-rong.Influence of interblade phase angle on the flutter of rotor blades[J]. Journal of Areospace Power,2010,25(2):412-416.

[19]楊炳淵,宋偉力.前緣激波脫體的大迎角翼面顫振工程計算方法[J].振動與沖擊,2002,21(4):100-103.

YAN Bing-yuan, SONG Wei-li. Engineering approximation of high attack flutter for wing with bow wave on leading edge[J]. Journal of Vibration and Shock,2002,21(4):100-103.

[20]Thirstrup Petersen J, Aagaard Madsen H, Bj?rck A, et al. Prediction of dynamic loads and induced vibrations in stall[M]. Roskilde:Riso National Laboratory,1998.