基于整機剛柔耦合的風力機葉片氣彈穩(wěn)定性分析

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

（中原工學院 機電學院，河南 鄭州450007）

0 引言

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

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

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

1 風力機多體動力學模型

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

1.1 柔性體扭轉變形的表達方法

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

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

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

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

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

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

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

1.2 整機動力學模型

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

圖2 風力機整機動力學模型Fig.2 Dynamicmodel ofwind turbine

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

在完成整機坐標系的構建后，風力機多體動力學系統(tǒng)的自由度參數(shù)可表述為

建立好風力機整機的剛柔耦合多體動力學模型后，采用凱恩方法對風力機動力學方程進行求解。

2 風力機多體動力學仿真方法

本文的研究對象為某內(nèi)陸商用2MW風力機，葉片長度為54m，塔筒高度為85m，額定風速為9.5 m/s。基于多體動力學軟件SIMPACK，本文所構建的風力機氣-彈-控聯(lián)合仿真模型及求解流程如圖3所示。

圖3 風力機氣彈動力學模型與分析流程圖Fig.3 Aerodynamicmodel and analysis flow chartof wind turbine

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

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

2.2 塔筒部件

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

2.3 空氣動力學模塊

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

2.4 控制模塊

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

3 結果對比分析

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

圖4 輪轂高度處的風速變化Fig.4Wind speed variation athub height

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

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

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

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

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

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

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

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

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

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

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

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

為了進一步分析葉片的氣彈穩(wěn)定性，在兩種工況下，對葉尖在擺振和扭轉方向的變形進行頻域分析，分析結果如圖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)生較大的扭轉和擺振變形，表明葉片在該工況下產(chǎn)生了明顯的氣彈失穩(wěn)現(xiàn)象[8]，[12]。

5 結論

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