風力機復合材料柔性葉片的顫振分析

任勇生，杜向紅，楊樹蓮

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

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

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

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

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

1 基本理論

1.1 位移場和應變場

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

葉片坐標系和主軸坐標系之間，存在下列關系:

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

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

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

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

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

1.2 運動方程

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

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

外力的虛功:

其中:Fx，Fy，Fz分別是作用于x，y，z軸方向的分布氣動力，Mx是繞彈性軸的氣動力矩。

其中:

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

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

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

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

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

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

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

1.3 氣動載荷

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

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

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

1.4 求解方法

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

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

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

葉片的彎曲變形和扭轉變形可以表示如下:

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

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

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

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

2 計算結果與討論

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

在下面的葉片顫振數值計算中，選取葉片的幾何尺寸為:截面寬度0.95 m，高度0.4 m，薄壁厚度h=0.047 5 m，單層厚度h/6，鋪層數為6層;葉片長度R=19 m。對于CAS構型葉片，所采用鋪層方式為:上壁［θ］6，下壁［-θ］6，左壁和右壁［θ/-θ］3。葉片復合材料選取石墨/環氧樹脂，其性能參數為: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。

葉片其他參數取為:

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

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

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

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

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

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

3 結論

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

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

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

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