旋轉風力機葉片氣動彈性偏移變形的數值模擬研究

王旭東 陳 進 王立存 陸群峰

1．重慶工商大學廢油資源化技術與裝備教育部工程研究中心，重慶，400067 2．重慶大學機械傳動國家重點實驗室，重慶，400030

0 引言

風力機系統作為一個多自由度系統，在對其進行動力學特性研究時就必須研究系統的穩定性。在進行風力機系統設計時，為了保證這個多自由度系統的穩定性，首先需要對每個部件的載荷和穩定性進行研究［1－2］。

風力機葉片作為風力機的關鍵部件，其良好的設計、可靠的質量和優越的性能是保證機組正常穩定運行的決定因素。隨著風力機向大功率發展，葉片的載荷及振動偏移特性的研究就更突出其重要性，而當風力機塔架的剛度很大時，風力機系統的動力學特性主要取決于風力機葉片動力學特性，可以說，葉片作為風力機的主要部件，其振動偏移和穩定性研究直接決定了風力機系統設計的可靠性［3－5］。因此，對風力機葉片的振動偏移和穩定性進行研究就顯得尤為必要。本文從結構動力學理論出發，基于風力機旋轉風輪的氣動彈性模型，提出了旋轉葉片偏移變形特性的數值模擬計算理論和方法。

1 旋轉風輪的載荷特性模型

作用在風輪上的空氣動力學載荷包括各個葉片拍打方向所承受的推力和擺振方向所承受的剪力，推力和剪力的計算采用動量葉素理論，利用Prandtl葉尖損失修正因子和Shen葉尖損失修正模型，對計算過程中葉片每個葉素中的法向力和切向力系數 （Cn，Ct）進行修正［6］。那么，葉片上推力和剪力的計算公式就可以表示為

它所產生的力矩的計算公式為

式中，B為葉片個數；λ為葉尖速比；c為葉片弦長。

當風輪旋轉時，將會對葉片產生離心力，并且方向沿葉片向外。同時，在風輪運轉的過程中，作用在葉片上的空氣動力學載荷所產生的推力使得葉片偏離旋轉平面，進而使得作用在葉片上的離心力方向發生了改變，如圖1所示。此時葉片所承受的離心力可以表示為

式中，R為風輪半徑；Ω為風輪轉速；mi為風輪葉素質量。

圖1 葉片的離心力示意圖

從圖1可以看出，葉片離心力在拍打方向的分量為

式中，δ為葉片偏離旋轉平面的角度；θ為風輪的錐角。

可以看出，由于葉片離心力在拍打方向的分量Fcn與葉片所承受的推力Tx方向相反，因此，該離心力的存在能夠減小拍打彎矩帶來的葉片偏離距離。

2 風輪葉片偏移變形的動力學模型

根據虛功原理，基于結構動力學理論建立葉片的動力學方程：

式中，M為質量矩陣；C為阻尼矩陣；K為剛度矩陣；Fg為外載荷，即空氣動力學載荷和離心力；x為位移向量。考慮葉片前三階振動模態的影響，取11個自由度，即

式中，w為風輪彈性位移；φ為葉片的角度位置；xji為葉片i在第j階模態時所產生的位移；i，j＝1，2；上標f、e分別表示揮舞方向和擺動方向。

通過對式（7）的迭代求解得到式（8）中葉片在各自由度方向的位移。對于葉片在外載荷作用下產生的位移、速度和加速度，采用葉片前三階振動模態的線性關系表達為［7］

式中，d為葉片產生的位移；u1f為一階揮舞方向模態單位特征值；u1e為一階擺振方向模態單位特征值；u2f為二階揮舞方向模態單位特征值。

3 旋轉風輪葉片載荷的數值模擬

根據上述建立的風輪載荷計算模型，在MATLAB軟件中采用數值迭代法對某2MW風力機風輪的載荷進行了計算分析。該風輪由3個葉片組成，每個葉片長30．56m，設計額定風速為14m／s，葉尖速比λ＝5，風輪的設計旋轉速度為2．3rad／s。以葉素理論為基礎，將葉片沿展向劃分成20段，通過對風輪誘導因子的迭代求解，確定風輪各個葉片的空氣動力學載荷［8］。圖2所示為額定風速14m／s，風輪旋轉12s時，即旋轉了4．43r后，所求得的葉片軸向和周向誘導因子在葉片展向的分布。圖3所示為該時刻葉片的空氣動力學載荷法向力和切向力分布，其中切向力主要表現為風輪的輸出功率，法向力主要表現為風輪所承受的推力。風輪在該推力的作用下，在垂直旋轉平面的方向會發生拍打位移，該位移的產生使得葉片的離心力方向發生變化，進而使得葉片的離心力在拍打和擺振方向產生了分量，由于拍打方向的離心力和推力方向相反，因此使得葉片的推力減小，圖4所示為相同時刻條件下計算得到的葉片離心力在拍打和擺振方向的分布。

圖2 葉片誘導因子沿展向的分布

圖3 葉片空氣動力學載荷的分布

圖4 葉片離心力的分布

4 風輪偏移變形特性的算例分析

為了很好地研究風力機葉片的偏移變形特性和驗證該模型的可靠性，對風輪在不同風速和不同位置時刻的葉片的偏移變形特性進行了計算分析。同樣針對上述2MW風力機，在已得的載荷基礎上，在MATLAB軟件中編寫程序對葉片在該載荷下的偏移變形位移進行了數值模擬。葉片在載荷作用下的偏移變形與其所處的運轉位置有關，本文以圖5所示各葉片所處的實線位置為風輪的起始旋轉位置，風輪以圖中箭頭方向即沿逆時針方向旋轉。圖6所示為風輪在額定風速14m／s下，運轉8s時，風輪三個葉片沿展向的不同位移偏移變形。因運轉時間為8s，即剛好旋轉了3r，那么該時刻各葉片所處的位置和起始位置相同。圖7所示為相同風速下，風輪在旋轉12s以后，即旋轉了4．43r，各個葉片處于圖5中虛線位置時，葉片的展向偏移變形。和圖6對比，盡管來流風速和風輪旋轉速度相同，但是由于風輪旋轉時間不同，各葉片計算時刻所處的位置不一樣，因此，三個葉片的偏移變形和圖6所示的偏移變形有著不同。圖8所示為該風輪運轉12s時，風速為10m／s下，風輪各個葉片的偏移變形分布。可以看出，圖8中每個葉片的偏移變形規律和圖7一致，這是由于兩種工況的運轉時間相同，即計算時刻每個葉片所處的位置相同，但是由于圖8風速為10m／s，比圖7中風速14m／s要小，即圖8葉片所承受的載荷要小些，因此葉片展向所產生的偏移變形也相對較小。

圖5 風輪運轉時葉片所處的位置

圖6 額定風速14m／s下風輪旋轉3轉時各個葉片的偏移變形分布

圖7 額定風速14m／s下風輪旋轉4．43轉時各個葉片的偏移變形分布

圖8 額定風速10m／s下風輪旋轉4．43轉時各個葉片的偏移變形分布

5 結論

（1）應用動量葉素理論對風力機風輪旋轉條件下的空氣動力學載荷和離心力進行了理論模型的建立和數值模擬。在充分考慮該載荷特性條件下，基于風輪的結構動力學方程，提出了葉片的偏移變形特性計算數學模型。

（2）針對建立的葉片偏移變形計算數學模型，采用迭代法對某2MW風力機在不同工況下的葉片展向不同位置的偏移變形進行了數值模擬，得到了該風輪在不同風速和不同運轉時刻的各個葉片的偏移變形分布。研究結果對風力機葉片的振動特性和疲勞壽命預測具有理論指導意義。

［1］ 賀德馨．風工程與工業空氣動力學［M］．北京：國防工業出版社，2006．

［2］ Maheri A，Noroozi S，Vinney J．Combined Analytical／FEA－based Coupled Aero Structure Simulation of a Wind Turbine with Bend-twist Adaptive Blades［J］．Renewable Energy，2007，32：916－930．

［3］ Larsen J W，Nielsen S R K．Nonlinear Parametric Instability of Wind Turbine Wings［J］．Journal of Sound and Vibration，2007，299：64－82．

［4］ Sarkar S，Bijl H．Nonlinear Aeroelastic Behavior of an Oscillating Airfoil During Stall－induced Vibration［J］．Journal of Fluids and Structures，2008，24：757－777．

［5］ Turhan O，Bulut G．On Nonlinear Vibrations of a Rotating Beam［J］．Journal of Sound and Vibration，2009，322：314－335．

［6］ Shen W Z，Mikkelsen R，Sorensen J N．Tip Loss Correction for Wind Turbines Computations［J］．Wind Energy，2005，8：457－475．

［7］ Hansen M O L．Aerodynamics of Wind Turbines［M］．London：James ＆ James （Science Publishers）Ltd，2000．

［8］ 陳進，王旭東，沈文忠，等．風力機葉片的形狀優化設計［J］．機械工程學報，2010，46（3）：131－134．