任勇生，杜向紅，楊樹蓮

(1.山東科技大學 機械電子工程學院，山東青島 266510;2.山東工商學院，山東煙臺 264005)

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

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

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

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

1 基本理論

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

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

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

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

葉片上的任意點沿坐標系(x，y，z)的三個坐標軸方向的位移分量為［9］:

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

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

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

1.2 運動方程

為了導出葉片的氣動彈性方程，利用Hamilton原理建立葉片運動方程:

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

外力的虛功:

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

其中:

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

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

式中，kij的表達式見文獻［9］。

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

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

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

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

1.3 氣動載荷

根據(jù)文獻［3］的建議，基于準定常葉素理論的葉片的分布氣動力和力矩的表達式，可以導出如下:

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

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

1.4 求解方法

其中:下標s和d分別為葉片的靜力和動力分量。

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

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

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

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

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

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

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

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

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

在下面的葉片顫振數(shù)值計算中，選取葉片的幾何尺寸為:截面寬度0.95 m，高度0.4 m，薄壁厚度h=0.047 5 m，單層厚度h/6，鋪層數(shù)為6層;葉片長度R=19 m。對于CAS構(gòu)型葉片，所采用鋪層方式為:上壁［θ］6，下壁［-θ］6，左壁和右壁［θ/-θ］3。葉片復合材料選取石墨/環(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 與文獻［11］固有頻率的實驗結(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ī)律，并且還顯示了葉片預扭轉(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 可見，當預扭轉(zhuǎn)角 β0為零時，在所給定的入流比參數(shù)變化范圍內(nèi)，由于對應(yīng)的特征值的實部位于零附近，故葉片的一階和二階揮舞模態(tài)和擺振模態(tài)都是臨界穩(wěn)定的;而前兩階扭轉(zhuǎn)模態(tài)由于對應(yīng)的特征值實部小于零，因此是嚴格穩(wěn)定的。這說明，作用于葉片扭轉(zhuǎn)模態(tài)的氣彈阻尼相對較大，具有較強的抑制扭轉(zhuǎn)模態(tài)顫振的能力。當β0增加時，葉片揮舞模態(tài)明顯地由臨界穩(wěn)定變?yōu)閲栏穹€(wěn)定的，β0的存在對揮舞模態(tài)的穩(wěn)定效果非常明顯，而擺振模態(tài)則由臨界穩(wěn)定變?yōu)椴环€(wěn)定的，并且其不穩(wěn)定性隨著β0的增加而增加，這很可能是由于作用于擺振方向上的氣動阻力非常小，并且同時也沒有考慮結(jié)構(gòu)阻尼的緣故;與此同時，還可以看到，β0對扭轉(zhuǎn)模態(tài)的顫振性能的作用卻非常微弱，不會使扭轉(zhuǎn)模態(tài)的顫振性能產(chǎn)生本質(zhì)的變化。

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

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

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

3 結(jié)論

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

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

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

［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.