李　亮，李映輝，楊鄂川，3

（1.安徽理工大學　理學院力學系，安徽　淮南232001；2.西南交通大學　力學與工程學院，成都610031；3.重慶理工大學　車輛工程學院，重慶400054）

風力機葉片揮舞
—擺振氣彈穩(wěn)定性分析

李亮1，李映輝2，楊鄂川2，3

（1.安徽理工大學理學院力學系，安徽淮南232001；2.西南交通大學力學與工程學院，成都610031；3.重慶理工大學車輛工程學院，重慶400054）

根據(jù)Euler-Bernoulli梁理論和粘彈性材料的Kelvin-Voigt理論建立風力機葉片揮舞—擺振耦合非線性動力學方程。將位移視為靜態(tài)位移和動態(tài)位移的疊加，進而將非線性動力學方程線性化為動態(tài)位移的線性方程，得到葉片耦合振動特征方程。使用基于加權(quán)殘值的Galerkin方法求解特征方程，分析葉片氣彈穩(wěn)定性，討論風速、安裝角、耦合效應(yīng)和材料阻尼對葉片顫振穩(wěn)定性和非線性自激振動行為的影響。結(jié)果表明：擺振方向易出現(xiàn)不穩(wěn)定振動，通過設(shè)置安裝角，利用揮舞—擺振耦合可以控制不穩(wěn)定振動，但當安裝角太大時，揮舞—擺振耦合會引起不穩(wěn)定振動。

振動與波；風力機葉片；氣彈穩(wěn)定性；顫振；自激振動

風力機運轉(zhuǎn)過程中，氣動力、彈性力和慣性力相互耦合，易引起葉片不穩(wěn)定振動。已有對氣彈穩(wěn)定性的研究多基于葉片翼型截面模型：林學海和任勇生［1，2］采用ONERA氣動力模型對葉片揮舞—擺振動力失速穩(wěn)定性進行分析，討論葉片剛度、錐角、結(jié)構(gòu)阻尼等因素對氣彈穩(wěn)定性的影響；劉廷瑞等［3］基于Beddoes-Leishman氣動力模型，使用擬合氣彈系數(shù)的方法簡化計算，研究了葉片動力失速穩(wěn)定性；王偉等［4］建立葉片三維葉型模型和風力機整機模型，采用振型疊加法研究葉片顫振穩(wěn)定性，討論葉片結(jié)構(gòu)參數(shù)對顫振穩(wěn)定性的影響。為獲得較高的風能利用率，葉片通常設(shè)計為變翼型結(jié)構(gòu)。因此，相對于翼型截面模型，連續(xù)系統(tǒng)模型能夠更好地描述葉片整體結(jié)構(gòu)特征：虞心田和崔爾杰［5］建立葉片揮舞—扭轉(zhuǎn)耦合動力學方程，然后將連續(xù)模型簡化為等效的剛體模型，研究葉片的顫振邊界問題；李本立和安玉華［6］考慮葉片揮舞—擺振—扭轉(zhuǎn)耦合非線性振動，使用模態(tài)法求解微分方程，分析葉片氣彈穩(wěn)定性；任勇生和張明輝［7］基于葉片揮舞—擺振—扭轉(zhuǎn)耦合振動偏微分方程和擴展的Theodorsen氣動力模型，使用Galerkin法研究均勻葉片氣彈穩(wěn)定性問題，分析彎曲剛度、扭轉(zhuǎn)剛度、預(yù)錐角等因素對氣彈穩(wěn)定性的影響。

本文研究變截面葉片揮舞—擺振氣彈穩(wěn)定性，討論揮舞—擺振耦合、材料阻尼等已有研究較少考慮的因素對顫振穩(wěn)定性的影響。為考慮非線性影響，將位移分解為靜態(tài)位移和動態(tài)位移，進而將非線性方程線性化為動態(tài)位移的線性方程，得到葉片振動特征方程，使用Galerkin方法分析葉片振動特性，討論葉片顫振穩(wěn)定性和自激振動行為。

1　動力學方程

葉片構(gòu)型如圖1所示，擺振為旋轉(zhuǎn)平面內(nèi)的彎曲，揮舞為垂直于旋轉(zhuǎn)平面的彎曲。旋轉(zhuǎn)坐標系（xyz）用來描述葉片變形，其原點在葉根截面剪切中心，x軸沿彈性軸指向葉尖，y軸沿擺振方向，z軸沿揮舞方向。扭角θ=θ0+θt是弦線與旋轉(zhuǎn)平面的夾角，θ0是安裝角，θt=θt（x）是截面預(yù)扭角。e是葉根到輪轂中心的偏置量，ex和ey是其沿x軸、y軸的分量。假設(shè)葉片截面質(zhì)心、形心和剪切中心三者重合，將細長葉片簡化為Euler-Bernoulli型懸臂梁，考慮其彎彎耦合振動，在彈性梁模型［6-8］的基礎(chǔ)上，通過Kelvin-Voigt理論（σ=Eε+ηdε/dt，σ為應(yīng)力，E為楊氏模量，ε為應(yīng)變，η為粘性系數(shù)）引入材料阻尼，得到動力學方程。

圖1　葉片構(gòu)型及坐標系

其中EI1=EI1（x）和EI2=EI2（x）是第一、第二彎曲剛度，EIw=EIw（x）和EIv=EIv（x）分別是揮舞、擺振剛度，近似表示為

w和v分別是揮舞和擺振位移，撇和點分別表示對x和t求導，F(xiàn)T=∫AEεdA是軸力∫AσdA的彈性部分，A=A（x）是截面面積，m=m（x）是單位長度的質(zhì)量，Ω是葉片旋轉(zhuǎn)角速度，βp是錐角，L是葉片長度，和是材料阻尼力，其表達式如下

（Fy，F(xiàn)z）和（My，Mz）分別是沿（y，z）軸的分布力和力矩，氣動力和重力可作為分布荷載加入方程（1）—（2）右端。選用文獻［6］給出的非定常氣動力模型，其表達式如下式中ρ為空氣密度，a為升力曲線斜率，c=c（x）為弦長，R為旋轉(zhuǎn)半徑，λ為入流速度比，其表達式

v0為輪轂處的風速，λ0=v0/（ΩR）為輪轂處的入流速度比

為轉(zhuǎn)子實度，nb為葉片數(shù)，CDO為阻力系數(shù)，rx=rx（x）是截面x的徑向位置，eA=eA（x）為氣動中心到剪切中心的偏置量（氣動中心在前時為正）。重力表達式

式中g(shù)為重力加速度。

2　非定常氣動力作用下葉片自由振動特性

將位移看作靜態(tài)位移κs與動態(tài)位移κd（κ=w，v）的疊加

將式（8）代入方程（1）—（2），得靜態(tài)位移方程

方程組（9）為關(guān)于x的4階非線性變系數(shù)常微分方程組，其解可由數(shù)值方法（如Euler方法、Runge-Kutta法）求出。動態(tài)位移滿足非線性偏微分方程組

先討論自由振動特性，在（10）中去掉強迫項和動態(tài)位移的非線性項，得到如下線性偏微分方程組

將動態(tài)位移分離變量

式中ω為復數(shù)頻率，其實部為固有頻率，虛部為阻尼系數(shù)，γ為模態(tài)函數(shù)。將式（12）代入式（11），約掉時間項，可得特征方程：

式（13）為4階變系數(shù)常微分方程組，其精確解無法直接求出，本文采用基于加權(quán)殘值的Galerkin法［7，9］求解。將模態(tài)函數(shù)表示為滿足邊界條件且線性無關(guān)的試函數(shù)的線性組合

其中，Nw和Nv分別表示所選揮舞和擺振試函數(shù)的個數(shù)。將（14）代入（13），然后在方程Mμ=0兩端同乘）后對x從0到L積分，得到如下代數(shù)方程組

式中M和K是正定的質(zhì)量和剛度陣，C是阻尼陣，T是矩陣轉(zhuǎn)置，0（Nw+Nv）×1是（Nw+Nv）×1零矩陣。系數(shù)a1，···，aNw，b1，···，bNv不能同時為零，所以矩陣-ω2M+ iωC+K非滿秩，即有

由方程（14）-（16）可求出復數(shù)頻率ω和模態(tài)函數(shù)，阻尼比ζ=Im（ω）/［2Re（ω）］。當ζ為負數(shù)時，葉片出現(xiàn)自激振動和發(fā)散不穩(wěn)定。

3　非定常氣動力作用下葉片非線性自激振動

葉片動能主要集中在前兩階模態(tài)（第1階揮舞和第1階擺振模態(tài)）［10］，使用Galerkin方法對連續(xù)系統(tǒng)1階截斷，設(shè)

式中γf和γl分別是揮舞和擺振第1階模態(tài)函數(shù)，qf和ql是廣義位移（其對時間的一、二次導分別為廣義速度和廣義加速度）。將（17）代入方程（1）-（2）得離散系統(tǒng)

式（18）為非線性常微分方程組，為討論葉片非線性自激振動，在（18）中去掉周期外激勵，得

式（19）可由數(shù)值方法［11-12］或攝動理論進行求解，本文采用數(shù)值積分求解。

4　數(shù)值結(jié)果

選取某MW級風力機葉片為計算模型，葉片長度L=48 m，楊氏模量E=30 GPa，氣動中心到彈性軸距離eA=c/4，葉片密度ρb=1.8×103kg/m3，輪轂高度h0=2 L，錐角βp=2°，偏置量ex=2 m，ey=0 m，其它分布參數(shù)見圖2。空氣密度ρ=1.25 kg/m3，升力曲線斜率a=2 π，阻力系數(shù)CDO=0.016，風速梯度常數(shù)p=0.167，轉(zhuǎn)速Ω=15 r/min。使用4階Runge-Kutta法求解靜態(tài)位移方程，采用Galerkin法求解特征方程，其中揮舞、擺振試函數(shù)均選取為｛（x/L）2，···，（x/L）9｝。

圖3給出了揮舞、擺振第1階阻尼比隨安裝角的變化。若忽略揮舞—擺振耦合效應(yīng)，則揮舞阻尼比隨安裝角絕對值的增大而減小，這是因為截面扭角對揮舞頻率影響較大，頻率隨安裝角幅值的增大而增大，但安裝角對揮舞阻尼系數(shù)影響相對較小。由于預(yù)扭角的存在，揮舞最大阻尼比出現(xiàn)在-10°附近。安裝角對擺振阻尼系數(shù)影響較大，阻尼系數(shù)隨安裝角增加而增大，但安裝角對擺振頻率影響較小，因此擺振阻尼比隨安裝角的增大而增大。忽略耦合效應(yīng)的結(jié)果顯示：安裝角為負（逆時針旋轉(zhuǎn)的葉片，前緣在旋轉(zhuǎn)面前面時為負）且幅值較大時，擺振方向易出現(xiàn)氣彈性不穩(wěn)定性。

考慮揮舞—擺振耦合效應(yīng)影響，由于安裝角和截面預(yù)扭角的存在，擺振方向的能量傳遞到揮舞方向，揮舞阻尼較大，在幅值較大的負安裝角時出現(xiàn)的擺振不穩(wěn)定性被控制。當風速較低時（如v0=8 m/s），揮舞阻尼比變化趨勢與忽略耦合效應(yīng)時相同。當風速較高時（如v0=32 m/s），揮舞阻尼比隨安裝角增大而增大，當安裝角為正且幅值較大時，擺振方向出現(xiàn)不穩(wěn)定振動。兩種情況（忽略和考慮耦合效應(yīng)）的對比結(jié)果顯示：通過設(shè)置安裝角，利用揮舞—擺振耦合能改進葉片氣彈穩(wěn)定性，但當安裝角很大時，耦合效應(yīng)反而會引起氣彈不穩(wěn)定振動。

圖2　葉片幾何和物理參數(shù)

圖3阻尼比隨安裝角的變化

圖4—圖8揭示葉片非線性自由振動行為，計算時選取風速v0=32 m/s，圖4—圖7不考慮材料阻尼。氣彈不穩(wěn)定性出現(xiàn)在擺振方向（見圖3），因此只考慮擺振響應(yīng)。圖4給出安裝角-2°時，忽略耦合效應(yīng)的擺振相圖，初始值選取為［10-50］，結(jié)果顯示擺振響應(yīng)出現(xiàn)發(fā)散不穩(wěn)定現(xiàn)象，響應(yīng)幅值隨時間逐漸增大。圖5給出安裝角20°時，忽略耦合因素的擺振相圖，此時葉片運動穩(wěn)定，擺振響應(yīng)隨時間逐漸收斂于穩(wěn)定的平衡點。圖6給出安裝角-2°時，考慮耦合效應(yīng)的擺振相圖，系統(tǒng)運動穩(wěn)定，擺振相軌跡收斂到穩(wěn)定的平衡點，圖4中出現(xiàn)的發(fā)散不穩(wěn)定通過耦合效應(yīng)被控制。圖7給出安裝角20°時，考慮耦合效應(yīng)的擺振相軌跡，同圖5相比，耦合效應(yīng)使系統(tǒng)運動行為發(fā)生變化，擺振響應(yīng)由穩(wěn)定的平衡點變成大幅值的擬周期運動。圖4—圖7驗證了圖3的結(jié)論。圖8給出安裝角20°時，考慮耦合效應(yīng)的擺振響應(yīng)隨材料阻尼的分岔圖，結(jié)果顯示擺振響應(yīng)隨材料阻尼的增大由大幅值的擬周期運動變?yōu)榉€(wěn)定的平衡點，因此材料阻尼改善了葉片氣彈穩(wěn)定性。

圖4　忽略耦合的擺振相軌跡（θ0=-2°）

5　結(jié)語

圖5　忽略耦合的擺振相軌跡（θ0=20°）

圖6　考慮耦合時擺振相軌跡（θ0=-2°）

基于Euler-Bernoulli梁模型和Galerkin法分析葉片氣彈穩(wěn)定性，結(jié)果表明：通過設(shè)置安裝角，利用揮舞—擺振耦合效應(yīng)可以改善葉片氣彈穩(wěn)定性，但當安裝角太大時，耦合效應(yīng)反而引起葉片不穩(wěn)定振動；材料阻尼可以改善葉片氣彈穩(wěn)定性。

圖7　考慮耦合時擺振時間歷程和龐加萊映射（θ0=20°）

圖8　擺振位移對r=η/E的分岔圖

［1］林學海，任勇生.基于ONERA氣動力模型的風力機葉片顫振時域分析［J］.山東科技大學學報（自然科學版），2009，28（3）：56-60.

［2］任勇生，林學海.風力機葉片揮舞/擺振的動力失速非線性氣彈穩(wěn)定性研究［J］.振動與沖擊，2010，29（1）：121-124.

［3］劉廷瑞，任勇生，楊興華.基于氣彈系數(shù)的氣彈穩(wěn)定性分析［J］.太陽能學報，2010，31（4）：513-516.

［4］王偉，廖明夫，Martin Kühn.風力機葉片三維模型顫振問題［J］.噪聲與振動控制，2013，33（5）：90-96.

［5］虞心田，崔爾杰.分析水平軸風力機葉片氣彈性穩(wěn)定性的簡單方法［J］.太陽能學報，1990，11（1）：52-59.

［6］李本立，安玉華.風力機氣動彈性穩(wěn)定性的研究［J］.太陽能學報，1996，17（4）：314-320.

［7］任勇生，張明輝.水平軸風力機葉片的彎扭耦合氣彈穩(wěn)定性研究［J］.振動與沖擊，2010，29（7）：196-200.

［8］Li L，Li Y H，Lv H W，et al.Flapwise dynamic response of a wind turbine blade in super-harmonic resonance［J］. Journal of Sound and Vibration，2012，331（17）:4025-4044.

［9］李映輝，李中華.超音速下粘彈性夾層壁板顫振分析［J］.力學季刊，2012，33（3）：449-455.

［10］J?rgensen K H，Nielsen S R K.System reduction in multibody dynamics of wind turbines［J］.Multibody System Dynamics，2009，21:147-165.

［11］呂海煒，李映輝，李中華，等.超聲氣流下粘彈性夾層壁板非線性顫振分析［J］.航天器環(huán)境工程，2013，30（1）：40-48.

［12］覃文源，張志誼，華宏星.摩擦力作用下推進軸系彎—扭耦合振動特性分析［J］.噪聲與振動控制，2013，33（5）：4-9.

Flap/Lead-lagAeroelastic Stability Analysis for Wind Turbine Blades

LI Liang1，LI Ying-hui2，YANG E-chuan2，3
（1.Department of Mechanics，College of Science，Anhui University of Science and Technology，Huainan 232001，Anhui China；2.School of Mechanics and Engineering，Southwest Jiaotong University，Chengdu 610031，China；3.ChongqingAutomobile Institute，Chongqing University of Technology，Chongqing 400054，China）

The nonlinear partial differential equations which govern the coupled flap-lead/lag vibration of wind turbine blades were established based on the Euler-Bernoulli beam theory and the Kelvin-Voigt theory for cohesive elastic composite materials.By decomposing blade displacement into static displacement and dynamic displacement，the nonlinear governing equations were linearized to the linear equations for the dynamic displacement.And then，the characteristic equation of the coupled vibration was obtained.The Galerkin method based on the method of weighted residuals was employed to solve the characteristic equation and analyze the aeroelastic stability of the blades.The influence of wind speed，installation angle，coupling effect and material damping on the aeroelastic stability and nonlinear self-excited oscillation of the blades was discussed.Results show that the coupling effect between flap and lead-lag may be employed to improve the aeroelastic stability of the blades，but very strong coupling effect between them due to the large installation angle may bring about vibration instability.

vibration and wave；wind turbine blade；aeroelastic stability；flutter；self-excited oscillation

O317；O322；O323

ADOI編碼：10.3969/j.issn.1006-1335.2015.05.005

1006-1355（2015）05-0030-05

2014－12－16

國家自然科學基金項目（11072204，11372257）；中央高校專項基金（SWJTU11ZT15）；2013年四川省青年科技創(chuàng)新團隊項目；上海市力學在能源工程中的應(yīng)用重點實驗室資助項目；教育部留學回國人員科研啟動基金項目（2012-42）；安徽省自然科學基金項目（1308085MA13）

李亮（1984－），男，山東省曲阜市人，講師，目前從事非線性動力學研究。

E-mail:lli@my.swjtu.edu.cn