風(fēng)力機(jī)復(fù)合材料柔性葉片的顫振分析

任勇生，杜向紅，楊樹蓮

(1.山東科技大學(xué) 機(jī)械電子工程學(xué)院，山東青島 266510;2.山東工商學(xué)院，山東煙臺(tái) 264005)

為了最大限度地提高發(fā)電功率，降低發(fā)電成本，現(xiàn)代風(fēng)力發(fā)電機(jī)葉片的尺寸顯著增大。隨著風(fēng)力機(jī)葉片長度的增加以及復(fù)合材料在葉片制作中的廣泛使用，葉片的剛度也明顯降低，表現(xiàn)為更加柔性化的趨勢。在風(fēng)載荷的作用下，葉片的運(yùn)動(dòng)會(huì)影響周圍風(fēng)場從而改變作用于葉片的風(fēng)載荷，存在著空氣動(dòng)力、慣性力和彈性力的相互作用，因此，在葉片模態(tài)之間會(huì)發(fā)生振動(dòng)耦合，形成氣動(dòng)彈性不穩(wěn)定——顫振。葉片的顫振問題屬于氣動(dòng)彈性力學(xué)的研究范疇，是現(xiàn)代風(fēng)力機(jī)系統(tǒng)設(shè)計(jì)必須解決的主要隱患之一。近年來，在氣動(dòng)彈性力學(xué)的理論框架內(nèi)揭示葉片顫振行為發(fā)生的機(jī)理，并用于指導(dǎo)風(fēng)力機(jī)的氣動(dòng)彈性設(shè)計(jì)，已經(jīng)成為風(fēng)力機(jī)結(jié)構(gòu)動(dòng)力學(xué)研究中的一個(gè)新的重要內(nèi)容，日益受到風(fēng)能領(lǐng)域的重視［1］。

在風(fēng)力機(jī)葉片的顫振研究中，基于連續(xù)分布參數(shù)模型與典型截面集中參數(shù)模型相比，由于能夠反映葉片的剖面、變形以及慣性沿葉片展向的變化，因此可以更好地描述葉片的顫振行為。經(jīng)典的葉片分布參數(shù)模型是具有彎扭耦合變形梁的偏微分方程組，稱之為Hodges-Dowell方程［2］。文獻(xiàn)［3］基于 Hodges-Dowell方程和準(zhǔn)定常葉素氣動(dòng)力模型，取單模態(tài)近似表示葉片的位移，采用特征值法分析葉片揮舞/擺振/扭轉(zhuǎn)的穩(wěn)定性;文獻(xiàn)［4］引入重力、調(diào)距作用和轉(zhuǎn)速變化等因素的影響，將Hodges-Dowell方程進(jìn)行了擴(kuò)展;文獻(xiàn)［5］忽略葉片截面翹曲的影響，基于Euler梁理論并且采用模態(tài)疊加法，研究準(zhǔn)定常氣動(dòng)力作用下的風(fēng)力機(jī)葉片的受迫振動(dòng)問題。但是，在上述研究中，葉片均被視為由各向同性材料構(gòu)成的實(shí)心梁。事實(shí)上，典型的風(fēng)力機(jī)葉片具有復(fù)合材料空心結(jié)構(gòu)特點(diǎn)，通常需要采用薄壁閉室復(fù)合材料梁模型進(jìn)行描述，它的結(jié)構(gòu)動(dòng)力學(xué)建模，要求恰當(dāng)?shù)乇磉_(dá)由材料各向異性引起揮舞、擺振和扭轉(zhuǎn)變形之間的彈性耦合和橫截面翹曲等。文獻(xiàn)［6］雖然涉及薄壁復(fù)合材料葉片結(jié)構(gòu)的氣彈穩(wěn)定性問題，但沒有給出具體的求解過程;文獻(xiàn)［7］研究直升機(jī)旋翼復(fù)合材料薄壁葉片氣彈穩(wěn)定性問題，但所采用的是諧波平衡法(HBM)，計(jì)算過程比較復(fù)雜。

基于葉片連續(xù)分布參數(shù)模型的顫振問題研究，通常涉及高維偏微分方程組的求解，一般難以得到問題的精確解。雖然有限元方法作為一種數(shù)值近似方法，在復(fù)雜結(jié)構(gòu)的力學(xué)建模已經(jīng)得到廣泛的應(yīng)用，能夠提供各向異性復(fù)合材料葉片的結(jié)構(gòu)模型的詳細(xì)的解答，但是由于計(jì)算成本較高，對(duì)于葉片氣動(dòng)彈性問題的初步設(shè)計(jì)，缺乏實(shí)用性。而振型疊加法利用葉片的低階主要模態(tài)實(shí)現(xiàn)葉片顫振系統(tǒng)的模型降階，能夠極大地減少計(jì)算量，是風(fēng)力機(jī)氣動(dòng)彈性力學(xué)研究中經(jīng)常采用的一種有效的近似計(jì)算方法。

本文在現(xiàn)有研究結(jié)果的基礎(chǔ)上，提出一個(gè)風(fēng)力機(jī)葉片氣動(dòng)彈性力學(xué)理論模型，葉片結(jié)構(gòu)采用具有預(yù)扭轉(zhuǎn)的彎扭耦合閉合截面薄壁復(fù)合材料梁模型進(jìn)行描述，氣動(dòng)載荷采用葉素動(dòng)量理論和準(zhǔn)定常氣動(dòng)力理論進(jìn)行描述?；贖amilton原理并結(jié)合VAM［8］，進(jìn)行橫截面分析，導(dǎo)出旋轉(zhuǎn)復(fù)合材料葉片的結(jié)構(gòu)力學(xué)模型，在此基礎(chǔ)上，與葉素動(dòng)量理論和準(zhǔn)定常氣動(dòng)力理論進(jìn)行組合，建立葉片的顫振分析模型。將位移按廣義坐標(biāo)進(jìn)行模態(tài)展開，采用Galerkin法求解顫振分析模型。借助于特征值技術(shù)，對(duì)CAS構(gòu)型葉片的顫振性能進(jìn)行數(shù)值近似計(jì)算，揭示了入流比、預(yù)扭轉(zhuǎn)角和纖維鋪層角對(duì)風(fēng)力機(jī)葉片顫振性能的影響。

1 基本理論

1.1 位移場和應(yīng)變場

圖1表示具有任意截面形狀、長度為R的復(fù)合材料葉片，其變形前的軸線x相對(duì)旋轉(zhuǎn)平面傾斜小角度βp，稱之為預(yù)錐角，葉片以不變角速度Ω繞轉(zhuǎn)軸旋轉(zhuǎn)。為便于分析，建立三個(gè)坐標(biāo)系:葉片旋轉(zhuǎn)坐標(biāo)系(x，y，z)，坐標(biāo)原點(diǎn)o位于固定端;葉片局部坐標(biāo)系(x，η，ζ)，η和ζ是葉片橫截面的慣性主軸;葉片局部坐標(biāo)系(x，s，ξ)，s沿截面中線(即，薄壁葉片的中面與橫截面的交線)切向，ξ沿截面中線外法線方向。變形前的葉片截面，如圖2所示，假定截面質(zhì)心與彈性中心重合，彈性中心與氣動(dòng)中心的距離為xA，β是葉片橫截面的預(yù)扭轉(zhuǎn)角。

葉片坐標(biāo)系和主軸坐標(biāo)系之間，存在下列關(guān)系:

其中:β=β0(x/R)2;β0表示葉片尖端的預(yù)扭轉(zhuǎn)角。

葉片上的任意點(diǎn)沿坐標(biāo)系(x，y，z)的三個(gè)坐標(biāo)軸方向的位移分量為［9］:

其中:截面翹曲函數(shù)g(s，x)由位移場的連續(xù)性、環(huán)向力為零以及定常剪力的條件決定。

與位移場(2)相對(duì)應(yīng)的二階近似應(yīng)變場為:

其中:G、g1、g2、g3分別是與扭轉(zhuǎn)、軸向拉伸、繞z軸彎曲和y軸彎曲相關(guān)的翹曲分量;()'表示對(duì)x求偏導(dǎo)。

1.2 運(yùn)動(dòng)方程

為了導(dǎo)出葉片的氣動(dòng)彈性方程，利用Hamilton原理建立葉片運(yùn)動(dòng)方程:

其中，U和T分別為應(yīng)變能和動(dòng)能，可由下式確定:

外力的虛功:

其中:Fx，F(xiàn)y，F(xiàn)z分別是作用于x，y，z軸方向的分布?xì)鈩?dòng)力，Mx是繞彈性軸的氣動(dòng)力矩。

其中:

葉片截面內(nèi)力(矩)可以表示如下:

其中cij表示葉片截面剛度系數(shù)，它們與橫截面主軸的剛度系數(shù)kij之間滿足下列關(guān)系:

式中，kij的表達(dá)式見文獻(xiàn)［9］。

將式(11)的第一式代入方程組(7)的第一個(gè)方程，假定沿葉片軸向的氣動(dòng)力Fx=0，從x到R做定積分，得:

如果將方程式(8)、(9)和(12)代入方程組(7)，令Ω=0和β=0，且不計(jì)外載荷的作用，則可以導(dǎo)出與文獻(xiàn)［9］中的方程(27)一致的結(jié)果。

將式(11)和式(12)代入方程組(7)的后三個(gè)方程，得:

在本文的算例中，僅討論具有CAS構(gòu)型的薄壁復(fù)合材料葉片的情形，此時(shí)有:θ(y)=- θ(-y)，θ(z)=-θ(-z)，其中θ表示由正的s軸進(jìn)行度量的纖維鋪層角，于是得:k12=k13=k14=k24=k34=0，由此可以推出:C12=C13=C14=0，但是本文模型引入了葉片的預(yù)扭轉(zhuǎn)角β，故C24=C314≠0。因此可以看出，本文模型比文獻(xiàn)［9］的結(jié)果更具有普遍性。

1.3 氣動(dòng)載荷

根據(jù)文獻(xiàn)［3］的建議，基于準(zhǔn)定常葉素理論的葉片的分布?xì)鈩?dòng)力和力矩的表達(dá)式，可以導(dǎo)出如下:

其中:ρA為空氣密度;a為升力曲線斜率;b為無量綱半弦長=xA/R;αe=β+tan-1(Up/UT)+φ為有效攻角;CD0為形狀阻力系數(shù);Up、UT分別為葉片彈性軸沿?fù)]舞和擺振方向的相對(duì)速度分量;是葉片截面位置x和時(shí)間t的函數(shù)。不計(jì)陣風(fēng)并且略去變形高階量，則有:

其中:λ=(Vm-vi)/(RΩ)為入流比，Vm和vi分別為平均風(fēng)速和誘導(dǎo)速度。

1.4 求解方法

其中:下標(biāo)s和d分別為葉片的靜力和動(dòng)力分量。

將式(14)、式(15)和式(16)代入方程組(13)，導(dǎo)出氣彈方程組如下:

為了對(duì)氣彈方程組(17)進(jìn)行無量綱化，引入下列變換:

葉片的彎曲變形和扭轉(zhuǎn)變形可以表示如下:

根據(jù)Galerkin近似求解方法，利用假設(shè)振型(19)消除薄壁梁的空間位置變量，并結(jié)合變換關(guān)系(18)，可以將葉片氣動(dòng)彈性方程組(17)，轉(zhuǎn)化為關(guān)于時(shí)間的常微分方程如下:

其中:［M］，［C］，［K］分別表示關(guān)于廣義坐標(biāo){X}={Vj，Wj，φj}T的質(zhì)量矩陣、阻尼矩陣和剛度矩陣，為了節(jié)省篇幅，有關(guān)這些矩陣元素的具體計(jì)算公式，不再列出。

葉片的氣彈穩(wěn)定性取決于下列6N×6N矩陣A的特征值

在葉片的其它參數(shù)給定的情況下，矩陣A的特征值可以視為λ、β0和θ的函數(shù)。顯然，當(dāng)特征值的實(shí)部小于零，對(duì)應(yīng)的振動(dòng)模態(tài)是穩(wěn)定的;當(dāng)特征值的實(shí)部大于零，對(duì)應(yīng)的振動(dòng)模態(tài)是不穩(wěn)定的;而實(shí)部等于零，對(duì)應(yīng)的振動(dòng)模態(tài)是臨界穩(wěn)定的。

2 計(jì)算結(jié)果與討論

為了驗(yàn)證本文建立的模型及其近似計(jì)算方法的正確性，圖3給出CAS構(gòu)型復(fù)合材料箱型截面薄壁懸臂梁固有頻率的近似計(jì)算結(jié)果，并且與文獻(xiàn)［10］的固有頻率的精確解進(jìn)行了比較，結(jié)構(gòu)的幾何參數(shù)和材料參數(shù)均取自文獻(xiàn)［10］?？梢钥闯?，二者符合得很好。在計(jì)算過程中，我們發(fā)現(xiàn)，選取模態(tài)項(xiàng)數(shù)N=5，即可得到收斂的第一階ω1、二階揮舞為主的固有頻率ω3和第一階扭轉(zhuǎn)為主的固有頻率ω5。本文計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的對(duì)比，如表1所示。其中的實(shí)驗(yàn)數(shù)據(jù)取自文獻(xiàn)［11］中的圖4，圖6，圖7，圖10和圖11，實(shí)驗(yàn)用懸臂復(fù)合材料箱型梁的幾何尺寸、鋪層方式和材料特性，詳見文獻(xiàn)［11］的表1。

在下面的葉片顫振數(shù)值計(jì)算中，選取葉片的幾何尺寸為:截面寬度0.95 m，高度0.4 m，薄壁厚度h=0.047 5 m，單層厚度h/6，鋪層數(shù)為6層;葉片長度R=19 m。對(duì)于CAS構(gòu)型葉片，所采用鋪層方式為:上壁［θ］6，下壁［-θ］6，左壁和右壁［θ/-θ］3。葉片復(fù)合材料選取石墨/環(huán)氧樹脂，其性能參數(shù)為:E11=142 GPa，E22=E33=9.8 GPa，G12=G13=6.0 GPa，G23=4.83 GPa，ν12= ν13=0.42，ν23=0.50，ρ=1 601.1 kg/m3。

葉片其他參數(shù)取為:

表1 與文獻(xiàn)［11］固有頻率的實(shí)驗(yàn)結(jié)果的比較Tab.1 Comparison of frequencies(Hz)with experimental results from Ref.［11］

圖3 CAS構(gòu)型懸臂梁的第一、二階揮舞為主和第一階扭轉(zhuǎn)為主固有頻率隨鋪層角變化規(guī)律Fig.3 Natural frequencies vs.ply angle CAS cantilevered beam

圖4表示葉片顫振性能隨入流比λ的變化規(guī)律，并且還顯示了葉片預(yù)扭轉(zhuǎn)角的影響。其中，圖4(a)、圖4(c)表示葉片前兩階揮舞模態(tài)的情形(ζi，i=1，3)，圖4(b)、圖4(e)表示葉片前兩階擺振模態(tài)的情形(ζi，i=2，5)，圖4(d)、圖4(f)表示葉片前兩階扭轉(zhuǎn)模態(tài)的情形(ζi，i=4，6)。由圖 4 可見，當(dāng)預(yù)扭轉(zhuǎn)角 β0為零時(shí)，在所給定的入流比參數(shù)變化范圍內(nèi)，由于對(duì)應(yīng)的特征值的實(shí)部位于零附近，故葉片的一階和二階揮舞模態(tài)和擺振模態(tài)都是臨界穩(wěn)定的;而前兩階扭轉(zhuǎn)模態(tài)由于對(duì)應(yīng)的特征值實(shí)部小于零，因此是嚴(yán)格穩(wěn)定的。這說明，作用于葉片扭轉(zhuǎn)模態(tài)的氣彈阻尼相對(duì)較大，具有較強(qiáng)的抑制扭轉(zhuǎn)模態(tài)顫振的能力。當(dāng)β0增加時(shí)，葉片揮舞模態(tài)明顯地由臨界穩(wěn)定變?yōu)閲?yán)格穩(wěn)定的，β0的存在對(duì)揮舞模態(tài)的穩(wěn)定效果非常明顯，而擺振模態(tài)則由臨界穩(wěn)定變?yōu)椴环€(wěn)定的，并且其不穩(wěn)定性隨著β0的增加而增加，這很可能是由于作用于擺振方向上的氣動(dòng)阻力非常小，并且同時(shí)也沒有考慮結(jié)構(gòu)阻尼的緣故;與此同時(shí)，還可以看到，β0對(duì)扭轉(zhuǎn)模態(tài)的顫振性能的作用卻非常微弱，不會(huì)使扭轉(zhuǎn)模態(tài)的顫振性能產(chǎn)生本質(zhì)的變化。

圖5表示葉片顫振性能隨纖維鋪層角θ的變化規(guī)律，并且也顯示了葉片預(yù)扭轉(zhuǎn)角的影響。由圖5可見，總體而言，纖維鋪層角對(duì)揮舞模態(tài)和擺振模態(tài)的顫振的影響較大，而對(duì)扭轉(zhuǎn)模態(tài)顫振的影響較小。在預(yù)扭轉(zhuǎn)角不為零時(shí)，增加鋪層角對(duì)第一階揮舞模態(tài)和擺振模態(tài)的穩(wěn)定性，分別產(chǎn)生增強(qiáng)和削弱的作用。除了在鋪層角為10°附近，對(duì)第二階揮舞模態(tài)的穩(wěn)定性有所削弱之外，在整個(gè)鋪層角變化范圍內(nèi)，不會(huì)帶來明顯的影響，而鋪層角的變化對(duì)第二階擺振模態(tài)的穩(wěn)定性，幾乎不會(huì)產(chǎn)生影響。

圖4 入流比對(duì)葉片的穩(wěn)定性的影響(取不同的預(yù)扭轉(zhuǎn)角)Fig.4 Effect of inflow on blade stability for different pre-twist angle

圖6表示葉片顫振性能隨預(yù)扭轉(zhuǎn)角β0的變化規(guī)律，同時(shí)也顯示了葉片纖維鋪層角的影響。由圖6可以清楚地看到，纖維鋪層角對(duì)葉片的高階揮舞與擺振模態(tài)的顫振性能的影響很小，而預(yù)扭轉(zhuǎn)角對(duì)葉片扭轉(zhuǎn)模態(tài)的顫振性能的影響很小。這個(gè)結(jié)論與圖5的結(jié)果是相互對(duì)應(yīng)的。

3 結(jié)論

本文基于彎扭耦合閉合截面薄壁復(fù)合材料梁的分析理論，采用Hamilton原理并且借助于Galerkin法，建立了在準(zhǔn)定常氣動(dòng)力作用下的風(fēng)力機(jī)葉片顫振分析模型。研究結(jié)果表明，采用本文建立的模型及其近似計(jì)算方法所獲得的葉片的固有振動(dòng)頻率數(shù)值近似解，與已有文獻(xiàn)的精確解相當(dāng)一致。針對(duì)一個(gè)長度為18 m的石墨/環(huán)氧樹脂復(fù)合材料葉片，在轉(zhuǎn)速 Ω=4π/3(rads-1)的情況下的顫振分析結(jié)果表明:

(1)利用所建立的模型能夠?qū)︼L(fēng)力機(jī)葉片的顫振性能進(jìn)行參數(shù)分析，以確定影響葉片顫振性能的主要因素，為進(jìn)一步改善或提高葉片的氣動(dòng)彈性性能，提供有用的信息。

(2)葉片的預(yù)扭轉(zhuǎn)角對(duì)葉片的顫振性能會(huì)產(chǎn)生顯著的影響，但對(duì)于揮舞模態(tài)和擺振模態(tài)來說，預(yù)扭轉(zhuǎn)角作用效果是相反的。在實(shí)際應(yīng)用中，預(yù)扭轉(zhuǎn)角是影響風(fēng)力機(jī)葉片顫振性能的一個(gè)主要設(shè)計(jì)參數(shù)。以往的研究表明，葉片的預(yù)扭轉(zhuǎn)角可以借助于形狀記憶合金和壓電材料的驅(qū)動(dòng)而產(chǎn)生，因此，上述研究結(jié)果對(duì)于客觀地評(píng)價(jià)智能材料在葉片顫振抑制中的應(yīng)用效果，有一定的參考價(jià)值。

［1］ Hansen M O L ，Sorensen J N，Voutsinas S，et al.State of the art in wind turbine aerodynamics and aeroelasticity［J］.Progress in Aerospace Sciences，2006，42(4):285-330.

［2］Hodges D H，Dowell E H.Nonlinear equations of motion for the elastic bending and torsion of twisted nonuniform rotor blades［R］.Technical Report TN D-7818，NASA，December 1974.

［3］Kottapalli S B R，F(xiàn)riedmann P P.Aeroelastic stability and response of horizontal axis wind turbine blades［J］.AIAA Journal，1979，17(12):1381-1389.

［4］ Kallesoe B S.Equations of motion for a rotor blade，including gravity，pitch action and rotor speed variations［J］.Wind Energy，2007，10:209-230.

［5］Kim K T，Lee C W.Forced vibration analysis of flexible wind turbine blades by using assumed modes method［C］.Proc.of 6thInt.Vibration Engineering Conference，2008，192-201.

［6］Wilkie W K，Belvin W K.Aeroelastic analysis of helicopter rotorblades incorporating anisotropic piezoelectric twist actuation［C］.Technical Report:NASA-96-asme-wkw Year of Publication:1996.

［7］Kim T Y.Nonlinear large amplitude structural and aeroelastic behavior of composite rotor blades at large static deflection［D］.Dorctor Dissertation， Massachusetts Institute of Technology，1992.

［8］Berdichevsky V，Armanios E A，Badir A M.Theory of anisotropic thin-walled closed-cross-section beams［J］.Composite Engineering(Special Issue)，1992，2(5-7):411-432.

［9］Armanios E A，Badir A M.Free vibration analysis of anisotropic thin-walledclosed-section beams［J］.AIAA J，1995，33(10):1905-1910.

［10］ Dancila D S，Armanios E A.The influence of coupling on the free vibration of anisotropic thin-walled closed-section beams［J］.Int J Solids Struct，1998，35(23):3105-3119.

［11］ Chandra R，Chopra I.Experimental theoretical investigation of the vibration characteristics of rotating composite box beams［J］.J Aircraft，1992，29(4):657-664.