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

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

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

0 引言

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

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

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

1 風力機多體動力學模型

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

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

基于多體系統（MBS）理論，柔性體的葉片被看作是若干個截面的組合，在葉片的某截面A處建立當地局部參考坐標系Oc-XcYcZc，當葉片發生彈性變形時，截面A相對于其初始狀態會發生三維空間上的彈性移動和扭轉變形，在某一時刻其位置和姿態如截面A＇所示，截面A上的參考坐標系Oc-XcYcZc轉變為Oc＇-Xc＇Yc＇Zc＇（圖1）。

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

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

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

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

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

因此，截面A在某一時刻下的位置和形態為

1.2 整機動力學模型

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

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

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

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

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

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

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

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

2.1 葉片參數化表達

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

2.2 塔筒部件

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

2.3 空氣動力學模塊

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

2.4 控制模塊

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

3 結果對比分析

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

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

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

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

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

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

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

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

4 氣彈穩定性分析

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

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

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

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

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

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

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

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

5 結論

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