風力機柔性葉片振動變形對其氣動阻尼的影響分析

劉 雄，馬新穩，沈 世，陳 嚴

（汕頭大學能源研究所，廣東 汕頭 515063）

0 引言

風力機通常安裝在條件惡劣的野外環境中，其動力源為隨機性很強的自然風，所以運行時經常承受復雜的隨機載荷，在復雜的載荷作用下，導致風力機破壞的主要是失速顫振、經典顫振。失速顫振是由葉片失速誘發的振動，由于葉片失速導致的負氣動阻尼誘發葉片在擺振方向發生振動失穩，使風力機非正常停機或導致葉片損壞。研究風力機的失速特性，尤其是失速工況下的氣動阻尼特性是分析和處理系統穩定性和可靠性的關鍵途徑。

國外對風力機的氣動阻尼的研究開始于20世紀90年代，1994年Hasen A C率先通過翼型的法向力函數微分表達式分析了蹺蹺板約束的風力機葉片的穩態失速下的氣動阻尼［1］，同時基于周期內取平均值的方法分析了受動態入流、非線性靜態失速、旋轉效應影響的動態失速氣動阻尼特性。1997 年Rasumussen F，Petersen J T 等人重點分析了動態失速下的風力機的氣動阻尼特性［2］，通過翼型升力函數表達式定性分析了氣動阻尼受各種因素的影響機理，并引入了能量損失法來計算在多因素影響下的葉片的氣動阻尼特性，同時考慮了變槳振動對氣動阻尼的影響。1998年Risφ實驗室系統分析了風力機動態失速下的氣動阻尼［3］，并為了減輕氣動載荷，防止負阻尼的出現提出了對翼型的優化方法。1999年Bertagnolio F，Gaunaa M 等人使用CFD 的方法對基于納維——斯托克流體方程與基于半經驗模型的氣動阻尼的計算方法進行比較分析［4］。國內，上海交通大學通過流場測量和CFD 分析，從實驗測量和數值計算方面研究了穩態失速和動態失速問題。汕頭大學在對現代大型風力機進行了動態氣彈穩定性分析的基礎上［5］，針對防止負氣動阻尼產生，對風力機翼型進行了優化。但以上進行的研究都是基于風力機小變形假設的，沒有考慮葉片的振動變形對氣動阻尼的影響。

隨著新材料和新工藝的發展，風電機組的大型化、葉片柔性化已成為風力機技術發展的方向，把葉片及塔架等理想化為小變形剛體的模型已不再適用。考慮機組柔性對于動態氣動模型的影響及與之相對應的氣動載荷變化對氣動彈性分析的反饋是建立葉片真實氣動阻尼模型的關鍵因素。

1 葉片的氣動阻尼模型

1.1 葉片的有限單元模型

由于葉片的展弦比比較大，可以將葉片簡化成非均勻懸臂梁，采用二節點梁單元對葉片進行有限元離散，建立葉片的有限單元模型。對于揮舞和擺振方向，用三次Hermite多項式建立二節點梁單元的一致質量矩陣（式（1））和一致剛度矩陣（式（2））［6］，其中剛度矩陣由彈性剛度和幾何剛度兩部分組成，考慮了離心剛化作用。

式中，ρe為單元密度，l為單元長度，Ae為截面面積，Ee為彈性模量，Ie為截面慣性矩，Ne為單元所受軸向力。結構阻尼采用Rayleigh 阻尼模型［7］。

由已知的葉片的結構幾何參數，組集單元質量陣和剛度矩陣，根據瞬時最小勢能原理建立葉片在時變載荷作用下的多自由度運動方程：

式中，M為系統質量矩陣，C為系統阻尼矩陣，K為系統剛度矩陣為扭轉振動有限單元節點位移、速度、加速度陣列，為揮舞（擺振）有限單元節點位移、速度、加速度陣列。

1.2 葉片上的氣動載荷

本文主要考慮3個振動角對葉片氣動載荷的影響，即葉片截面扭角變化量Δθ，葉片的揮舞傾角Δβ，葉片擺振傾角σ（擺振傾角初值為0）。考慮葉片的方位角ψ，不考慮偏航角，機艙坐標系和塔架坐標系重合。通過坐標變換可以得到槳葉坐標系中風速分量是：

式中，β＝β0＋Δβ為葉片傾角，β0 為葉片傾角初值，u為來流風速；為了精確的得到葉片的氣動載荷，本文采用葉素動量理論計算葉片的氣動載荷并考慮葉尖損失和輪轂損失修正、以及葉柵理論和軸向誘導因子修正［8］。

圖1 葉片截面圖Fig.1 Blade cross－section schematic

單位長度上的氣動力如圖1所示：E為截面扭轉的剛心，本文統一取距前沿1／4處為各截面剛心，o為氣動中心，α為攻角，θ為扭角，L為升力，D為阻力，M為繞剛心扭矩，Δθ為葉片截面扭角變化量，葉片振動在揮舞和擺振方向的速度分量分別和單位長度翼型截面揮舞方向、擺振方向的氣動載荷分別為：

軸向誘導因子：

切向誘導因子：

式中，F為葉尖損失系數和輪轂損失系數的乘積，ρ為空氣密度，c為弦長，r為截面所在半徑，e為氣動中心到剛心的有效距離，θ0為截面扭角初值，Ω為風輪轉速，CL，CD，CM為升力系數、阻力系數、力矩系數。

1.3 葉片上的非氣動載荷

加入葉片的非氣動載荷可得揮舞和擺振作用力及扭矩分別為：

式中，揮舞方向的Fcflap為離心力載荷，Fgflap為重力載荷，對應的下標為edge的為對應的擺振方向的載荷；這些載荷的計算在相關文獻中可以找到［9］。

1.4 單翼型的氣動阻尼模型

當阻尼是線性粘滯型阻尼時，阻尼力對振動系統做功的功率（也可理解為由于阻尼引起的能量耗散的速率）為：

式中，m為振動系統的質量，為振動頻率，ω為固有頻率，q為振幅，c為阻尼系數，ξ為阻尼比。

當阻尼是非線性粘滯形式時上式不再成立，但每個周期內阻尼力對振動系統做的功是相等的［7］，可用一個等效功表示：

利用這一點可以計算等效的阻尼比：

等效阻尼系數為：

應用在求氣動阻尼時，氣動阻尼力對單位長度翼型在揮舞方向和擺振方向周期內做的功分別為［3］（取n個周期的平均值）：

由此可得單翼型截面的揮舞、擺振氣動阻尼系數：

對應的第n階模態的模態氣動阻尼比為：

式中，qflap為揮舞方向的振幅，flap為揮舞方向的振動頻率，ωnflap為揮舞方向的n階固有頻率，mnflap為揮舞方向的n階模態質量，相應的下標為edge為擺振方向的相應的物理量。

1.5 葉片氣動阻尼矩陣的建立

為了便于進行結構動力學分析，還需要建立全葉片的氣動阻尼矩陣。

一般情況下對振型進行絕對值最大項歸一化得到葉片的一種規范化振型。第n階模態振型中各個節點的振幅可以表示為：

式中，qn0為為第n階模態的最大振幅，qni（i＝1，2，…，N）為第i節點的振幅。

將上式中的振幅代入（25）式得：

根據模態氣動阻尼的定義，第n階模態氣動阻尼系數為：

模態氣動阻尼比為：

Wi為當風速為v時一個周期內氣動力對第i單元所做的功：

式中，γi為第i個單元首距輪轂的距離，γi＋1為第i個單元尾距輪轂的距離。通常將模態阻尼表示為對數衰減的形式，即：

式中fn為第n階振動頻率；同樣第i個單元的n階模態氣動阻尼比為：

對數衰減的形式為：

2 算 例

2.1 程序的整個迭代過程如下

迭代的終止條件由設置的收斂條件決定。并且在每個循環結束時保存相應變量。

圖2 程序流程圖Fig.2 Program flow chart

本文研究大型水平軸風力機柔性葉片，所以選用NREL 5MW［10］風力機為計算實例，這臺風力機的葉片參數如表1所示。

表1 葉片的主要數據Table 1 Main datas of the blade

2.2 計算結果及結論

2.2.1 不考慮反饋的氣動阻尼的計算

在穩態失速下氣動阻尼的大小在很大程度上依賴于翼型升力系數隨攻角變化曲線的斜率［2］，而且氣動負阻尼一般出現在低階模態，氣動阻尼系數用翼型函數表達式可以表示為：

本文采用能量損失法通過式（34）計算了關鍵節點的氣動阻尼沿不同風速的分布如圖3（a）所示，與翼型函數表達式預測的結果一致。

圖3 關鍵節點只考慮一種變形影響的一階揮舞和擺振阻尼對數衰減率Fig.3 Logarithmic decrement of 1st flap－wise and edge－wise mode damping for the important section considering one deformation

2.2.2 各種振動變形對氣動阻尼的影響分析

本文的計算結果包括以下7種情況，見表2。

表2 結果包含的七種情況Table 2 Seven conditions of the results

從圖3中可以看出：揮舞傾角對揮舞方向的氣動阻尼影響最大，扭轉傾角對其影響次之，擺振傾角對其影響最小，幾乎沒有什么影響；這三種振動變形都會使揮舞方向的氣動阻尼減小。擺振傾角對擺振方向的氣動阻尼影響最大，揮舞傾角對其影響次之，扭轉傾角對其影響最小；這三種振動變形都會使擺振方向氣動阻尼增大。

從圖4中可以看出：揮舞與扭轉方向振動變形的耦合作用在低風速時對氣動阻尼的影響不大，當風速增大時會使揮舞方向的氣動阻尼減小。擺振與扭轉變形的耦合作用同樣在低風速時對擺振方向的氣動阻尼影響不大，當風速增大時會使擺振方向的氣動阻尼增大。擺振—扭轉耦合對擺振方向的氣動阻尼的影響比揮舞—扭轉耦合對揮舞方向阻尼的影響明顯。

圖4 關鍵節點考慮彎曲扭轉耦合變形的一階揮舞和擺振阻尼對數衰減率Fig.4 Logarithmic decrement of 1st flap－wise and edge－wise mode damping for the important section considering deformation coupling

圖5 單葉片考慮變形影響的一階揮舞和一階擺振阻尼對數衰減率Fig.5 Logarithmic decrement of 1st flap－wise and edge－wise mode damping for the single blade considering deformation

從圖5中可以看出：不論是哪種振動變形只會影響各截面影響氣動阻尼的大小，可以近似認為振動變形不會影響氣動阻尼沿葉片的分布。

3 結論

本文基于能量損失法建立了大型風力機柔性葉片氣動阻尼的分析模型。通過把結構的振動變形反饋到氣動模型中，分析了其在揮舞、擺振、扭轉振動變形影響下的氣動阻尼特性。從得出的結論可以知道柔性葉片的氣動阻尼相對于剛性葉片有較大變化，各種振動變形對各個截面氣動阻尼大小產生了很大的影響。

［1］HASEN A C.Aerodynamic damping of blade flap motions at high angles of attack［J］.JournalofSolarEnergyEngineering，TransactionsoftheASME，1995，117（3）：194－199.

［2］RASMUSSEN F，PETERSEN J T，MADSEN H A.Dynamic stall and aerodynamic damping［R］.Denmark：RisφNational Laboratory，1997.

［3］PETERSEN J T，MADSEN H A，BJORCK A，et al.Prediction of dynamic loads and induced vibrations in stall［R］.Denmark：RisφNational Laboratory，1998.

［4］BERTAGNOLIO F，GAUNAA M，SRENSEN N N.Computation of modal aerodynamic damping using CFD［J］.ASME2003WindEnergySymposium，2003：115－125.

［5］王小虎.水平軸風力機氣動阻尼特性分析與優化設計［D］.汕頭大學，2010.

［6］王勖成.有限單元法［M］.北京：清華大學出版社，2003：56－72.

［7］CLOUGH R W，PENZIEN J.Dynamics of structures［M］.New York：McGraw Hill，1993：53－58.

［8］劉雄，陳嚴，葉枝全.水平軸風力機氣動性能計算模型［J］.太陽能學報，2005，26（6）：792－800.

［9］MARTIN O L HASEN.Aerodynamics of wind turbines［M］.James ＆James Ltd，2000：104－124.

［10］JOKANMA J S.Butterfield definition of a 5MW reference wind turbine for offshore system development［R］.NREL／TP 500－3806.National Renewable Energy Laboratory，Goldern，co，2009.