曲珍壯，夏鴻建，李德源，林偉豪

（廣東工業(yè)大學(xué)機(jī)電工程學(xué)院，廣東 廣州 510006）

0 引言

為了增強(qiáng)風(fēng)力機(jī)的風(fēng)能捕獲能力，降低風(fēng)力發(fā)電成本，風(fēng)力機(jī)葉片的尺寸正逐漸變大，致使風(fēng)力機(jī)葉片柔性增強(qiáng)，在慣性力、氣動(dòng)力的綜合影響下，大型風(fēng)力機(jī)葉片的彈性大變形和氣動(dòng)載荷耦合特性增強(qiáng)，氣彈失穩(wěn)導(dǎo)致的顫振問題突出。

近年來，對(duì)于大型風(fēng)力機(jī)葉片的顫振特性分析，國(guó)內(nèi)外學(xué)者進(jìn)行了大量的研究,主要方法可分為時(shí)域分析方法和頻域分析方法。時(shí)域分析方法通過判斷振動(dòng)位移的發(fā)散趨勢(shì)，確定葉片顫振邊界的風(fēng)速和轉(zhuǎn)速。如王曉東[1]采用FAST 軟件，在時(shí)域條件下對(duì)NREL 5 MW 風(fēng)力機(jī)葉片在運(yùn)轉(zhuǎn)和靜止時(shí)的顫振性能進(jìn)行了分析預(yù)測(cè)，并研究了扭轉(zhuǎn)剛度、質(zhì)心位置、結(jié)構(gòu)阻尼和槳距角對(duì)顫振臨界風(fēng)速的影響。時(shí)域分析方法雖然可以快速判斷顫振邊界，但難以獲得引起葉片顫振的結(jié)構(gòu)振動(dòng)模態(tài)和氣動(dòng)阻尼。頻域分析方法，通過對(duì)氣彈耦合方程線性化，并采用復(fù)特征值分析方法，判別顫振邊界，能獲得引起顫振的結(jié)構(gòu)模態(tài)和氣動(dòng)阻尼變化。如HANSEN[2]采用（Multi-Blade Coordinates，MBC）坐標(biāo)和模態(tài)法，結(jié)合葉素動(dòng)量理論和B-L 模型，分析了NREL 5 MW 風(fēng)力機(jī)顫振臨界轉(zhuǎn)速；李德源[3]采用多體動(dòng)力學(xué)方法和葉素動(dòng)量理論的方法建立葉片氣彈耦合方程，研究了NREL 5 MW 風(fēng)力機(jī)氣彈穩(wěn)定性問題；ZHOU[4]采用混合變分幾何精確梁模型和葉素動(dòng)量理論，分析了5 MW 風(fēng)力機(jī)葉片的彎扭耦合顫振特性。

然而，為獲得相對(duì)精確的氣彈顫振邊界，分析引起顫振的結(jié)構(gòu)模態(tài)耦合特點(diǎn)和氣彈耦合特性，需能高精度描述葉片大變形特點(diǎn)的非線性結(jié)構(gòu)模型，以及結(jié)構(gòu)振動(dòng)與氣動(dòng)耦合的非定常氣動(dòng)模型。傳統(tǒng)采用歐拉梁結(jié)構(gòu)模型[5-6]或模態(tài)法[7]構(gòu)建葉片結(jié)構(gòu)模型，難以表達(dá)柔性葉片的幾何非線性變形，以及轉(zhuǎn)動(dòng)產(chǎn)生的動(dòng)力剛化現(xiàn)象。對(duì)轉(zhuǎn)動(dòng)葉片非線性結(jié)構(gòu)建模，主要有幾何精確梁模型[4,8]、絕對(duì)節(jié)點(diǎn)坐標(biāo)模型和超級(jí)單元模型[3,9]等。其中，幾何精確梁方法由于具有梁模型的自由度少、精度高，并且能準(zhǔn)確表示非線性大變形的優(yōu)點(diǎn)，在直升機(jī)葉片和風(fēng)力機(jī)葉片獲得廣泛應(yīng)用。然而，傳統(tǒng)幾何精確梁模型由于存在截面姿態(tài)角度插值和動(dòng)力學(xué)線性化相對(duì)繁瑣等問題，多數(shù)應(yīng)用于大型葉片的時(shí)域響應(yīng)分析。另一方面，大型柔性葉片振動(dòng)會(huì)影響葉片氣動(dòng)載荷，氣動(dòng)載荷又反饋?zhàn)饔糜诮Y(jié)構(gòu)，從而形成氣彈耦合特性。傳統(tǒng)的準(zhǔn)靜態(tài)的葉素動(dòng)量理論采用靜態(tài)氣動(dòng)系數(shù)，計(jì)算氣動(dòng)載荷，忽略了葉片振動(dòng)引起的葉片氣動(dòng)非定常特性。為有效描述葉片的非定常特性，學(xué)者們分別提出了Boeing-Vertol 模型，ONERA 模型，?ye模型以及半經(jīng)驗(yàn)B-L（Beddoes-Leishman）[10]等動(dòng)態(tài)失速模型。由于B-L 模型具有求解方便，并且容易線性化的特點(diǎn)，獲得廣泛應(yīng)用。

為此，本文針對(duì)大型風(fēng)力機(jī)葉片結(jié)構(gòu)非線性以及氣動(dòng)載荷的非定常特點(diǎn)，采用幾何精確梁模型與集成B-L（Beddoes-Leishman）動(dòng)態(tài)失速模型的葉素動(dòng)量方法，建立風(fēng)力機(jī)葉片的氣彈耦合動(dòng)力學(xué)方程。為進(jìn)一步提高氣彈模型的求解精度，便于對(duì)葉片結(jié)構(gòu)方程線性化，由弱形式的求積元方法，高效構(gòu)建葉片的氣彈耦合動(dòng)力學(xué)離散平衡方程；再基于變分原理，綜合線性化氣彈耦合動(dòng)力學(xué)方程，由復(fù)特征值分析法，計(jì)算葉片的氣彈頻率和氣彈阻尼比，辨識(shí)葉片氣彈顫振區(qū)域；并以DTU-10MW 風(fēng)力機(jī)葉片，研究大型風(fēng)力機(jī)葉片的氣彈顫振特性。

1 幾何精確梁變形的描述

幾何精確梁采用平截面假設(shè)，梁的位形由梁軸線的位移和相應(yīng)截面的轉(zhuǎn)動(dòng)決定。如圖1 有全局坐標(biāo)系（不考慮葉片預(yù)彎預(yù)扭，以及截面偏心等因素），不變形位置坐標(biāo)系（考慮葉片真實(shí)幾何形狀）和變形位置坐標(biāo)系（考慮氣彈耦合后葉片的變形）。

圖1 梁的轉(zhuǎn)動(dòng)坐標(biāo)系

在變形位置，梁上任意一點(diǎn)的位置向量可以寫為：

式中：

u——剪切中心在全局坐標(biāo)系下的坐標(biāo)；

z——截面在梁軸線上的坐標(biāo)，m；

x——截面在x軸的坐標(biāo)，m；

y——截面在y軸的坐標(biāo)，m；

Λ——截面轉(zhuǎn)動(dòng)矩陣；

——截面坐標(biāo)系下的位置向量。

三維旋轉(zhuǎn)矩陣Λ有9 個(gè)元素，包含6 個(gè)約束關(guān)系，難以直接作為姿態(tài)變量，歐拉四元數(shù)是描述截面轉(zhuǎn)動(dòng)非奇異的最小參數(shù)表達(dá)，本文采用歐拉四元數(shù)描述截面的轉(zhuǎn)動(dòng)[11]。

式中：

q0——四元數(shù)的標(biāo)量部分；

q——四元數(shù)的向量部分。

四元數(shù)和轉(zhuǎn)動(dòng)矩陣Λ的關(guān)系為：

式中：

I——單位矩陣。

2 結(jié)構(gòu)動(dòng)力學(xué)方程的建立

2.1 幾何精確梁本構(gòu)關(guān)系

幾何精確梁理論基于Reissner 應(yīng)變，根據(jù)參考文獻(xiàn)[12-13]，定義Reissner 應(yīng)變?chǔ)槠絼?dòng)應(yīng)變向量，應(yīng)變K為轉(zhuǎn)動(dòng)應(yīng)變向量，截面等效應(yīng)力N，M和應(yīng)變?chǔ)?，K之間的關(guān)系為：

式中：

S——截面剛度矩陣。

1.3 樣本采集及病原菌鑒定 觀察病害癥狀，用數(shù)碼相機(jī)拍照，填寫樣方調(diào)查表，采集典型樣本壓制成標(biāo)本，帶回實(shí)驗(yàn)室采用科赫法則對(duì)病原菌進(jìn)行驗(yàn)證和鑒定。對(duì)白粉病、銹病等病癥明顯的病害，直接標(biāo)本制片或以徒手切片鏡檢，參照相關(guān)資料[10-19]進(jìn)行鑒定；對(duì)病癥不明顯的病害，進(jìn)行病原菌常規(guī)組織分離[20]，獲得純培養(yǎng)的病原菌，再進(jìn)一步鑒定并做致病性測(cè)定。

風(fēng)力機(jī)葉片由于采用復(fù)合材料鋪層和復(fù)雜的截面形狀，截面耦合特性明顯。本文采用專業(yè)的截面分析軟件VABS 計(jì)算截面剛度矩陣S。

應(yīng)變與變形的關(guān)系：

式中：

Γ0——梁的初始位置；

——梁的初始曲率參數(shù)；

（）'——在全局坐標(biāo)中對(duì)截面在Z 軸上坐標(biāo)的導(dǎo)數(shù)；

（）*——四元數(shù)的共軛；

‘°’——四元數(shù)乘法。

2.2 動(dòng)力學(xué)平衡方程

首先，建立風(fēng)力機(jī)葉片弱形式動(dòng)力學(xué)平衡方程，然后通過數(shù)值積分和微分方法離散求解。

對(duì)內(nèi)力沿單元軸線積分，得到內(nèi)力虛功為：

式中：

δ——變分符號(hào)；

e——葉片單元；

le——葉片單元長(zhǎng)度。

對(duì)氣動(dòng)力沿軸線積分，重力進(jìn)行體積積分，得到單元的外力虛功：

式中：

δua——軸線氣動(dòng)節(jié)點(diǎn)的虛位移；

δθa——軸線節(jié)點(diǎn)氣動(dòng)點(diǎn)的虛轉(zhuǎn)動(dòng)；

fa——軸線節(jié)點(diǎn)的氣動(dòng)力向量；

ma——軸線節(jié)點(diǎn)的氣動(dòng)彎矩向量；

g——重力加速度向量；

ρ——葉片體密度，kg/m3。

各點(diǎn)的速度進(jìn)行體積積分，得到單元的慣性力虛功：

式中：

——位置向量對(duì)時(shí)間求二階導(dǎo)。

根據(jù)虛功原理，單元內(nèi)力，外力和慣性力虛功滿足：

為求解平衡方程，需要通過數(shù)值積分與數(shù)值微分方法，轉(zhuǎn)換為節(jié)點(diǎn)求和形式，本文采用Gauss-Lobatto 積分[15]進(jìn)行數(shù)值離散，Gauss-Lobatto 積分一般形式為：

式中：

積分區(qū)域需轉(zhuǎn)換到標(biāo)準(zhǔn)區(qū)域[-1，1]；

ξ——標(biāo)準(zhǔn)化坐標(biāo)；

f（ξ）——被積函數(shù)；

n——積分點(diǎn)個(gè)數(shù)；

ξi——積分點(diǎn)；

ωi——對(duì)應(yīng)的積分權(quán)系數(shù)。平衡方程中還存在對(duì)葉片軸線坐標(biāo)z 的導(dǎo)數(shù)項(xiàng)，需要采用微分求積法[15]，轉(zhuǎn)換為單元節(jié)點(diǎn)加權(quán)求和形式，微分求積的一般形式為：

式中：

f（ξj）——函數(shù)在離散點(diǎn)上的值；

——微分求積權(quán)系數(shù)。

通過Gauss-Lobatto 積分和微分求積，得到基于節(jié)點(diǎn)虛位移表達(dá)的內(nèi)力，外力和慣性力虛功[8]：

節(jié)點(diǎn)k虛位移表示為：

葉片的一個(gè)單元位移δde=［δd1δd2···δdK］。本文風(fēng)力機(jī)葉片每個(gè)單元包含21 個(gè)節(jié)點(diǎn)。

把公式（13）代入公式（10），并在全局坐標(biāo)系下進(jìn)行組裝，得到等效節(jié)點(diǎn)內(nèi)力Rint、外力Rext和慣性力Rinert的離散平衡方程：

3 動(dòng)力學(xué)方程線性化

為了分析葉片的振動(dòng)頻率和模態(tài)特性，需要對(duì)動(dòng)力學(xué)平衡方程進(jìn)行線性化處理，分為2 個(gè)部分：結(jié)構(gòu)方程線性化和氣動(dòng)力線性化。

3.1 結(jié)構(gòu)動(dòng)力學(xué)線性化

分別對(duì)節(jié)點(diǎn)內(nèi)力，外力和慣性力矩陣線性化得到[8]：

將公式（16）-（18）代入到平衡方程（15）可以得到：

式中：

K=Kint-Kext+Kinert，M=Minert，C=Cinert；

μ——弱形式求積長(zhǎng)度系數(shù)；

B——微分求積定位矩陣。

3.2 氣動(dòng)力部分線性化

根據(jù)葉素動(dòng)量定理，作用在葉片第false個(gè)截面的氣動(dòng)力和氣動(dòng)力矩可以寫為：

式中：

c——截面弦長(zhǎng)，m；

ρa(bǔ)——空氣密度，kg/m3；

W——相對(duì)來流風(fēng)速，m/s。

C dyn=（CL，CD，CM）T；其 中CL、CD、CM分 別是升力系數(shù)，阻力系數(shù)和扭矩系數(shù)，葉片翼型截面如圖2 所示。

圖2 葉片翼型截面

考慮葉片振動(dòng)和入流角的影響，相對(duì)來流風(fēng)速W和攻角α可以寫為：

式中：

Vn——垂直旋轉(zhuǎn)平面的來流風(fēng)速，m/s；

——擺振方向振動(dòng)速度，m/s；

——揮舞方向振動(dòng)速度，m/s；

R——翼型截面剪心到輪轂中心的距離，m；

Ω——葉片旋轉(zhuǎn)角速度，rad/s。

β——槳距角，（°）；

θ——葉片扭角，（°）。

對(duì)氣動(dòng)力求變分

式中：

系數(shù)矩陣Tx，Tc，Tα，Tα.的具體表達(dá)可以參考文獻(xiàn)[5,16]；

——攻角對(duì)時(shí)間的導(dǎo)數(shù)；

α3/4——葉片翼型截面3/4 弦線處的攻角，求變分得：

式中：

L0——翼型截面3/4 弦點(diǎn)到彈性中心的距離。

為了考慮渦旋脫落和尾緣流動(dòng)分離對(duì)氣動(dòng)特性的影響，引入基于狀態(tài)變量xa=｛x1x2x3x4｝T的B-L動(dòng)態(tài)失速模型[16]，修正葉片各翼型截面處的升力系數(shù)及其斜率。

對(duì)狀態(tài)變量線性化得到

式中：

Ms，Ks，Cs——狀態(tài)變量線性化引入的質(zhì)量陣，剛度陣和阻尼陣；

Ad——系數(shù)矩陣可以參考文獻(xiàn)[16]。

將式（24）-（26）代入（23）中得到氣動(dòng)載荷的線性化表示：

Ka，Ma，Ca——?dú)鈩?dòng)部分線性化得到的氣動(dòng)質(zhì)量陣，剛度陣和阻尼陣。

聯(lián)立式（19），（28）和（29），氣彈線性平衡方程可以寫為：

4 特征值復(fù)模態(tài)分析

根據(jù)式（30）建立風(fēng)力機(jī)顫振分析的特征值平衡方程：

利用Matlab 軟件提供的eig函數(shù)計(jì)算特征值和特征向量，得到葉片的氣彈頻率、氣彈阻尼比和振型。其中特征值的虛部為氣彈頻率，特征值的實(shí)部為氣彈阻尼比。

5 計(jì)算結(jié)果

5.1 結(jié)構(gòu)頻率驗(yàn)證

不考慮風(fēng)力機(jī)葉片結(jié)構(gòu)阻尼和氣動(dòng)力作用的特征方程為：

式中：

K——公式（19）中的結(jié)構(gòu)質(zhì)量矩陣；

M——公式（19）中的剛度矩陣。

通過求解特征方程可以獲得葉片的固有頻率ωn和模態(tài)。

為了驗(yàn)證所建立的葉片結(jié)構(gòu)的動(dòng)態(tài)精度，以DTU-10 MW 風(fēng)力機(jī)葉片為研究對(duì)象，計(jì)算了葉片的固有頻率，并與文獻(xiàn)[17] 中發(fā)表的結(jié)果進(jìn)行了比較。對(duì)比結(jié)果如表1 所示，最大誤差為3.90%，說明本文建立的葉片模型是正確的。圖3 還繪制了第一階揮舞振型、第一階擺振振型和第一階扭轉(zhuǎn)振型模態(tài)圖。

表1 風(fēng)力機(jī)葉片固有頻率對(duì)比

圖3 葉片的模態(tài)圖

5.2 運(yùn)行工況下氣彈頻率和氣彈阻尼比

根據(jù)參考文獻(xiàn)[17]提供的葉片數(shù)據(jù)和運(yùn)行工況參數(shù)，計(jì)算了DTU-10MW 風(fēng)力機(jī)葉片在運(yùn)行工況下的低階氣彈頻率和氣彈阻尼比，并與HANSEN[18]的結(jié)果進(jìn)行對(duì)比。從圖4 中可以看出，在運(yùn)行工況下，一階擺振頻率基本保持不變，和HANSEN 的結(jié)果也比較吻合。

圖4 運(yùn)行工況下一階擺振氣彈頻率和阻尼比

一階擺振氣彈阻尼比和HANSEN 相比，在風(fēng)速4～11 m/s 區(qū)間趨勢(shì)基本相同，都處于緩慢上升趨勢(shì)，風(fēng)速11 m/s 之后本文的結(jié)果在略微下降后有比較明顯的上升趨勢(shì)，由于和HANSEN 采用的模型不同，以及截面耦合矩陣S的計(jì)算軟件的不同，所以在結(jié)果上有一定差異。從圖5可以看出，在運(yùn)行工況下，一階揮舞氣彈阻尼，在風(fēng)速4～7 m/s的保持穩(wěn)定，在風(fēng)速7～12 m/s的處于上升趨勢(shì)，之后趨勢(shì)保持平穩(wěn)，此時(shí)葉片達(dá)到額定轉(zhuǎn)速9.6 r/min，結(jié)果和HANSEN的較為接近，揮舞方向氣彈頻率在到達(dá)額定轉(zhuǎn)速后，也保持穩(wěn)定。

圖5 運(yùn)行工況下一階揮舞氣彈頻率和阻尼比

5.3 定葉尖速比工況下顫振特性分析

為了研究風(fēng)力機(jī)葉片的顫振極限，對(duì)風(fēng)力機(jī)葉片在大變形和高轉(zhuǎn)速顫振特性進(jìn)行分析。將風(fēng)力機(jī)葉片運(yùn)行工況設(shè)為定葉尖速比10，即風(fēng)力機(jī)葉尖速度為來流風(fēng)速的10 倍，并將葉片的槳距角設(shè)為0。分析結(jié)果如圖6，在葉片轉(zhuǎn)速0.6～2.2 rad/s 之間，揮舞方向氣彈頻率和和擺振方向氣彈頻率都在增大，但揮舞方向變化率明顯比擺振方向大，分別為106%和9%，這種現(xiàn)象也和之前對(duì)5 MW 風(fēng)力機(jī)顫振特性的研究結(jié)果相似[3-4]。揮舞方向氣彈阻尼比為較大的正值，擺振方向氣彈阻尼比下降到較小的正值后保持穩(wěn)定，扭轉(zhuǎn)方向氣彈頻率先緩慢上升，然后呈下降趨勢(shì)，扭轉(zhuǎn)方向氣彈阻尼比在轉(zhuǎn)速1.77 rad/s 時(shí)出現(xiàn)負(fù)值，對(duì)應(yīng)的氣彈頻率為8 Hz，且隨轉(zhuǎn)速增加持續(xù)下降，說明此時(shí)風(fēng)力機(jī)葉片發(fā)生顫振失穩(wěn)現(xiàn)象。

圖6 定葉尖速比工況下氣彈頻率和氣彈阻尼比

6 結(jié)論

本文采用非線性幾何精確梁模型和葉素動(dòng)量定理建立風(fēng)力機(jī)葉片的氣彈模型，在B-L（Beddoes-Leishman）動(dòng)態(tài)失速模型的基礎(chǔ)上對(duì)氣彈模型線性化，構(gòu)建特征平衡方程，通過特征值分析，計(jì)算風(fēng)力機(jī)葉片的氣彈頻率和氣彈阻尼比，主要結(jié)論如下。

在運(yùn)行工況下，一階擺振氣彈頻率保持穩(wěn)定,受葉片模型和耦合等因素影響較小。一階擺振氣彈阻尼比變化趨勢(shì)為：在額定轉(zhuǎn)速位置上升到最高值，經(jīng)歷短暫的下降后繼續(xù)上升。本文的結(jié)果和HANSEN相比，上升趨勢(shì)明顯，表明一階擺振氣彈阻尼比受模型和截面耦合影響較大。一階揮舞氣彈頻率和氣彈阻尼比都是在額定轉(zhuǎn)速位置上升到最高值，然后保持穩(wěn)定趨勢(shì)。

在定葉尖速比工況下，隨著轉(zhuǎn)速的增加，風(fēng)力機(jī)葉片揮舞方向氣彈頻率快速增加，擺振方向氣彈頻率變化較小，擺振方向氣彈頻率緩慢下降到較小的正值后保持平穩(wěn)，沒有出現(xiàn)負(fù)值。一階扭轉(zhuǎn)方向氣彈頻率在轉(zhuǎn)速1.77 rad/s 時(shí)，出現(xiàn)負(fù)值，并且繼續(xù)下降，風(fēng)力機(jī)葉片出現(xiàn)顫振失穩(wěn)現(xiàn)象。

本文建立的氣彈模型，還可以用來分析風(fēng)力機(jī)葉片其它的氣彈問題。