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

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

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

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

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

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

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

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

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

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

1　動力學(xué)方程

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

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

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

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

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

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

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

式中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）中去掉強(qiáng)迫項(xiàng)和動態(tài)位移的非線性項(xiàng)，得到如下線性偏微分方程組

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

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

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

其中，Nw和Nv分別表示所選揮舞和擺振試函數(shù)的個(gè)數(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不能同時(shí)為零，所以矩陣-ω2M+ iωC+K非滿秩，即有

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

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

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

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

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

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

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

選取某MW級風(fēng)力機(jī)葉片為計(jì)算模型，葉片長度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，風(fēng)速梯度常數(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)，則揮舞阻尼比隨安裝角絕對值的增大而減小，這是因?yàn)榻孛媾そ菍]舞頻率影響較大，頻率隨安裝角幅值的增大而增大，但安裝角對揮舞阻尼系數(shù)影響相對較小。由于預(yù)扭角的存在，揮舞最大阻尼比出現(xiàn)在-10°附近。安裝角對擺振阻尼系數(shù)影響較大，阻尼系數(shù)隨安裝角增加而增大，但安裝角對擺振頻率影響較小，因此擺振阻尼比隨安裝角的增大而增大。忽略耦合效應(yīng)的結(jié)果顯示：安裝角為負(fù)（逆時(shí)針旋轉(zhuǎn)的葉片，前緣在旋轉(zhuǎn)面前面時(shí)為負(fù)）且幅值較大時(shí)，擺振方向易出現(xiàn)氣彈性不穩(wěn)定性。

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

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

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

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

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

5　結(jié)語

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

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

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

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

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

［1］林學(xué)海，任勇生.基于ONERA氣動力模型的風(fēng)力機(jī)葉片顫振時(shí)域分析［J］.山東科技大學(xué)學(xué)報(bào)（自然科學(xué)版），2009，28（3）：56-60.

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

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

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

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

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

［7］任勇生，張明輝.水平軸風(fēng)力機(jī)葉片的彎扭耦合氣彈穩(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］.力學(xué)季刊，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ìn)軸系彎—扭耦合振動特性分析［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

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

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

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