張 穎， 安利強(qiáng)， 王璋奇

(1.華北電力大學(xué) 能源動力與機(jī)械工程學(xué)院，河北 保定 071003；2.河北大學(xué) 建筑工程學(xué)院，河北 保定 071002)

0 引 言

風(fēng)力發(fā)電機(jī)葉片隨著尺寸的增加，輸出功率越來越大。風(fēng)機(jī)的風(fēng)輪直徑從幾m到170多m，實(shí)現(xiàn)發(fā)電量從kW級別增加到MW級別?，F(xiàn)代商用5 MW風(fēng)機(jī)葉片長度超過60 m，目前研發(fā)的額定功率為8～10 MW的風(fēng)力發(fā)電機(jī)葉片長度超過80 m。增加葉片長度可以捕獲更多的風(fēng)能，同時也帶來了新的挑戰(zhàn)。本文研究這類細(xì)長的柔性葉片，隨著長度的增加，可能發(fā)生氣動彈性不穩(wěn)定性問題[1-3]。未來的風(fēng)力發(fā)電機(jī)重要的工程問題是如何考慮葉片的柔性對氣動彈性不穩(wěn)定性的影響[4]。

大型風(fēng)力發(fā)電機(jī)葉片是具有高縱橫比的柔性葉片。柔性葉片及其周圍流動之間的高度非線性相互作用導(dǎo)致葉片出現(xiàn)各種動態(tài)不穩(wěn)定性[5，6]。

20 kW風(fēng)力發(fā)電機(jī)剛性葉片的顫振速度比其正常運(yùn)行速度提高幾倍，而1.5 MW柔性葉片的顫振速度大約是正常運(yùn)行速度的1.75倍[7]。由于葉片長度的增加，顫振速度更接近葉片的正常運(yùn)行速度，因此在運(yùn)行過程中更容易出現(xiàn)顫振。Hansen[8]回顧了現(xiàn)代商業(yè)風(fēng)力發(fā)電機(jī)可能發(fā)生的兩類氣動彈性問題，失速型風(fēng)機(jī)的失速顫振和變槳距風(fēng)機(jī)的經(jīng)典顫振。指出經(jīng)典顫振可能是在附著流狀態(tài)下運(yùn)行的變槳距調(diào)節(jié)變速風(fēng)機(jī)葉片的動態(tài)不穩(wěn)定性。顫振可能發(fā)生在承受氣動載荷的柔性結(jié)構(gòu)中，并且由兩種結(jié)構(gòu)模式的結(jié)合引起的。已有文獻(xiàn)風(fēng)機(jī)葉片運(yùn)行中，經(jīng)典的氣彈顫振不是風(fēng)力發(fā)電機(jī)設(shè)計的首要因素，但是隨著葉片越來越大和越來越柔，顫振在實(shí)際風(fēng)機(jī)葉片中已經(jīng)逐漸成為設(shè)計的考慮因素。為了研究顫振發(fā)生時刻的特征，對于線彈性模型，隨著葉片長度增加，預(yù)測葉片顫振速度與風(fēng)輪運(yùn)行轉(zhuǎn)速比越來越低，葉片在運(yùn)行期間發(fā)生顫振的可能性增加[8，9]。Pourazarm[10]等人研究了扭轉(zhuǎn)和揮舞的固有頻率對臨界顫振速度的影響。研究表明，在扭轉(zhuǎn)固有頻率不變的情況下，僅僅增加揮舞自振頻率能夠降低葉片的顫振頻率。然而，扭轉(zhuǎn)頻率對于氣動彈性穩(wěn)定性的影響更為顯著。隨著扭轉(zhuǎn)固有頻率的降低，顫振甚至在低于額定速度的情況下發(fā)生。

一些研究人員使用各種分析工具量化MW級風(fēng)機(jī)葉片的穩(wěn)定性極限，Lobitz[7]基于NASTRAN對靜止空氣中的水平軸風(fēng)力發(fā)電機(jī)葉片進(jìn)行顫振分析，該模塊最初是為分析垂直軸風(fēng)力發(fā)電機(jī)葉片開發(fā)的。他研究預(yù)測WindPact 1.5 MW葉片臨界顫振速度為42.3 r/m(是額定轉(zhuǎn)速的兩倍)，顫振頻率為6.234 Hz。Lobitz[11]進(jìn)一步研究縮比結(jié)構(gòu)對風(fēng)力發(fā)電機(jī)葉片穩(wěn)定性的影響，對比研究長度為35 m的Windpact 1.5 MW 葉片和縮比后長度為9.86 m的GX-100葉片，縮比系數(shù)為0.271 8，設(shè)計相同的尖速比和展弦比，得出兩支葉片單位轉(zhuǎn)速下的顫振速度接近，進(jìn)一步指出了Windpact 1.5 MW葉片在降低扭轉(zhuǎn)頻率與第二階揮舞頻率比和將葉片質(zhì)心在彈性軸后移兩種方式下，在較低的轉(zhuǎn)速下發(fā)生顫振。Hansen[8]使用氣動伺服彈性工具HAWCStab穩(wěn)定性工具，研究了美國可再生能源實(shí)驗(yàn)室NREL5MW 葉片[12]的氣動不穩(wěn)定性，預(yù)測葉片顫振速度為24 r/m,是額定轉(zhuǎn)速的兩倍，進(jìn)行了風(fēng)輪顫振速度分析結(jié)果相似。Stablein[13]等研究了DTU 10 MW 風(fēng)力發(fā)電機(jī)的模態(tài)特性和穩(wěn)定性，他們使用氣動伺服彈性工具HAWCSTAb2 ,通過穩(wěn)態(tài)平衡下的特征值分析，研究了揮舞-擺振和扭轉(zhuǎn)耦合的葉片。結(jié)果表明，當(dāng)擺振-扭轉(zhuǎn)順槳耦合時，第一階擺振模態(tài)的阻尼增加，而在擺振-扭轉(zhuǎn)失速耦合時阻尼降低。揮舞-扭轉(zhuǎn)順槳耦合時，揮舞模態(tài)的頻率增加，阻尼降低，然而對于揮舞-扭轉(zhuǎn)失速耦合時，揮舞模態(tài)的頻率降低，阻尼增加。Owens 和Griffith[14]采用Blast設(shè)計工具預(yù)測Windpact 1.5 MW和SNL100葉片的顫振，分別在40.6 r/m(額定轉(zhuǎn)速的兩倍)和13.05 r/m(額定轉(zhuǎn)速的1.75倍)觀察到顫振現(xiàn)象。在我國某風(fēng)電場，6 MW機(jī)組風(fēng)輪最大直徑達(dá)到百米級別，發(fā)現(xiàn)某長度大于70 m的葉片由于振動發(fā)生損傷，發(fā)現(xiàn)損傷葉片時間處于掛機(jī)后尚未并網(wǎng)投運(yùn)前，或并網(wǎng)投運(yùn)最初三個月內(nèi)。位置集中于葉身中段后緣區(qū)域，主梁和后緣完好，但蒙皮區(qū)域發(fā)生了類似失穩(wěn)表象的破壞，葉片最終破壞形式是解體式的破壞，解體后葉片各部位保持相對的完整。以往文獻(xiàn)中僅僅研究一兩個參數(shù)，未對葉片的結(jié)構(gòu)參數(shù)深入研究。

本文主要研究NREL 5 MW葉片揮舞剛度、扭轉(zhuǎn)剛度以及葉片的質(zhì)量矩、重心偏移等幾個參數(shù)的變化對葉片顫振的影響規(guī)律。利用哈密頓原理推導(dǎo)彎扭耦合連續(xù)梁偏微分方程，結(jié)合Theodorson非定常氣動荷載，求解復(fù)特征值方程，得到葉片顫振發(fā)生時刻的顫振頻率和顫振速度。按照不同翼型族將風(fēng)機(jī)葉片結(jié)構(gòu)參數(shù)分成三個區(qū)域，通過數(shù)值分析，研究區(qū)域化結(jié)構(gòu)參數(shù)變化對葉片耦合模態(tài)顫振的影響規(guī)律。

1 風(fēng)力發(fā)電機(jī)葉片的氣動彈性模型

氣動彈性模型由兩部分組成：結(jié)構(gòu)模型和氣動模型。在結(jié)構(gòu)模型中，歐拉-伯努利梁，因細(xì)長的結(jié)構(gòu)廣泛用于風(fēng)力發(fā)電機(jī)建模，利用哈密頓原理導(dǎo)出相應(yīng)的控制微分方程，采用西奧道森理論求解非定常氣動荷載的計算。

1.1 結(jié)構(gòu)模型

基于Hodges[15]給出的懸臂旋轉(zhuǎn)歐拉-伯努利梁的運(yùn)動微分方程，假設(shè)梁是各向同性，在恒定角速度下的揮舞彎曲，擺振彎曲，拉伸和扭轉(zhuǎn)運(yùn)動微分方程，坐標(biāo)系(XYZ)表示剛性連接到葉片根部的主體坐標(biāo)系，X軸對應(yīng)于未變形葉片的彈性軸，Y軸位于旋轉(zhuǎn)平面內(nèi)，Z軸垂直于旋轉(zhuǎn)平面，如圖1所示。點(diǎn)O為未變形彈性軸，O’為變形彈性軸，位移分別為u,v,w，其中，u,v,w分別代表軸向，揮舞和擺振方向位移。它繞變形的彈性軸旋轉(zhuǎn)，角度β表示葉片的幾何槳距角，是扭轉(zhuǎn)和預(yù)扭轉(zhuǎn)角的總和。

圖1 葉片坐標(biāo)系及其彈性位移Fig. 1 Coordinates of blade and its elastic displacements

與變形葉片的橫截面正交的局部坐標(biāo)系是ξ,η,ζ，其中ξ軸與偏轉(zhuǎn)的彈性軸相切，η和ζ是截面的主軸。

如圖1所示，其中，X,Y,Z是參考軸，旋轉(zhuǎn)槳葉的建模考慮了揮舞彎曲和扭轉(zhuǎn)兩個方向的運(yùn)動。運(yùn)動方程是由哈密頓原理導(dǎo)出

(1)

式中：δU為應(yīng)變能的變分；δT為動能的變分δW是做的虛功。Hodges和Dowell針對考慮各向同性風(fēng)輪推導(dǎo)了δU和δT。Hong和Chopra修正了考慮復(fù)合材料正交各向異性的方程。僅僅考慮揮舞彎曲和扭轉(zhuǎn)的復(fù)合材料風(fēng)機(jī)葉片控制微分方程如下：

(2)

(3)

邊界條件為

x=0,w=0,w′=0,φ=0x=L,w″=0,

w?=0,φ′=0

(4)

式中：ω和φ分別代表揮舞和扭轉(zhuǎn)位移；GJ是扭轉(zhuǎn)剛度；EI是揮舞彎曲剛度；Km是繞彈性軸的極回轉(zhuǎn)半徑；Km1和Km2分別是繞主中性軸和垂直于彈性軸的軸質(zhì)量回轉(zhuǎn)半徑；m為單位長度的葉片質(zhì)量；ρ為葉片的密度；A為橫截面積；e為質(zhì)心偏移量，彈性軸和質(zhì)心的距離；Ω為風(fēng)輪轉(zhuǎn)速；L和M為氣動載荷作用在葉片上產(chǎn)生的氣動力和變槳距。

1.2 氣動載荷

結(jié)合Theodosen理論[16]，非穩(wěn)態(tài)氣動載荷在2D翼型上作用的變槳距M和升力L如下式：

(5)

(6)

2 風(fēng)力發(fā)電機(jī)葉片顫振速度計算方法

將方程(5)和(6)代入到(2)和(3)中，結(jié)合有限元方法，將葉片看作變截面線性歐拉-伯努利梁，每個節(jié)點(diǎn)自由度為3個，分別為w,δw/δx,φ。方程(2)，(3)的變分形式要求單元的插值函數(shù)連續(xù)，并且φ一階導(dǎo)數(shù)非零，w二階導(dǎo)數(shù)非零。結(jié)構(gòu)區(qū)域離散為1維2節(jié)點(diǎn)單元，每個節(jié)點(diǎn)自由度的φ(x) 和w(x)分別假設(shè)為線性函數(shù)和三次多項(xiàng)式。

(7)

代入(2)和(3)，利用 Galerkin 方法[17]，在單元上加權(quán)積分，得到如下微分方程。

(8)

式中：Ms和Ma結(jié)構(gòu)和氣動質(zhì)量矩陣，Ca為氣動阻尼矩陣，Ks和Ka是結(jié)構(gòu)和氣動剛度矩陣。方程(5)和(6)是顫振頻率ff的函數(shù)，方程(8)通過迭代方法求解復(fù)特征值問題，得到顫振頻率和顫振阻尼。阻尼比計算利用復(fù)特征值的負(fù)實(shí)部除以復(fù)特征值的模得到，阻尼比變?yōu)榱惚硎景l(fā)生顫振。顫振發(fā)生的頻率和速度定義為顫振頻率和顫振速度。通常認(rèn)為阻尼比為負(fù)值時，容易發(fā)生顫振，因?yàn)槿~片要不斷從空氣中吸收能量，積累到一定程度會產(chǎn)生氣動彈性不穩(wěn)定[18]。顫振頻率發(fā)生在第一階扭轉(zhuǎn)和第三階彎曲模態(tài)耦合的時刻，此時，顫振速度對應(yīng)氣動阻尼比為零。

3 5 MW風(fēng)力發(fā)電機(jī)葉片顫振分析

3.1 5 MW風(fēng)力發(fā)電機(jī)葉片顫振速度

利用文獻(xiàn)[12]的葉片翼型參數(shù)和截面特性，編制MATLAB有限元顫振分析程序，計算結(jié)果如表1和表2所示。表1中的自振頻率為固有頻率，和已有文獻(xiàn)對比，表明程序設(shè)計正確，表2中為顫振頻率與顫振速度對比，與Pourazarm[10]計算結(jié)果接近。

如圖2所示，葉片隨著轉(zhuǎn)速的增加，顫振頻率和氣動阻尼的變化規(guī)律，在轉(zhuǎn)速21.4 r/m時，一階扭轉(zhuǎn)模態(tài)和三階揮舞彎曲模態(tài)耦合，負(fù)阻尼產(chǎn)生時對應(yīng)的顫振頻率為4.7 Hz。

表1 葉片自振頻率Tab.1 Comparison of natural frequencies

表2 葉片顫振速度和顫振頻率Tab.2 Comparison of flutter speed and frequencies

圖2 5 MW葉片顫振頻率和阻尼比Fig. 2 Damping ratio and aeroelastic frequency of NERL 5 MW wind turbine blade

3.2 葉片結(jié)構(gòu)參數(shù)對氣動彈性穩(wěn)定性影響

風(fēng)力發(fā)電機(jī)葉片的氣動彈性不穩(wěn)定現(xiàn)象屬于模態(tài)耦合，當(dāng)葉片的第一階扭轉(zhuǎn)模態(tài)和某一階揮舞模態(tài)發(fā)生耦合，則出現(xiàn)經(jīng)典顫振。因此，如果結(jié)構(gòu)的參數(shù)發(fā)生變化，即使相同的風(fēng)速和攻角作用下，葉片的顫振頻率將產(chǎn)生變化。大型風(fēng)機(jī)葉片沿著展向，具有不同的翼型截面，每個翼型在梁帽、腹板、前緣和尾緣都具有不同的復(fù)合材料鋪層。按照翼型特性，圖3所示，將葉片分為三個不同區(qū)域，分別為葉根，葉中，葉尖區(qū)域。假設(shè)葉片的幾何形狀不變，作用在葉片上的氣動力不變，本文對比研究不同區(qū)域葉片的結(jié)構(gòu)參數(shù)變化，初步分析耦合模態(tài)顫振規(guī)律。

圖3 葉片的不同區(qū)域Fig. 3 Different airfoil regions for 5MW blade

考慮到葉片的翼型截面特殊特性，主要從葉片抗彎剛度、抗扭剛度、質(zhì)量矩，質(zhì)心偏移量的變化等幾個方面，分析結(jié)構(gòu)參數(shù)變化引起的顫振頻率及顫振速度變化規(guī)律。

(1) 揮舞剛度變化

(9)

葉片簡化為懸臂結(jié)構(gòu)，按照懸臂梁理論，則有關(guān)系式(10)成立。揮舞固有頻率和抗彎剛度EI平方根成正比，與單位長度的質(zhì)量平方根、長度的平方成反比。假定其它因素不變，研究揮舞剛度變化引起的頻率和轉(zhuǎn)速。

(10)

根據(jù)式(10)，降低葉片剛度EI或增大葉片質(zhì)量和長度，降低葉片的揮舞自振頻率，表明葉片的揮舞方向更柔，相應(yīng)圖5的轉(zhuǎn)速變化規(guī)律，圖4中，當(dāng)EI*<0.4,葉片的尖部對顫振頻率影響較大，剛度減小至30%，葉尖顫振頻率變化率最大約為71%。0.41三者接近，葉尖和中部區(qū)域?qū)︻澱耦l率的影響較大。

圖4 顫振頻率隨揮舞剛度的變化圖Fig. 4 Flutter frequency versus dimensionless flapwise stiffness

對于顫振速度，隨著抗彎剛度的減小，對應(yīng)的揮舞模態(tài)頻率減小，和扭轉(zhuǎn)模態(tài)耦合的時刻延后，因此顫振速度相對基準(zhǔn)速度增大，當(dāng)EI*<0.4，葉尖區(qū)域顫振速度變化為基準(zhǔn)速度的1.15倍，葉根和葉中區(qū)域顫振速度保持在1.34倍。EI*>0.4，葉根區(qū)域變化引起的顫振速度迅速減小至基準(zhǔn)速度。葉尖和葉中區(qū)域的變化引起的顫振速度從基準(zhǔn)速度的1.3倍逐漸減小至基準(zhǔn)顫振速度，當(dāng)剛度繼續(xù)增加，顫振速度減小緩慢，趨于基準(zhǔn)速度的0.9倍。

圖5 顫振轉(zhuǎn)速隨揮舞剛度的變化圖Fig. 5 Flutter speed versus dimensionless flapwise stiffness

(2)扭轉(zhuǎn)剛度變化

定義無量綱扭轉(zhuǎn)剛度為

(11)

扭轉(zhuǎn)固有頻率與扭轉(zhuǎn)剛度的平方根成正比，與單位長度的質(zhì)量平方根、長度和回轉(zhuǎn)半徑成反比，由式(12)可得

(12)

扭轉(zhuǎn)剛度減小，扭轉(zhuǎn)頻率減小，不再是扭轉(zhuǎn)和三階揮舞模態(tài)耦合，而是扭轉(zhuǎn)和二階揮舞或一階揮舞模態(tài)耦合產(chǎn)生的顫振.如圖6，隨著扭轉(zhuǎn)剛度的減小，顫振頻率逐漸減小，葉中和葉尖區(qū)域?qū)︻澱耦l率的影響程度接近。如圖7所示，葉尖和葉中區(qū)域的顫振轉(zhuǎn)速隨著扭轉(zhuǎn)剛度減小迅速減小，葉根區(qū)域減小緩慢。葉尖區(qū)域的扭轉(zhuǎn)剛度對顫振速度的影響程度最大。

圖6 顫振頻率隨扭轉(zhuǎn)剛度的變化圖Fig. 6 Flutter frequency versus dimensionless torsional stiffness

圖7 顫振轉(zhuǎn)速隨扭轉(zhuǎn)剛度的變化圖Fig. 7 Flutter speed versus dimensionless torsional stiffness

(3)回轉(zhuǎn)半徑變化

葉片質(zhì)量矩和回轉(zhuǎn)半徑平方成正比。因此，定義葉片對主中性軸和過彈性軸且垂直于弦的坐標(biāo)軸的無量綱回轉(zhuǎn)半徑

(13)

公式(12)可知，扭轉(zhuǎn)頻率與極慣性矩Km成反比，回轉(zhuǎn)半徑Km1和Km2的大小反映了橫截面的結(jié)構(gòu)特性，體現(xiàn)彈性軸的質(zhì)量分布特性，且有式(14)成立

(14)

分別改變對葉片截面的兩個質(zhì)量矩，發(fā)現(xiàn)對于葉片的截面的y軸的回轉(zhuǎn)半徑Km1*對顫振轉(zhuǎn)速和顫振頻率的影響較大。在Km1*=0.5時，達(dá)到最大值為基準(zhǔn)頻率的1.23倍。在Km1*>0.5時，隨著回轉(zhuǎn)半徑的增大，顫振頻率減小。說明對于截面的關(guān)于質(zhì)量分布離彈性軸越遠(yuǎn)，更容易發(fā)生顫振。在區(qū)間Km1*<0.4，葉片中部和葉片尖部的顫振頻率也隨著回轉(zhuǎn)半徑增大而減小，總體頻率均小于基準(zhǔn)顫振頻率，Km1*=0.4，葉中區(qū)域頻率最小為基準(zhǔn)頻率的0.76倍。

圖8 顫振頻率隨回轉(zhuǎn)半徑的變化圖Fig. 8 Flutter frequency versus radius of gyration

圖9 顫振轉(zhuǎn)速隨回轉(zhuǎn)半徑的變化圖Fig. 9 Flutter speed versus the radius of gyration

(4)質(zhì)心位置變化

定義無量綱質(zhì)心偏移量

(15)

圖10 顫振頻率隨質(zhì)心的變化圖Fig. 10 Flutter frequency versus mass offset

在葉片截面中，重心和彈性軸的距離e*，是一個影響葉片結(jié)構(gòu)整體的動態(tài)不穩(wěn)定性的參數(shù)。質(zhì)心向葉片后緣移動，才可能發(fā)生顫振。由圖10可見，隨著質(zhì)心偏移量e*增加，顫振頻率略有增加，葉尖部位的截面影響稍微大些，最大顫振頻率變化率僅僅2%。對于圖11的顫振速度，葉尖部位影響較大，當(dāng)減小e*，質(zhì)心向尾緣后移，顫振速度減小，更容易發(fā)生顫振。

圖11 顫振速度隨質(zhì)心的變化圖Fig. 11 Flutter speed versus mass offset

4 結(jié) 論

本文通過對不同區(qū)域葉片參數(shù)變化，利用復(fù)特征值方法求解葉片顫振特征方程，對NREL 5 MW 風(fēng)力發(fā)電機(jī)葉片的耦合模態(tài)顫振參數(shù)進(jìn)行研究，驗(yàn)證了葉片第一階扭轉(zhuǎn)模態(tài)和某一階揮舞彎曲模態(tài)耦合能引起顫振的現(xiàn)象。研究結(jié)果表明：

(1)通過分析葉尖、葉根、葉中三塊區(qū)域的結(jié)構(gòu)參數(shù)變化，葉尖區(qū)域的揮舞剛度對顫振頻率影響較大，當(dāng)揮舞剛度減小至30%，葉尖顫振頻率變化率最大約為71%，顫振速度減小為基準(zhǔn)速度的1.15倍，而葉根和葉中區(qū)域顫振速度保持在1.34倍。當(dāng)揮舞剛度變化率從40%開始繼續(xù)增加，葉根區(qū)域引起的顫振速度迅速減小至基準(zhǔn)速度,葉尖和葉中區(qū)域的變化引起的顫振速度從基準(zhǔn)速度的1.3倍逐漸減小至基準(zhǔn)顫振速度，三個區(qū)域最終趨于基準(zhǔn)速度的0.9倍。

(2)扭轉(zhuǎn)剛度與顫振速度和顫振頻率關(guān)系線性相關(guān)。葉尖區(qū)域的扭轉(zhuǎn)剛度對轉(zhuǎn)速影響較大，減小至30%出現(xiàn)轉(zhuǎn)速突變?；剞D(zhuǎn)半徑Km1*對顫振頻率影響較大，葉尖區(qū)域的變化率為50%時，為基準(zhǔn)頻率的1.23倍，葉中區(qū)域變化率為40%時，為基準(zhǔn)頻率的0.76倍。回轉(zhuǎn)半徑Km2*變化與顫振速度和顫振頻率近似成反比。Km1*在葉尖處存在一段增長緩慢區(qū)域達(dá)到基準(zhǔn)速度的1.32倍迅速減小至基準(zhǔn)速度。葉尖部位的質(zhì)心偏移對葉片的顫振速度影響較大。

(3)對于MW級風(fēng)機(jī)，葉片顫振問題，扭轉(zhuǎn)頻率和揮舞頻率將成為未來設(shè)計的必要考慮因素。如何保持較高的扭轉(zhuǎn)頻率，避免和某一階彎曲模態(tài)發(fā)生耦合，成為避免發(fā)生顫振的關(guān)鍵因素。