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

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

(武漢大學(xué)動(dòng)力與機(jī)械學(xué)院， 武漢　430072)

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

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

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

(武漢大學(xué)動(dòng)力與機(jī)械學(xué)院， 武漢430072)

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

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

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目(51376140)

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

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

中圖分類(lèi)號(hào):Tp62

文獻(xiàn)標(biāo)志碼：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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.1顫振預(yù)測(cè)

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

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

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

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

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

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

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

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

2.2壓力分布

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

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

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

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

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

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

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

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

(1)

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

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

(2)

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

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

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

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

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

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

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

3結(jié)論

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

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

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

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

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

參考文獻(xiàn)

[1]劉萬(wàn)琨.汽輪機(jī)末級(jí)葉片顫振設(shè)計(jì)[J].東方電氣評(píng)論,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]李宇峰,任大康,黃鋼,等.空冷汽輪機(jī)低壓末級(jí)變工況設(shè)計(jì)[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]王家楹.改善汽輪機(jī)低壓缸末級(jí)葉片調(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ù)估葉片失速顫振的變形激盤(pán)法研究[D].北京航空航天大學(xué)，1984.

[6]陳佐一，劉紅，王繼宏,等. 汽輪機(jī)末級(jí)葉片失速顫振的全三維粘性流數(shù)值分析[J].中國(guó)電機(jī)工程學(xué)報(bào),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]施永強(qiáng).三維葉片顫振與葉片設(shè)計(jì)關(guān)聯(lián)性究[D].西安：西北工業(yè)大學(xué),2006.

[8]張興國(guó),劉鋒,張陳安等.基于CFD技術(shù)的葉片顫振分析[J].航空計(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]弓三偉,任麗蕓,劉火星,等.汽輪機(jī)末級(jí)轉(zhuǎn)子葉片流固耦合計(jì)算[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ī)葉片流固耦合及顫振分析[J].航空動(dòng)力學(xué)報(bào),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]張小偉,王延榮.葉間相位角對(duì)葉片顫振的影響[J].航空動(dòng)力學(xué)報(bào),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ì)算方法[J].振動(dòng)與沖擊,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.