郭小鋒，黃鑫祥，車江寧，張 忍

（中原工學(xué)院 機(jī)電學(xué)院，河南 鄭州450007）

0 引言

為了提高風(fēng)輪對風(fēng)能的吸納能力，降低發(fā)電成本，現(xiàn)代風(fēng)力機(jī)葉片的長度、重量和極限載荷較以往有明顯的增加。早期的2MW Tj?reborg風(fēng)力機(jī)葉片長度僅為29.1 m，重量為6.8 t[1]?，F(xiàn)代的2 MW風(fēng)力機(jī)葉片長度一般在50～60m，重量為10～15 t[2]。風(fēng)輪葉片重量及長度的增加使葉片的氣彈扭轉(zhuǎn)變形分析及氣彈失穩(wěn)問題成為風(fēng)力機(jī)設(shè)計(jì)中的關(guān)鍵。

對于風(fēng)力機(jī)的氣彈動力學(xué)問題，傳統(tǒng)的分析方法包含模態(tài)分析法、有限元法（FEM）和多體動力學(xué)方法[3]～[5]。模態(tài)分析的基本原理是將塔筒和風(fēng)輪葉片的變形描述為若干形函數(shù)的線性組合，然后采用虛功法建立整個風(fēng)力機(jī)的動力學(xué)方程。有限元方法一般采用Timoshenko梁或非線性梁理論將葉片離散成多個單元，根據(jù)結(jié)構(gòu)動力學(xué)理論和哈密頓原理建立葉片的動力學(xué)方程。多體動力學(xué)方法是將風(fēng)力機(jī)看作由若干剛性或柔性部件的組合，這些部件通過力元素（如彈簧）或運(yùn)動學(xué)約束（如關(guān)節(jié)）相互連接，然后采用拉格朗日方程或牛頓-歐拉方程構(gòu)建整機(jī)的動力學(xué)模型。

風(fēng)力機(jī)的來流風(fēng)況具有復(fù)雜多變性，風(fēng)輪葉片在氣動力、彈性力及變速變槳控制的耦合作用下呈現(xiàn)出復(fù)雜的動態(tài)響應(yīng)?，F(xiàn)有的分析方法和工具主要側(cè)重于葉片變形及載荷的研究，針對當(dāng)前風(fēng)力機(jī)葉片大型化、柔性化的發(fā)展趨勢，研究風(fēng)力機(jī)葉片在復(fù)雜工況下氣彈穩(wěn)定性分析的新方法，具有重要的意義。

1 風(fēng)力機(jī)多體動力學(xué)模型

在工作狀態(tài)下，風(fēng)力機(jī)葉片和塔筒存在明顯的彈性變形，所以葉片和塔筒需作為柔性體處理，而機(jī)艙、主軸及輪轂等部件在工作時的彈性變形很小，可以看作為剛體。因此，風(fēng)力機(jī)整機(jī)模型是典型的剛?cè)狁詈隙囿w系統(tǒng)[6]。

1.1 柔性體扭轉(zhuǎn)變形的表達(dá)方法

基于多體系統(tǒng)（MBS）理論，柔性體的葉片被看作是若干個截面的組合，在葉片的某截面A處建立當(dāng)?shù)鼐植繀⒖甲鴺?biāo)系Oc-XcYcZc，當(dāng)葉片發(fā)生彈性變形時，截面A相對于其初始狀態(tài)會發(fā)生三維空間上的彈性移動和扭轉(zhuǎn)變形，在某一時刻其位置和姿態(tài)如截面A＇所示，截面A上的參考坐標(biāo)系Oc-XcYcZc轉(zhuǎn)變?yōu)镺c＇-Xc＇Yc＇Zc＇（圖1）。

圖1 柔性葉片氣彈變形示意圖Fig.1 Generalmodel of a flexible body

截面A＇相對于A的彈性扭轉(zhuǎn)變形可以用三維向量θc表示，并采用模態(tài)疊加法求得。

在葉片的氣彈變形分析中，截面的氣彈扭轉(zhuǎn)變形是通過歐拉變換進(jìn)行表達(dá)的。考慮截面A相對于葉根坐標(biāo)系存在氣動扭角β的變換，因此，將O-XYZ坐標(biāo)系平移至Oc點(diǎn)，繞Z軸旋轉(zhuǎn)β，即得到Oc-XcYcZc坐標(biāo)系。同樣，將Oc-XcYcZc坐標(biāo)系平移至Oc＇點(diǎn)，先繞Zc軸旋轉(zhuǎn)ψ，再繞新的Yc軸旋轉(zhuǎn)θ，最后繞新的Zc軸旋轉(zhuǎn)φ，得到Oc＇-Xc＇Yc＇Zc＇。

綜上所示，ψ+β為截面A相對于葉根坐標(biāo)系的角扭轉(zhuǎn)角位移，θ為截面A在揮舞方向的扭轉(zhuǎn)角，而φ為截面A在擺振方向的扭轉(zhuǎn)角。

截面A＇相對于A的彈性位移變形可以用三維向量uc表示，同樣采用模態(tài)疊加法求得[7]。

因此，截面A在某一時刻下的位置和形態(tài)為

1.2 整機(jī)動力學(xué)模型

本文將風(fēng)力機(jī)模型進(jìn)行合理簡化，只包括葉片、輪轂、塔筒、機(jī)艙等主要部件，而不考慮齒輪箱、發(fā)電機(jī)等部件。偏航法蘭與塔筒頂端固接，用于連接機(jī)艙，雖然機(jī)艙中不創(chuàng)建主機(jī)架、齒輪箱、發(fā)電機(jī)等模型，但是，需要把忽略部件的質(zhì)量和轉(zhuǎn)動慣量算入機(jī)艙中。低速軸與機(jī)艙為轉(zhuǎn)動連接，輪轂與低速軸進(jìn)行固接，葉片與輪轂由變槳法蘭連接。所構(gòu)建的風(fēng)力機(jī)整機(jī)模型如圖2所示。

圖2 風(fēng)力機(jī)整機(jī)動力學(xué)模型Fig.2 Dynamicmodel ofwind turbine

圖2中：塔筒底部定義整機(jī)基礎(chǔ)坐標(biāo)系eb；機(jī)艙質(zhì)心處定義en坐標(biāo)系；主軸質(zhì)心和輪轂中心分別定義了es和eh坐標(biāo)系；在葉片根部定義了ebi坐標(biāo)系。

在完成整機(jī)坐標(biāo)系的構(gòu)建后，風(fēng)力機(jī)多體動力學(xué)系統(tǒng)的自由度參數(shù)可表述為

建立好風(fēng)力機(jī)整機(jī)的剛?cè)狁詈隙囿w動力學(xué)模型后，采用凱恩方法對風(fēng)力機(jī)動力學(xué)方程進(jìn)行求解。

2 風(fēng)力機(jī)多體動力學(xué)仿真方法

本文的研究對象為某內(nèi)陸商用2MW風(fēng)力機(jī)，葉片長度為54m，塔筒高度為85m，額定風(fēng)速為9.5 m/s?；诙囿w動力學(xué)軟件SIMPACK，本文所構(gòu)建的風(fēng)力機(jī)氣-彈-控聯(lián)合仿真模型及求解流程如圖3所示。

圖3 風(fēng)力機(jī)氣彈動力學(xué)模型與分析流程圖Fig.3 Aerodynamicmodel and analysis flow chartof wind turbine

2.1 葉片參數(shù)化表達(dá)

首先采用EXCEL表格對葉片的外形、鋪層結(jié)構(gòu)進(jìn)行參數(shù)化設(shè)計(jì)；其次編寫MATLAB程序讀取處理EXCEL表格中的葉片參數(shù)數(shù)據(jù)，調(diào)用PreComp軟件計(jì)算得到SIMPACK葉片建模中所需葉片各截面的屬性參數(shù)（如截面質(zhì)心、彈性中心、剪切中心、扭轉(zhuǎn)剛度、線密度等）；最后對上一步輸出的數(shù)據(jù)進(jìn)行整理，生成Rotorblade Generation模塊的輸入文件，完成風(fēng)輪葉片參數(shù)化模型的構(gòu)建。

2.2 塔筒部件

本文采用SIMPACK中的Flexible（model）模塊進(jìn)行柔性塔筒建模。在已知塔筒內(nèi)外徑、材料等參數(shù)的情況下，借助ANSYS自身的參數(shù)化語言APDL實(shí)現(xiàn)塔筒有限元模型的構(gòu)建。由于塔筒有限元模型的網(wǎng)格規(guī)模龐大、自由度多，為了節(jié)省SIMPACK多柔體數(shù)值計(jì)算的仿真時間，須要對有限元塔筒模型進(jìn)行模型縮減，即對模型進(jìn)行子結(jié)構(gòu)分析。最后由ANSYS計(jì)算輸出包含有限元塔筒模型節(jié)點(diǎn)、單元、材料等信息的.cdb模型文件以及包含塔筒縮減模型質(zhì)量矩陣和剛度矩陣等信息的.sub文件。Flexible模塊調(diào)用這些文件，最后生成塔筒模型。

2.3 空氣動力學(xué)模塊

為了實(shí)現(xiàn)對風(fēng)力機(jī)在復(fù)雜工況下的氣彈動力學(xué)分析，本文采用SIMPACK與AeroDyn聯(lián)合仿真的方法。由ICwind或Turbsim[8]軟件生成相應(yīng)的風(fēng)模型文件，SIMPACK多體動力學(xué)求解模塊通過241號力元接口調(diào)用AeroDyn計(jì)算當(dāng)前時刻的氣動數(shù)據(jù)，根據(jù)氣動數(shù)據(jù)對風(fēng)力機(jī)進(jìn)行動力學(xué)方程求解，計(jì)算得到各部件的變形量和速度變化。

2.4 控制模塊

風(fēng)電機(jī)組采用變速變槳控制策略，即低于額定風(fēng)速時采用轉(zhuǎn)矩控制方式，以發(fā)電機(jī)轉(zhuǎn)矩需求為輸入，以發(fā)電機(jī)測量轉(zhuǎn)速為輸出建立風(fēng)力發(fā)電機(jī)動態(tài)數(shù)學(xué)模型。高于額定風(fēng)速時采用槳距控制方式，以葉片槳距需求為輸入，以發(fā)電機(jī)測量轉(zhuǎn)速為輸出建立風(fēng)力發(fā)電機(jī)動態(tài)數(shù)學(xué)模型。本文采用Fortran基于PI控制原理編寫風(fēng)力機(jī)轉(zhuǎn)速及槳距角控制的動態(tài)鏈接庫文件[9]。

3 結(jié)果對比分析

平均風(fēng)速為額定風(fēng)速的湍流風(fēng)況是風(fēng)力機(jī)極為惡劣的工況，在這種工況下風(fēng)力機(jī)葉片的葉根載荷和葉尖變形量最大，葉片也容易產(chǎn)生氣彈穩(wěn)定性[7]～[9]。因此本文選取平均風(fēng)速為9.5 m/s的湍流風(fēng)對葉片進(jìn)行動態(tài)載荷分析，湍流風(fēng)時間為600 s，其湍流強(qiáng)度縱向?yàn)?0.61%、橫向?yàn)?7.56 %、垂直方向?yàn)?1.04 %。SIMPACK的仿真時間設(shè)為100 s，選擇SODASRT_2求解器對風(fēng)力機(jī)模型進(jìn)行離線積分求解，求解步長設(shè)為0.01 s。選取50～100 s時間段進(jìn)行數(shù)據(jù)對比分析，該時間段內(nèi)輪轂高度處的風(fēng)速變化如圖4所示。

圖4 輪轂高度處的風(fēng)速變化Fig.4Wind speed variation athub height

風(fēng)力機(jī)開始運(yùn)行時處于不穩(wěn)定階段，所以選擇圖4時間段內(nèi)的分析結(jié)果與行業(yè)普遍認(rèn)可的GH-Bladed軟件的分析結(jié)果進(jìn)行對比分析，以驗(yàn)證本文方法的可靠性。

在揮舞、擺振及扭轉(zhuǎn)方向，2MW葉片葉尖截面（距葉根53m處）變形的分析結(jié)果對比見圖5。

圖5 湍流風(fēng)下葉尖動態(tài)氣彈變形Fig.5 Dynamic aerodynamic deformation of blade tip under turbulentwind

由圖5（a）可知：SIMPACK計(jì)算的葉尖平均揮舞位移為5.1 m，最大揮舞位移為7.4 m，在該時間段內(nèi)，葉尖揮舞位移的變化趨勢與圖4中X方向的風(fēng)速變化趨勢一致，說明葉尖揮舞變形主要受X方向的風(fēng)速影響；GH-Bladed計(jì)算的葉片平均揮舞位移為5.3 m，最大揮舞位移為5.8 m。兩種軟件計(jì)算的揮舞位移變化趨勢一致，且相差不大。

由圖5（b）可知，葉尖擺振變形范圍在-1～0.8 m，兩種軟件計(jì)算的擺振位移很接近。

圖5（c）為葉尖在扭轉(zhuǎn)方向的變形量，其中葉片扭轉(zhuǎn)方向的變形量是葉片在工作中各截面變形后相對于葉根坐標(biāo)系歐拉變換角φ，SIMPACK的風(fēng)機(jī)模塊可以直接輸出這個計(jì)算結(jié)果，從而可以非常直接地評估葉片各截面的氣彈扭轉(zhuǎn)變形。通過分析發(fā)現(xiàn)，在風(fēng)力機(jī)運(yùn)行穩(wěn)定后，葉片截面扭轉(zhuǎn)變形通常維持在某一角度，并隨著風(fēng)速、風(fēng)向的變化上下波動。由圖5（c）可知：葉尖截面的平均扭轉(zhuǎn)角度為3.5 °，且在50～100 s時間段內(nèi)扭轉(zhuǎn)角度波動幅度為0～7°；而GH-Bladed計(jì)算的平均扭轉(zhuǎn)角度為3.3 °，在該時間段內(nèi)扭轉(zhuǎn)角度的波動幅度也為0～7°。

4 氣彈穩(wěn)定性分析

在實(shí)際工作中，大型風(fēng)力機(jī)葉片截面會產(chǎn)生較大的彎扭彈性變形，特別是葉尖部位的截面。葉片截面的彈性變形會引起截面翼型攻角的變化，翼型攻角的變化使得葉片截面的氣彈載荷發(fā)生改變，進(jìn)而對其彈性變形產(chǎn)生影響。在極端工況下，這種氣彈耦合響應(yīng)會使葉尖部位產(chǎn)生失速或氣彈失穩(wěn)現(xiàn)象，因此對葉尖翼型的氣彈穩(wěn)定性分析對大型風(fēng)電葉片的設(shè)計(jì)具有重要的意義。

葉尖失速是葉片失穩(wěn)的重要特征之一。對于葉尖翼型，首先通過空氣動力學(xué)軟件AeroDyn分析葉片在旋轉(zhuǎn)過程中的攻角變化情況，并計(jì)算攻角對應(yīng)的升力系數(shù)；最后，定義1 s為一個周期來研究葉尖翼型的攻角-升力系數(shù)變化情況，在一個周期內(nèi)，同一攻角可對應(yīng)多個升力系數(shù)值，取這些升力系數(shù)中的最大值和最小值，這兩個值之間的差值越大，說明葉片在這個周期內(nèi)經(jīng)歷了較大的氣動載荷變化。

對于本文研究的2MW風(fēng)力機(jī)葉片，在前文所述的湍流風(fēng)況下，葉尖截面翼型（距葉根53m處）在控制系統(tǒng)正常和控制系統(tǒng)失效兩種工況下的攻角-升力系數(shù)關(guān)系如圖6所示。

圖6 葉尖截面翼型的攻角—升力系數(shù)曲線圖Fig.6 Angle of attack—lift coefficient curve of blade tip section

由圖6可知：在控制系統(tǒng)正常的情況下，葉尖截面在同一攻角下對應(yīng)的升力系數(shù)比較接近，升力系數(shù)之差的最大值為0.046；而在控制系統(tǒng)失效的情況下，最大升力系數(shù)與最小升力系數(shù)相差0.2，此時葉尖截面呈現(xiàn)較典型的失速狀態(tài)[10]，[11]。

為了進(jìn)一步分析葉片的氣彈穩(wěn)定性，在兩種工況下，對葉尖在擺振和扭轉(zhuǎn)方向的變形進(jìn)行頻域分析，分析結(jié)果如圖7所示。

圖7 葉尖變形的頻域分析Fig.7 Frequency domain analysis of blade tip deformation of the pre-bend blades

葉片一階擺振的固有頻率為0.962 Hz，在控制系統(tǒng)失效的情況下，圖7中葉片在0.962 Hz頻率處產(chǎn)生較大的扭轉(zhuǎn)和擺振變形，表明葉片在該工況下產(chǎn)生了明顯的氣彈失穩(wěn)現(xiàn)象[8]，[12]。

5 結(jié)論

本文以某2MW風(fēng)力機(jī)為研究對象，基于SIMPACK軟件了構(gòu)建了整機(jī)剛?cè)狁詈系臍?彈-控聯(lián)合仿真模型和方法，通過與行業(yè)內(nèi)公認(rèn)的風(fēng)力機(jī)動力學(xué)分析軟件GH-Blaed相比，本文提出的風(fēng)力機(jī)氣彈性能分析方法和程序可準(zhǔn)確分析風(fēng)輪葉片的載荷和氣彈變形量。相比于GH-Bladed軟件，本文采用與AeroDyn聯(lián)合氣動仿真的方法，可直觀地分析葉片在一個周期內(nèi)的升力系數(shù)變化情況，進(jìn)而實(shí)現(xiàn)對葉片葉尖截面失速失穩(wěn)現(xiàn)象的預(yù)測分析。