大型水平軸風力機葉片二自由度氣彈穩定性分析

朱由鋒，任勇生

(山東科技大學 機械電子工程學院，青島 266510)

目前，對大型風力機葉片的結構模型主要有彈性鉸鏈模型、連續分布參數葉片模型[1-7]。Kalles?e[8]在Hodges-Dowell葉片模型的基礎上，引入新的葉片坐標系統，建立了各向同性材料葉片的機械振動方程，但沒有對整個葉片的氣彈穩定性進行深入研究。本文將Kalles?e模型下的實心各向同性材料的大型水平軸風力機葉片的動力學方程進行了簡化，利用Greenberg氣動力理論，建立了葉片線性定常系統氣彈方程，并對葉片的氣彈穩定性的影響參數進行了分析。

目前對大型水平軸風力機葉片的顫振研究的算法有兩種:模態疊加法和直接積分法[9-10]。對比例阻尼或模態阻尼系統，可以用實模態分析的方法使方程解耦，用模態疊加的方法求得系統的響應。對一般阻尼不存在模態對阻尼的正交性的問題，可以化為狀態方程把方程降為一階，用解廣義特征值和特征向量的方法，組成變換矩陣、即采用復模態分析方法而使方程解耦，通過復模態疊加求得系統的解。本文所采用的狀態空間分析方法無需坐標變換使方程解耦，適用于一般的機械振動系統，也便于在計算機上對系統進行時域和頻域分析。

1 系統方程建立

1.1 運動微分方程

Kalles?e結構模型如圖1所示，它最大的特點就是葉片切面不垂直于葉片的中心軸，這極大的簡化了擺振、揮舞方向無阻尼振動方程，但增加了扭轉方向的復雜性。由于扭轉剛度較擺振、揮舞剛度高，所以可以用二自由度代替之。在公式中保留前二階量，忽略高階小量，可以得到x，y方向無氣動載荷的線性運動方程分別為:

圖1 水平軸風力機模型Fig.1 The model of horizontal axis turbine blade

式中，m為葉片等效線質量密度kg/m;Ω為葉片的旋轉速度(rad/min);R為葉片的半徑(m);E為葉片等效彈性模量(Pa);Ιζ為葉片 ζ向的慣性矩(kg·m2);Ιη為葉片η向的慣性矩(kg·m2);θ為葉片的預扭角(deg);β為葉片的槳矩角(rad)。

1.2 空氣動力方程

利用Greenberg氣動力理論[11]可以得到氣動力的無環量升力、有環量升力和阻力。利用本文介紹算法進行驗證，氣動力導致的慣性和剛度的變化對系統的穩定性幾乎不受影響，可忽略。根據坐標轉換和去高階簡化，可以導出小攻角時風力機葉片x，y方向分布氣動力為:

式中，ρ為空氣密度(kg/m3);α0為升力曲線斜率;c為弦長，(m);νin為風速(m/s)。

2 狀態空間方程的建立

2.1 有限差分離散

在差分離散之前，必須對不同類型風力機葉片的技術參數進行差值求解或者曲線擬合，得到不同節點的數值。式(1)～式(4)中的參數，m=m(s)，θ=θ(s)，Ιζ=Ιζ(s)，Ιη=Ιη(s)，C=C(s)。葉片的邊界條件為:

首先對X，Y方向的四階導數進行差分處理，由于X，Y方向的邊界條件相同，可以X方向的四階導數為例進行研究。根據泰勒的展開形式，運用二階中心差分后，可以得到X方向各次導數的一般差分形式。根據一次差分公式和邊界條件可以得到:

根據二次、三次差分公式和邊界條件:

求解后得到:

最后計算得到四階導數的差分為:

同理，一階導數的差分為:

二階導數的差分為:

三階導數的差分為:

將差分公式代入式(1)～式(4)中的微分項后，對氣動分布力進行合成，最后可得到風機葉片的矩陣表達式:

式中:M為(2N×2N)質量矩陣;C為(2N×2N)阻尼矩陣;K 為(2N ×2N)剛度矩陣;q 為 q=[μ1，μ2，μ3，…，μn，ν1，ν2，ν3，…，νn]。

不考慮葉片的旋轉耦合和氣動力的影響，M為對角實矩陣。由于進行了氣動力的簡化，C為一般粘性阻尼形式。K為不對稱剛度矩陣，即使對q進行矩陣行變換也無法解耦，所以不能用傳統的疊加法進行求解。

2.2 狀態空間轉換

根據前面進行的簡化，風力機葉片成為線性定常系統問題，狀態空間[12]的表達式為:

根據現代控制理論和公式(13)構建系統的狀態空間方程為:

式中:F1為(2N×2N)零矩陣;IE為(2N×2N)單位矩陣;F2為(1×(4N-1))零矩陣;F3為(1×2N))零矩陣;

3 數值計算及分析

3.1 參數計算與選取

為了進行數值計算與分析，選取風力機的典型參數[13]作為葉片系統輸入值。通過單個參數的變化分析系統的穩定性，選取葉片截面橢圓形狀，變量參數計算如下:

式中:ρ為葉片等效密度(kg/m3);b為葉片厚度(b=0.25 c);

常量參數選擇:R=60，E=0.8 ×1010，β =2，ρ=1.237，α0=0.3，νin=8。

3.2 穩定性分析

狀態空間全部極點落在左平面，系統穩定。為方便計算機運算和研究參數變化和穩定性的關系，引用極點實部最大值概念。若極點實部的最大值小于零，則認為系統穩定，反之不穩定。變化曲線為:

圖2 葉片參數對穩定性的影響Fig.2 Stability of variable turbine blade parameters

通過圖2我們可以得到以下驗證:

(1)彈性模量和葉片密度是葉片的結構參數，采用優化的空心加強筋設計和高強度復合材材料可以提高風力機的性能。

(2)空氣密度和風速是風力機的外部環境參數，是不可控制因素，適當控制葉片的轉速和槳矩角的大小，可以進行有效的顫振抑制。

(3)從模型計算參數可知，優化葉片表面形狀和復合材料的鋪層方式也能提高葉片的性能。

(4)以上分析的實用性可以反證模型和計算方法的可行性。

3.3 時域頻域分析

利用MATLAB振動工具箱和構建的狀態空間，對系統進行時域頻域分析，可以得到系統的伯德圖、乃奎斯特圖以及時間響應曲線。根據風力機葉片的差分離散公式，N＞5時才不會使系統失真。由于篇幅關系，只選取了 N=6，風速分別為 νin=8，νin=11.6，νin=14時的階躍時間響應進行驗證分析。

圖3 風機葉片的階躍時間響應Fig.3 Step response of turbine blade

圖3 中依次反映了 νin=8，νin=11.6，νin=14 時系統的收斂、顫振和發散階躍時間響應。ln(1)-ln(6)依次反映了6個節點的擺振，ln(7)-ln(12)依次反映了6個節點的揮舞。從圖3中可以得出:

(1)系統階躍時間響應曲線驗證了風速等穩定性曲線的正確性。也側證了模型和算法的可行性。

(2)收斂后的振幅反映了各個節點的變形量，與無阻尼自由振動系統的主振型相吻合。

(3)系統收斂速度和發散速度也可以在圖中體現。

4 結論

本文研究獲得以下主要結論:

(1)對Kalles?e坐標模型進行了簡化和推廣，根據Greenberg非定常氣動理論，推導出適合模型和差分計算的氣動力矩計算公式。

(2)利用狀態空間法解決了無法解耦的振動系統問題，可以推廣到一般振動系統;提出了極點最大實部的判據方法，適于系統參數的穩定性研究和計算機求解。

(3)用系統參數的穩定性曲線和葉片系統的階躍時間響應對模型和算法進行驗證。

[1]Hodges D H，Dowell E H.Nonlinear equation of motion for the elastic bending and torsion of twisted nonuniform rotor blades[R].NASA Tecnical Note，1974:19-25.

[2]Hodges D H，Robert A.Stability of elastic bending and torsion of uniform cantilever rotor blades in hover with variable structural couping[R].NASA Technical Note，1976:15-21.

[3]Gebhardt C G，Preidikman S.Numerical simulations of the aerodynamic behavior of large horizontal axis wind turbines[J].International Journal of Hydrogen Energy，2010，35:6005-6011.

[4]Maalawi K Y.A direct method for evaluating performance of horizontal axis wind turbines[J].Renewable and Sustainable Energy Reviews，2001，5:175-190.

[5]Maheri A，Noroozi S.WTAB，a computer program for predicting the performance of horizontal axis wind turbines with adaptive blades[J].Renewable Energy，2006，31:1673-11685.

[6]Pesmajoglou S D，Grahamb J M R.Prediction of aerodynamic forces on horizontal axis wind turbines in free yaw and turbulence[J].Journal of Wind Engineering and Industrial Aerodynamics，2000，6:1-14.

[7]Chmaissem W.Numerical study of the boussinesq model of natural convection in an annular space:having a horizontal axis bounded by circular and elliptical isothermal cylinders[J].Applied Thermal Engineering，2002，22:1013-1025.

[8]Kalles?e B S.Equations of motion for a rotor blade，including gravity，pitch action[J].Wind Energy，2007，10:211-216.

[9]Hansen M H.Aeroelastic stability analysis of wind turbines using an eigenvalue approach[J].Wind Energy，2004，7:133-143.

[10]任勇生.水平軸風力機葉片的彎扭耦合氣彈穩定性研究[J].振動與沖擊，2010，29(7):197-201.

[11]李本立，宋憲耕，賀德馨等.風力機結構動力學[M].北京:北京航空航天大學出版社，1999，102-113.

[12]林化方.多度系統振動問題的狀態方程解法[J].廣東機械學院學報，1994，12(2):70-74.

[13]Berry D，Jackson K.Parametric study for large wind turbine blades[R].Sand 2002-2519 Unlimited Release Printed，2002:12-18.