第一作者徐磊男，碩士生，1988年生

通信作者李德源男，博士后，教授，1965年生

基于非線性氣彈耦合模型的風力機柔性葉片隨機響應分析

徐磊1，李德源1，莫文威1，呂文閣1，劉雄2(1.廣東工業(yè)大學機電工程學院，廣州510006; 2.汕頭大學能源研究所，汕頭515063)

摘要：為準確模擬柔性葉片在紊流風速下的氣彈響應，研究葉片振動對氣動載荷與氣彈響應的反饋,構(gòu)建包含多體系統(tǒng)動力學模型與氣動模型的柔性葉片非線性氣動彈性力學模型。在將細長柔性葉片離散為多剛體系統(tǒng)基礎上，運用計算多體系統(tǒng)動力學理論和Roberson-Wittenburg的建模方法，結(jié)合葉素動量理論，采用Kaimal模型模擬脈動風速，建立風力機柔性葉片的氣彈耦合方程。算例以美國可再生能源實驗室(NREL)研制的5MW近海風力機為研究對象，分析葉片的振動和葉根的揮舞與擺振力矩，研究柔性葉片振動對氣動載荷影響。結(jié)果表明,葉片達到一定長度后，模擬風力機氣彈響應問題時，其振動影響不可忽略。

關鍵詞：風力機；氣彈響應；氣彈耦合；振動

基金項目：國家自然科學基金(51276043)；教育部高等學校博士點科研

收稿日期：2014-01-08修改稿收到日期：2014-11-10

中圖分類號:TK83文獻標志碼：A

Random response analysis for flexible blade of a wind turbine based on nonlinear aero-elastic coupled model

XULei1,LIDe-yuan1,MOWen-wei1,LüWen-ge1,LIUXiong2(1. School of Electromechanical Engineering, Guangdong University of Technology, Guangzhou 510006, China;2.Institute of Energy Science, Shantou University, Shantou 515063, China)

Abstract:In order to exactly model the aeroelastic response of a blade under turbulent flow and study blade vibration’s feedback on aerodynamic load and aeroelastic response, here a nonlinear aeroelastic mechanical model containing flexible blade multibody system’s dynamic model and an aerodynamic model was constructed. Through discreting a flexible slender blade into a multi-rigid-body system, the aeroelastic coupled equations of the blade were obtained by using the theory of computational dynamics of multibody system and Roberson-Wittenburg modeling method. In the coupled equations, the blade element-momentum theory (BEM) was used to calculate the aeroelastic load of the wind turbine blade and Kaimal model was used to calculate random wind speed. The time domain response of a 5-MW offshore wind turbine blade of the US. Renewable Energy Laboratory was calculated. The response included the vibration of the blade and the flapwise and edgewise bending moments at the blade root. The effect of the vibration of the flexible blade on the aerodynamic load was also studied. The results revealed that the effect of slender blade vibration can not be ignored when simulating the aeroelastic response of a wind turbine.

Key words:wind turbine; aeroelastic response; aeroelastic coupled; vibration

大型風力機葉片一般采用復合材料蒙皮，其外形尺寸不斷增大，作為主要承受空氣動力的構(gòu)件，細長葉片的有限變形和振動對于風力機組氣動力的影響以及由此引起的與葉片結(jié)構(gòu)響應的耦合成為不能忽略的問題。尋求包含非線性變形效應的柔性葉片動力特性與氣動彈性耦合的數(shù)值分析方法成為了風力機氣彈分析的重要基礎[1]。

對于風力機系統(tǒng)的力學建模，目前國內(nèi)外常用的建模方法主要有多體系統(tǒng)方法 (Multibody Systems，MBS)、有限元方法(Finite Element Systems，F(xiàn)ES)和連續(xù)系統(tǒng)(Continuous Systems，COS)方法等[2-4]。對于風力機葉片來說，COS方法所建立的偏微分方程組僅在特殊的、簡單的幾何結(jié)構(gòu)和載荷下才可能求解，不太適合復雜系統(tǒng)的時域仿真；而FES方法具有較多的自由度，計算和分析成本較高，僅適合于靜態(tài)載荷分析和微動力學分析；MBS方法是將實際的機械構(gòu)件視為剛體，用有限的自由度導出系統(tǒng)動力學微分方程(組)，但當系統(tǒng)中柔性構(gòu)件的變形對系統(tǒng)動力學行為產(chǎn)生較大影響時，該方法模擬的精度有限?；旌隙囿w系統(tǒng)(Hybrid Multibody System，HMBS)方法則考慮了系統(tǒng)中構(gòu)件的柔性，適合于包含剛體和柔性體的復雜系統(tǒng)的力學建模，也適合于風力機葉片或整機的力學建模。

對于風力機葉片這種具有空間運動和彈性變形的構(gòu)件，Rauh等[5-7]引入“超級單元”(Superelement)對其進行離散，以反映其彈性變形,將LW-50/750風力機中變形相對較大的構(gòu)件(如風輪葉片和塔架)處理為柔性體，將變形相對較小的構(gòu)件(如輪轂、機艙)處理為剛體，采用超級單元(無扭轉(zhuǎn)自由度)對風力機的柔性構(gòu)件進行離散，各剛體之間通過力元(彈簧和阻尼器)和鉸連接。這樣風力機系統(tǒng)可以用一離散的剛、柔混合多體系統(tǒng)來表示。采用超級單元(含扭轉(zhuǎn)自由度)對600 kW風力機的柔性葉片及塔架進行離散處理，基于虛功原理導出了整機混合多體系統(tǒng)的拉格朗日方程的顯式表達式并分析了整機的固有模態(tài)。使用HMBS方法建立了NM80風力機的力學模型，利用哈密頓原理導出了動力學方程；分析了超級單元個數(shù)及計算彈簧系數(shù)不同的方法對柔性件的固有頻率和阻尼系數(shù)計算精度的影響。

對于葉片的空氣動力計算，目前通過數(shù)值求解Navier-Stokes 方程的計算流體力學(Computational Fluid Dynamics, CFD)方法[8]，計算精度較高，但對于風力機葉片這種較大尺寸且旋轉(zhuǎn)運動的構(gòu)件，計算量較大，要實現(xiàn)實時的氣彈耦合還有一定困難。在風力機氣動載荷計算中，葉素動量理論(Blade Element-Momentum Theory, BEM)[9]計算結(jié)果較合理且求解快速，還可與其它非定常氣動模型結(jié)合，已在GH BLADED等大多數(shù)風力機性能與氣動分析軟件中廣泛應用。其它用于氣動載荷分析的方法如廣義動態(tài)尾流理論(Generalized Dynamic Wake Theory ,GDW)等[10]也逐漸受到關注。相比于BEM理論，GDW除了計算旋轉(zhuǎn)平面到尾流的軸向誘導速度還考慮了尾跡效應，氣動模型相對更復雜。

對于風力機葉片或整機的非線性氣彈分析，近年來國內(nèi)外學者也進行大量的研究工作。任勇生等[11-12]基于葉片的典型截面模型和ONERA非定常氣動模型對風力機葉片揮舞/擺振非線性氣彈系統(tǒng)的動力失速穩(wěn)定性進行了研究。劉廷瑞等[13-14]運用連續(xù)梁模型，通過方程的線性化和Galerkin數(shù)值求解方法，對葉片的非線性氣彈問題進行了研究。李映輝等[15]同樣通過建立連續(xù)梁的偏微分方程，施加外部和內(nèi)部的正弦激勵，通過Galerkin數(shù)值分析方法，研究了葉片的非線性振動。Holierhoek[7]利用WOBBE程序研究了由失速和附著流引起的葉片揮舞-擺振不穩(wěn)定性的可能性。Meng[10]則通過開發(fā)結(jié)構(gòu)和氣動接口程序?qū)崿F(xiàn)了基于現(xiàn)有的非線性多剛體動力學程序(MBDyn[16])和空氣動力學模塊(AeroDyn[17])的耦合求解。Ahlstr?m[18]則使用MSC.Marc建立Alsvik水平軸風力機的有限元模型用于計算其動力學響應；采用現(xiàn)有的基于BEM理論的空氣動力學程序AERFORCE計算氣動載荷，同時使用SOSIS-W軟件來模擬不同的風場條件，實現(xiàn)了風力機的氣彈耦合計算。Hansen[1]采用CFD方法計算了某風力機的風輪空氣動力學，研究了上游方向風輪的尾渦形成過程及其對下游方向風輪的影響。

本文應用計算多體動力學理論，將柔性葉片離散為若干通過運動副聯(lián)接的多體系統(tǒng)，基于動力學普遍方程導出了該系統(tǒng)的動力學方程的一般形式，編制了通用的動力學建模與仿真程序。該程序能根據(jù)用于描述系統(tǒng)拓撲構(gòu)型的輸入數(shù)據(jù)自動建立和數(shù)值求解受約束的系統(tǒng)動力學方程(微分-代數(shù)方程組)。通過對系統(tǒng)的脈沖響應與系統(tǒng)的傳遞函數(shù)分析，得到系統(tǒng)的動力特性。在葉片多體系統(tǒng)模型的基礎上，結(jié)合BEM氣動模型導出氣彈耦合模型，在數(shù)值積分的每一時間步中,氣動模型也能實時地根據(jù)葉片各剛體的廣義坐標和廣義速度計算出作用在葉片各剛體上的氣動載荷，并將其反饋到葉片動力學微分方程中，從而實現(xiàn)葉片的非線性氣彈耦合時域響應分析。

本文針對美國可再生能源實驗室(NREL)研制的5MW變速恒頻風力機葉片[19]，數(shù)值分析了葉片在紊流風速下?lián)]舞、擺振與扭轉(zhuǎn)振動和葉根力矩的時域響應，重點關注了柔性葉片的振動對其氣動載荷的影響。相關研究在大型風力機葉片設計與優(yōu)化以及保障風力機組的安全穩(wěn)定運行具有重要應用價值和意義。

1柔性葉片的多體系統(tǒng)動力學建模

大型葉片由于其較大的展弦比，且外殼一般采用玻璃鋼等復合材料制成，在氣動與機械載荷作用下，葉片在垂直來流和沿來流方向產(chǎn)生較大的彎曲變形并伴隨扭轉(zhuǎn)變形，表現(xiàn)出較大的柔性，與氣動載荷形成反饋。本文借助 “超級單元”的思想將葉片離散為通過帶彈簧和阻尼器的萬向節(jié)或轉(zhuǎn)動鉸相聯(lián)接的多剛體系統(tǒng)[7]，每個超級單元包含四個剛體。這樣，葉片的橫向彎曲變形與扭轉(zhuǎn)變形就可以用較少的自由度得到展現(xiàn)。結(jié)合現(xiàn)代計算技術，采用計算多體動力學中的R-W建模方法來建立葉片的多體動力學模型。

1.1超級單元模型和葉片的離散

對于葉片這種細長構(gòu)件，見圖1。將其離散為若干個超級單元。每個單元中，用B1和B2，B3和B4四個通過萬向節(jié)連接的剛體來反映其在相互垂直面內(nèi)的彎曲變形，其廣義坐標為θ1,θ2,θ4,θ5；而扭轉(zhuǎn)變形則由B2和B3間的旋轉(zhuǎn)鉸來反映，廣義坐標θ3。剛體間的相對運動通過彈簧與阻尼器約束(如圖中的Cy1～Cy3、Cz1～Cz3)。

圖1　超級單元模型 Fig.1 Superelement model

美國可再生能源實驗室NREL推出的5MW近海風力機葉片長61.5 m，風輪直徑125.0 m，本文將該葉片離散為4個超級單元，其中前后兩個超級單元中相鄰的剛體固結(jié)為一個剛體(如B4、B7和B10)，這樣，整個葉片的拓撲構(gòu)型為13個剛體組成的多體系統(tǒng)，共21個自由度，見圖2。慣性坐標系為OXYZ(固定坐標系)，為方便考察葉片的揮舞和擺振方向的變形，建立葉根隨動坐標系O′X′Y′Z′，OO′表示輪轂中心與葉根的距離，葉片繞固定坐標系的Z軸即風輪主軸的軸線轉(zhuǎn)動。

1.2柔性葉片的多體系統(tǒng)動力學方程

R-W方法的特點是應用圖論中的概念來描述多剛體系統(tǒng)，將系統(tǒng)中各對相鄰剛體間的相對定位參數(shù)作為描述系統(tǒng)位形的廣義坐標。如圖2中對葉片多剛體系統(tǒng)進行規(guī)則標號，共有13個剛體(Bi,i=1,…,13)和內(nèi)接鉸(Hi,i=1,…,13)。除鉸H1建立在定坐標系O點，其余內(nèi)接鉸Hi的鉸點均建立在各剛體質(zhì)心截面的氣動中心(弦長1/4處)。體B1繞Z軸以風輪轉(zhuǎn)速轉(zhuǎn)動，即存在運動約束；內(nèi)接鉸Hi的相對運動可通過固結(jié)在體Bi及其內(nèi)接體上的相對坐標系來描述。

圖2　剛體規(guī)則標號、慣性坐標系XYZ和葉片坐標系X′Y′Z′ Fig.2 Regular mark number of each rigid body, inertial coordinate system XYZ and blade coordinate system X′Y′Z′

如上所述，可得葉片多體系統(tǒng)的廣義坐標陣為

q(t)=(q1,…,q21)T

(1)

(2)

(3)

文獻[20]中給出了廣義質(zhì)量陣Z、廣義力陣z的具體表達形式和導出方法，由于篇幅所限，本文不再贅述。該方程組中，廣義坐標q描述了剛體間的相對轉(zhuǎn)動，由于柔性葉片較大的變形，這些量一般不為小量，不適合進行線性化，而且由于葉片的旋轉(zhuǎn)，廣義質(zhì)量陣等隨時間變化，即式(3)描述的是非線性、時變的系統(tǒng)。方程組的求解有一定困難，一般采用適當?shù)臄?shù)值積分方法結(jié)合相應的初始條件求得數(shù)值解。

1.3模型的相關參數(shù)與模型驗證

葉片離散為多剛體系統(tǒng)，各剛體之間通過帶彈簧和阻尼器的萬向節(jié)或旋轉(zhuǎn)鉸連接。彈簧和阻尼器的設置是為了約束剛體間的相對轉(zhuǎn)動，彈簧的剛度系數(shù)的取值應當使建立的葉片多剛體模型的力學特性與柔性葉片的盡量符合，一般應該滿足：①在相同的靜載荷作用下，葉片多剛體模型與原柔性葉片的彈性變形相同；②超級單元應與相同尺寸的剛性梁有相同的質(zhì)量和慣性性質(zhì)；③葉片多剛體模型的動力特性應盡量接近原柔性葉片的動力特性，如前若干階固有頻率、振型等相差滿足要求的精度。

分析計算中涉及的葉片相關參數(shù)如沿葉片軸的剛度、質(zhì)量分布，葉片截面幾何參數(shù)以及與氣動載荷計算相關的翼型翼型氣動特性等參數(shù)可參考NREL發(fā)布的文獻[19]。由于葉片的幾何形狀與力學特性可以用彈性梁來描述，參考文獻[7]中列出了各超級單元中彎曲和扭轉(zhuǎn)彈簧的剛度系數(shù)計算表達式。針對該5MW風力機葉片，計算得到的彈簧剛度系數(shù)見表1。在所建立的葉片多體模型基礎上，通過施加脈沖激勵，由頻響函數(shù)(FRF)求得了該葉片的前5階固有頻率，該結(jié)果與應用多體動力學軟件ADAMS對該葉片的分析結(jié)果以及NREL提供的分析結(jié)果進行了對比，結(jié)果見表2。

表1　NREL-5MW風力機葉片彈簧剛度系數(shù)/(Nm/rad)

表2　計算結(jié)果對比

注：單位為Hz,“—”表示沒有數(shù)據(jù)

從表2可以看出，本文計算的葉片前5階頻率與NREL提供的結(jié)果基本吻合，表明文中的葉片離散模型與實際葉片的力學特性基本相同，質(zhì)量、剛度分布也基本與實際相符；也證明了超級單元法能用較少的自由度準確地描述柔性葉片的彈性變形、氣動載荷和葉片旋轉(zhuǎn)運動間的耦合。結(jié)果也驗證了仿真程序的有效性與葉片動力學模型的正確性。

2帶約束項的柔性葉片多體動力學方程的建立及其數(shù)值求解

葉片在受力運動過程中，要受到各種幾何或運動約束，如：控制系統(tǒng)所加的對風輪轉(zhuǎn)速或主軸轉(zhuǎn)矩的控制等，使得描述葉片多體系統(tǒng)位形的坐標q(t)=(q1,…,qn)T(n=21)相互不獨立，若系統(tǒng)獨立的約束方程個數(shù)為s個，即

Φ=(Φ1,…,Φs)T=0

(4)

則系統(tǒng)的獨立的變量中只有δ=n-s個為獨立的。結(jié)合前面導出的動力學方程(3)，得到帶拉格朗日乘子葉片系統(tǒng)的動力學方程

(5)

式中：Φq為約束方程的雅克比s×n階矩陣;σ=(σ1,…,σs)T為與約束方程對應的拉格朗日乘子陣;ΦTqσ表示剛體間所有理想約束力的貢獻。此方程與約束方程(4)一起構(gòu)成封閉的微分-代數(shù)方程組。

采用違約穩(wěn)定方法求解(5)式與(4)式聯(lián)合的微分-代數(shù)方程組，將全部廣義坐標與拉格朗日乘子作為未知變量，方程組變?yōu)檩^大變量數(shù)的封閉方程。對違約現(xiàn)象進行修正，其思想在于對約束方程進行人工反饋[20]，方程為

(6)

式中:α,β為常數(shù)。適當選擇這兩個參數(shù)可保證Φ=0的解漸近穩(wěn)定。將式(6)與式(5)組合成矩陣形式

(7)

3葉片氣動載荷計算與氣彈耦合分析

作用在葉片上的隨機空氣動力的計算采用葉素-動量理論(BEM)，結(jié)合普朗特葉尖修正理論對氣流繞葉尖的損失[8]。假設某剛體Bi(i=1,…,13)所受氣動載荷沿葉片展向均勻分布，作用在葉片截面弦上距離前緣1/4即氣動中心處,見圖3。在葉素理論中，忽略氣流沿葉片展向的流動，此時，作用在葉片旋轉(zhuǎn)平面內(nèi)某一環(huán)形區(qū)域上的推力FT及轉(zhuǎn)矩Q為

FT=0.5BρcW2[CL(α)cosφ+CD(α)sinφ]L

(8)

Q=0.5BρcW2[CL(α)sinφ-CD(α)cosφ]Lr

(9)

進一步考慮氣流的三維流動，根據(jù)動量理論，此時，環(huán)形區(qū)域所受的推力FT及轉(zhuǎn)矩Q表達為

(10)

Q=4πr3ρU∞ω(1-a)a′FL

(11)

式中:c為Bi的質(zhì)心截面處的弦長；W為相對風速，表示剛體Bi質(zhì)心截面氣動中心處的轉(zhuǎn)速ωr與來流風速U∞的合速度；α為剛體Bi質(zhì)心截面處的攻角；L為剛體Bi的長度；ρ為空氣密度；CL、CD和CM分別是升力系數(shù)、阻力系數(shù)和力矩系數(shù)，由翼型特性與攻角確定；r為剛體Bi質(zhì)心截面氣動中心距離旋轉(zhuǎn)中心的距離；B為風輪葉片數(shù)；ω為葉片轉(zhuǎn)速；a、a′分別為軸向、切向誘導因子；F為普朗特葉尖損失因子。

圖3　剛體B i質(zhì)心截面上的氣動參數(shù)及氣動中心連體基 Fig.3 Aerodynamic parameters on a CM cross-section of rigid body B iand body-fixed coordinate system aerodynamic center

本文研究的重點之一在于分析葉片變形對氣動載荷的影響，這種影響體現(xiàn)在葉片揮舞與擺振速度對入流角φ及葉片扭轉(zhuǎn)對攻角α的影響

(12)

α=φ-θ

(13)

式中:φ為入流角(圖3)；θ為葉片截面的扭轉(zhuǎn)角；Ve-op、Ve-ip分別為氣動中心揮舞與擺振速度。

應用以上表達式進行葉片氣動載荷計算時，其中的軸向誘導因子a、切向誘導因子a′，以及入流角φ和攻角α需要根據(jù)翼型靜態(tài)升、阻力數(shù)據(jù)通過迭代計算，具體迭代步驟見文獻[21]。

葉片所受空氣動力與其彈性變形的相互影響即氣彈耦合分析流程見圖4。圖中顯示了葉片多體系統(tǒng)動力分析模塊與其氣動載荷計算模塊的耦合過程，在求解動力學微分方程數(shù)值積分的每一個時間步中，調(diào)用氣動模塊計算出該時刻的氣動載荷，然后，這些載荷返回到數(shù)值積分模塊，計算出下一時間步中各廣義坐標與廣義速度的值，這些參數(shù)再輸入到氣動模塊，如此循環(huán)直到時間設定的仿真時間tend。

在本文的分析中，假設葉片以勻角速度繞主軸轉(zhuǎn)動，為此建立約束方程

Φ=q1-ωet=0

(14)

式中:廣義坐標q1為葉根剛體相對于主軸的轉(zhuǎn)動;ωe為額定轉(zhuǎn)速。如果葉片轉(zhuǎn)速隨風速變化，如變轉(zhuǎn)速控制風力機，(14)式變?yōu)榉蔷€性約束方程。數(shù)值積分采用違約穩(wěn)定法，取α=β=10，積分精度可視要求由分析者設定，本文算例仿真分析時取為10-5。

圖4　氣彈耦合分析流程圖 Fig.4 Flow chart of aeroelastic coupling analysis

4紊流風場風速的模擬

從風力機風場大量實測風的數(shù)據(jù)可以知道，風速由兩部分組成。第一部分是長周期部分，是大氣宏觀整體運動形成的，其周期一般在10 min以上，方向一般為水平縱向，大小只與高度有關；第二部分是短周期部分，是局部的紊流運動形成的，其周期只有幾秒至幾十秒。那么對于實際風場，風力機風場模型可分為穩(wěn)態(tài)風和脈動風。

根據(jù)統(tǒng)計大量的風速實測樣本可知，脈動風是零均值平穩(wěn)的隨機過程?？梢圆捎米曰貧w(Auto-Regressive，AR)模型模擬風速時程，那么M個點空間相關脈動風速時程過程[V(t)]=[V1(t),…,VM(t)]T可由下式生成

(15)

式中:[V(t-kΔt)]=[V1(t-kΔt),…,VM(t-kΔt)]T；[ψk]為AR模型自回歸系數(shù)矩陣，M×M階方陣；k=1,…,p，p為AR 模型的階數(shù)；Δt為時間步長；N(t)為獨立隨機過程向量。

模擬風速時程的問題歸結(jié)為回歸系數(shù)矩陣[ψk]的求解、給定隨機方差[N(t)]的求解和AR模型階數(shù)的確定。具體分析過程見文獻[22]。

風速譜采用Kaimal譜。該譜考慮了脈動風速功率譜隨高度的變化，其縱向脈動風速功率譜表達式為

(16)

式中:f=nz/V(z)，V*=KV(z)/ln(z/z0),K≈0.4，V*為摩擦風速，V(z)為距地面高度z處的平均來流風速[22]。

頻率取樣點數(shù)取N=512，時間間隔Δt=0.1s，風速模擬時間為50 s，AR模型階數(shù)p=4。5 MW大型水平軸風力發(fā)電機輪轂高度為90 m，圖5為輪轂中心平均風速V(z)=11.4 m/s時中心處隨機紊流風速時程曲線。

圖5　輪轂處風速時間歷程曲線 Fig.5 The time history curves of wind speed at the hub

5算例分析

5.1模型及其參數(shù)的選取

基于以上的建?；A，對文獻[19]提供的5 MW風力機葉片在隨機風速下的氣彈響應進行了數(shù)值分析。計算過程在隨機風速及額定轉(zhuǎn)速(12.1 r/min)下進行。

考慮到此風輪直徑(126 m)較大，葉片上各剛體(葉素)在不同的高度風速存在風剪，引入風速隨高度變化的公式加以修正

(17)

式中:V,V0為距地面高度分別為H及H0處的平均風速。本文取V0=11.4 m/s，為輪轂高度H0=90 m處的風速；α為切變系數(shù)或粗糙度指數(shù)，是一個經(jīng)驗指數(shù)，其值為0.125～0.6，文中取0.16。

5.2葉片的氣彈響應分析

圖6中(a)和(b)分別為葉根的揮舞和擺振角加速度。由于葉片揮舞方向的彈簧剛度系數(shù)比擺振方向小，故葉片揮舞角加速度明顯比擺振角加速度大。雖然模擬風速的瞬時波動會迅速加劇葉片的振動(如圖6(b)中，擺振角加速度在35 s附近突然增大，這與圖5中的模擬風速在該時刻迅速增加有關)，但從整體趨勢來看，結(jié)構(gòu)阻尼的存在使得葉片在這兩個方向的角加速度慢慢地減少。

圖6　 葉片揮舞角加速度和 擺振角加速度的時間歷程曲線 Fig.6 The time history curves of flapwise and edgewise angular acceleration

圖7　葉根揮舞和擺振力矩時間歷程曲線 Fig.7 The time history curves of flapwise and edgewise bending moment at the blade root

圖7的(a)和(b)圖分別為葉片葉根處的揮舞力矩和擺振力矩隨時間變化曲線。在氣動推力的持續(xù)作用下，葉根處承受著較大的揮舞力矩，均值約為5×106N·m。而在擺振方向，葉根同樣承受著較大的力矩(約2×106N·m)，其中一部分來自于氣動切向力產(chǎn)生的氣動力矩，而主要部分則來自于葉片的自重。葉片在旋轉(zhuǎn)過程中，葉片重力的方向會周期性地改變，導致葉片的擺振力矩變化幅度比揮舞力矩大，所以葉根處很容易產(chǎn)生疲勞破壞。

圖8　葉尖軸向和切向誘導因子時間歷程曲線 Fig.8 The time history curves of axial and tangential inducing factors at the blade tip

圖9　剛體B 13(葉尖)質(zhì)心截面上攻角的時間歷程曲線 Fig.9 The time history curve of attack angle in the center-mass section of rigid body B 13 (blade tip )

圖8(a)和(b)分別為葉尖軸向誘導因子和切向誘導因子隨時間的變化情況。圖8(a)中軸向誘導因子在0.15～0.25左右變化，小于0.4；圖8(b)中切向誘導因子較小，僅在0.025左右變化。說明在風力機運行過程中，尾流較穩(wěn)定，其旋轉(zhuǎn)速度較小。

葉尖是葉片變形最大的位置，需要了解其氣動參數(shù)隨時間的變化情況。圖9為剛體B13(葉尖)質(zhì)心截面上的攻角的時間歷程曲線。由于葉片的各種振動的疊加作用，導致葉尖攻角產(chǎn)生波動，但大約保持在5°～9°范圍內(nèi)。

圖10(a)和(b)分別為葉片各剛體在葉片坐標系X′方向(揮舞)和Z′方向(擺振)的位移。剛體B1(葉根)代表葉片上與輪轂固結(jié)的剛體，其運動被約束，故X′方向和Z′方向位移均為零，距離風輪主軸越則位移越大。從圖10(a)中可以看出葉片在X′方向由于持續(xù)受到氣動推力作用而產(chǎn)生變形，葉尖處變形最大，約為5 m左右。當風力機進入到穩(wěn)定運行階段后，由于攻角的變化，引起氣動載荷的改變，導致葉片出現(xiàn)揮舞振動。圖10(b)中葉片在擺振方向受到了由氣動載荷在風輪旋轉(zhuǎn)平面內(nèi)產(chǎn)生的切向力的作用，故葉片在旋轉(zhuǎn)平面內(nèi)的振動較為頻繁。

圖10 各剛體揮舞與擺振位移(葉片坐標系) Fig.10 The flapwise and edgewise displacement of each rigid body (in blade coordinate system)

6結(jié)論

(1)本文針對大型風力機柔性葉片的氣彈耦合進行研究。由于葉片動力學微分方程的非線性與耦合的復雜性，文中應用超級單元法將柔性葉片離散為有限個剛體組成的多體系統(tǒng)，用較少的自由度來模擬這種耦合。應用計算多體動力學理論，建立柔性葉片多體系統(tǒng)的非線性代數(shù)-微分氣彈耦合方程，結(jié)合適當?shù)臄?shù)值求解方法，實現(xiàn)了柔性葉片的非線性氣彈耦合時域響應分析。算例對5 MW的風力機在紊流風速下的非線性氣彈響應進行了數(shù)值模擬。

(2)通過MATLAB仿真平臺，編制了柔性葉片氣彈耦合模型的構(gòu)建與數(shù)值求解的程序，使得葉片多剛體模型的構(gòu)建與求解的程序化，在葉片設計或優(yōu)化階段，僅需調(diào)整輸入?yún)?shù)，就能得到不同的優(yōu)化結(jié)果。對于變速或定轉(zhuǎn)速風力機，通過建立相應的轉(zhuǎn)速約束條件，可以實現(xiàn)相應的運動分析。除葉片的運動分析外，增加內(nèi)力分析模塊，可以分析葉片危險截面內(nèi)力和變形，實現(xiàn)在特殊風況條件下葉片的強度、剛度校核。在葉片適當位置施加脈沖激勵，通過傳遞函數(shù)分析，可以得出葉片的根跡圖，以及葉片動力特性如固有頻率、振型、氣彈阻尼等與氣彈穩(wěn)定性相關的分析。在輸入文件中，通過增加或減少超級單元的個數(shù)，可對分析精度或計算量進行調(diào)整。

(3)算例分析表明，長大葉片在隨機風載荷作用下，其揮舞與擺振等均為大幅值的非線性振動，且隨著風力機單機的功率不斷增大、葉片的柔性不斷增加，葉片的振動對葉片的氣動力有強烈的反饋，二者間相互影響即氣彈耦合越加明顯。

(4)大型風力機柔性葉片在紊流作用下的氣彈耦合時域響應分析，不僅在風力機葉片的設計階段，同時在風力機運行過程中，為保證其安全穩(wěn)定運行均有重要的應用價值，也是進行葉片氣彈穩(wěn)定性分析的重要基礎。本文的工作為進一步研究柔性葉片非線性氣彈失穩(wěn)機理與氣彈失穩(wěn)影響因素，構(gòu)建整機的氣彈耦合分析模型提供有參考價值的思路和平臺。

參考文獻

[1]Hansen M O L, S?rensen J N, Voutsinas S, et al. State of the art in wind turbine aerodynamics and aeroelasticity[J]. Progress in Aerospace Sciences,2006, 42: 285-330.

[2]Molenaar D P. Modeling and control of the NedFlex turbine-NedFlex: a flexible, variable rotational speed wind turbine[R]. Mechanical Engineering, Systems and Control Group, Delft University of Technology, The Netherlands, Technical Report TUD-WBMR-A-746,August,1996.

[3]Molenaar D P. Modeling the structural dynamics of the Lagerwey LW-50/750 turbine. Wind Engineering,1998,22(6):253-264.

[4]Molenaar D P. Cost-effective design and operation of variable speed wind turbines[D]. Ph.D. Thesis: Delft University of Technology,2003.

[5]Rauh J, Schiehlen W. Various approaches for the modeling of flexible robot arms[C]// Proceedings of the Euromech-Colloquium 219 on Refined Dynamical Theories of Beams, Plates, and Shells and their Applications, Kassel, 1986, Springer-Verlag, Berlin Heidelberg,1987:420-429.

[6]Zhao X Y, Peter M, Wu J Y. A new multibody modelling methodology for wind turbine structures using a cardanic joint beam element. Renewable Energy,2007, 32:532-546.

[7]Holierhoek J G. Aeroelasticity of large wind turbines[D]. Ph.D. Thesis: Delft University of Technology,2008.

[8]Meng F Z.Aero-elastic stability analysis for large-scale wind turbines[D]. Ph.D. Thesis:Delft University of Technology,2011.

[9]Hansen M O L. Aerodynamics of wind turbines(2nd ed)[M].UK: Earthscan,2008.

[10]Kim S, Sclavounos P D. Fully coupled response simulations of theme offshore structure in water depths of up to 10 000 feet[C]. In: Proceedings of 11th International Offshore and Polar Engineering Conference (ISOPE), Stavanger, Norway, 2005: 457-466.

[11]任勇生,林學海.風力機葉片揮舞/擺振的動力失速非線性氣彈穩(wěn)定性研究[J].振動與沖擊,2010,29(1): 121-124.

REN Yong-sheng, LING Xue-hai. Flap/lead-lag nonlinear aeroelastic stability of a wind turbine blade system during dynamic stall[J]. Journal of Vibration and Shock,2010,29(1): 121-124.

[12]任勇生,劉廷瑞.具有結(jié)構(gòu)阻尼的復合材料薄壁梁的動力失速非線性顫振特性[J].振動與沖擊,2013, 32(18): 146-152.

REN Yong-sheng, LIU Ting-rui. Stall nonlinear flutter behavior of a thin-walled composite beam with structural damping[J]. Journal of Vibration and Shock,2013, 32(18): 146-152.

[13]任勇生,劉廷瑞,楊樹蓮,等.風力機復合材料葉片的動力失速氣彈穩(wěn)定性研究[J].機械工程學報,2011, 47(12): 113-125.

PEN Yong-sheng, LIU Ting-rui, YANG Shu-lian. Aeroelastic Stability Analysis of Composite Wind Turbine Blade Dynamic Stall[J]. Journal of vibration engineering,2011, 47(12): 113-125.

[14]劉廷瑞,任勇生.基于復合材料薄壁結(jié)構(gòu)的轉(zhuǎn)子葉片非線性氣彈時域響應分析[J].太陽能學報,2012,33(1): 105-112.

LIU Ting-rui, REN Yong-sheng. Nonlinear aeroelastic response analysis of rotor blade modeled as composite thin-walled structure [J].Acta Energiae Solaris Sinica,2012,33(1):105-112.

[15]Li L, Li Y H, Liu Q K, et al. Flapwise non-linear dynamics of wind turbine blades with both external and internal resonances[C]. International Journal of Non-Linear Mechanics, 2014.

[16]Masarati P. Comprehensive multibody aero servo elastic analysis of inte-grated rotorcraft active controls[C]. Ph D thesis, Politecnio di Milano, Di-partimento di Ingegneria Aerospaziale, Jul,2000.

[17]Laino D J,Hansen A C. User’s guide to the wind turbine aerodynamics computer software. Technical report, Windward Engineering[D]. LC, Salt Lake City, UT 84117, December,2002.

[18]Ahlstr?m A. Aeroelastic simulation of wind turbine dynamics[D]. Ph.D. Thesis: Royal Institute of Technology; 2005.

[19]Jonkman J, Butterficle S, Msuial W, et al. Definition of a 5-MW reference wind turbine for offshore system development[R]. Springfield: U.S.National Rene-wable Energy Laboratory, 2009.

[20]洪嘉振.計算多體系統(tǒng)動力學[M]. 北京：高等教育出版社，1999.

[21]李春，葉舟，高偉,等.現(xiàn)代陸海風力機計算與仿真[M]. 上海：上海科學技術出版社，2012.

[22]袁波,應惠清,徐佳煒.基于線性濾波法的脈動風速模擬及其MATLAB程序的實現(xiàn)[J].結(jié)構(gòu)工程師，2007,23(4):55-61.

YUAN Bo, YING Hui-qing, XU Jia-wei. Simulation of turbulent wind velocity based on linear filter method and MATLAB program realization [J].Structural Engineers, 2007, 23(4):55-61.