劉 雄，馬新穩(wěn)，沈 世，陳 嚴(yán)

（汕頭大學(xué)能源研究所，廣東 汕頭 515063）

0 引言

風(fēng)力機(jī)通常安裝在條件惡劣的野外環(huán)境中，其動(dòng)力源為隨機(jī)性很強(qiáng)的自然風(fēng)，所以運(yùn)行時(shí)經(jīng)常承受復(fù)雜的隨機(jī)載荷，在復(fù)雜的載荷作用下，導(dǎo)致風(fēng)力機(jī)破壞的主要是失速顫振、經(jīng)典顫振。失速顫振是由葉片失速誘發(fā)的振動(dòng)，由于葉片失速導(dǎo)致的負(fù)氣動(dòng)阻尼誘發(fā)葉片在擺振方向發(fā)生振動(dòng)失穩(wěn)，使風(fēng)力機(jī)非正常停機(jī)或?qū)е氯~片損壞。研究風(fēng)力機(jī)的失速特性，尤其是失速工況下的氣動(dòng)阻尼特性是分析和處理系統(tǒng)穩(wěn)定性和可靠性的關(guān)鍵途徑。

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

隨著新材料和新工藝的發(fā)展，風(fēng)電機(jī)組的大型化、葉片柔性化已成為風(fēng)力機(jī)技術(shù)發(fā)展的方向，把葉片及塔架等理想化為小變形剛體的模型已不再適用。考慮機(jī)組柔性對(duì)于動(dòng)態(tài)氣動(dòng)模型的影響及與之相對(duì)應(yīng)的氣動(dòng)載荷變化對(duì)氣動(dòng)彈性分析的反饋是建立葉片真實(shí)氣動(dòng)阻尼模型的關(guān)鍵因素。

1 葉片的氣動(dòng)阻尼模型

1.1 葉片的有限單元模型

由于葉片的展弦比比較大，可以將葉片簡(jiǎn)化成非均勻懸臂梁，采用二節(jié)點(diǎn)梁?jiǎn)卧獙?duì)葉片進(jìn)行有限元離散，建立葉片的有限單元模型。對(duì)于揮舞和擺振方向，用三次Hermite多項(xiàng)式建立二節(jié)點(diǎn)梁?jiǎn)卧囊恢沦|(zhì)量矩陣（式（1））和一致剛度矩陣（式（2））［6］，其中剛度矩陣由彈性剛度和幾何剛度兩部分組成，考慮了離心剛化作用。

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

由已知的葉片的結(jié)構(gòu)幾何參數(shù)，組集單元質(zhì)量陣和剛度矩陣，根據(jù)瞬時(shí)最小勢(shì)能原理建立葉片在時(shí)變載荷作用下的多自由度運(yùn)動(dòng)方程：

式中，M為系統(tǒng)質(zhì)量矩陣，C為系統(tǒng)阻尼矩陣，K為系統(tǒng)剛度矩陣為扭轉(zhuǎn)振動(dòng)有限單元節(jié)點(diǎn)位移、速度、加速度陣列，為揮舞（擺振）有限單元節(jié)點(diǎn)位移、速度、加速度陣列。

1.2 葉片上的氣動(dòng)載荷

本文主要考慮3個(gè)振動(dòng)角對(duì)葉片氣動(dòng)載荷的影響，即葉片截面扭角變化量Δθ，葉片的揮舞傾角Δβ，葉片擺振傾角σ（擺振傾角初值為0）?？紤]葉片的方位角ψ，不考慮偏航角，機(jī)艙坐標(biāo)系和塔架坐標(biāo)系重合。通過坐標(biāo)變換可以得到槳葉坐標(biāo)系中風(fēng)速分量是：

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

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

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

軸向誘導(dǎo)因子：

切向誘導(dǎo)因子：

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

1.3 葉片上的非氣動(dòng)載荷

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

式中，揮舞方向的Fcflap為離心力載荷，F(xiàn)gflap為重力載荷，對(duì)應(yīng)的下標(biāo)為edge的為對(duì)應(yīng)的擺振方向的載荷；這些載荷的計(jì)算在相關(guān)文獻(xiàn)中可以找到［9］。

1.4 單翼型的氣動(dòng)阻尼模型

當(dāng)阻尼是線性粘滯型阻尼時(shí)，阻尼力對(duì)振動(dòng)系統(tǒng)做功的功率（也可理解為由于阻尼引起的能量耗散的速率）為：

式中，m為振動(dòng)系統(tǒng)的質(zhì)量，為振動(dòng)頻率，ω為固有頻率，q為振幅，c為阻尼系數(shù)，ξ為阻尼比。

當(dāng)阻尼是非線性粘滯形式時(shí)上式不再成立，但每個(gè)周期內(nèi)阻尼力對(duì)振動(dòng)系統(tǒng)做的功是相等的［7］，可用一個(gè)等效功表示：

利用這一點(diǎn)可以計(jì)算等效的阻尼比：

等效阻尼系數(shù)為：

應(yīng)用在求氣動(dòng)阻尼時(shí)，氣動(dòng)阻尼力對(duì)單位長(zhǎng)度翼型在揮舞方向和擺振方向周期內(nèi)做的功分別為［3］（取n個(gè)周期的平均值）：

由此可得單翼型截面的揮舞、擺振氣動(dòng)阻尼系數(shù)：

對(duì)應(yīng)的第n階模態(tài)的模態(tài)氣動(dòng)阻尼比為：

式中，qflap為揮舞方向的振幅，flap為揮舞方向的振動(dòng)頻率，ωnflap為揮舞方向的n階固有頻率，mnflap為揮舞方向的n階模態(tài)質(zhì)量，相應(yīng)的下標(biāo)為edge為擺振方向的相應(yīng)的物理量。

1.5 葉片氣動(dòng)阻尼矩陣的建立

為了便于進(jìn)行結(jié)構(gòu)動(dòng)力學(xué)分析，還需要建立全葉片的氣動(dòng)阻尼矩陣。

一般情況下對(duì)振型進(jìn)行絕對(duì)值最大項(xiàng)歸一化得到葉片的一種規(guī)范化振型。第n階模態(tài)振型中各個(gè)節(jié)點(diǎn)的振幅可以表示為：

式中，qn0為為第n階模態(tài)的最大振幅，qni（i＝1，2，…，N）為第i節(jié)點(diǎn)的振幅。

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

根據(jù)模態(tài)氣動(dòng)阻尼的定義，第n階模態(tài)氣動(dòng)阻尼系數(shù)為：

模態(tài)氣動(dòng)阻尼比為：

Wi為當(dāng)風(fēng)速為v時(shí)一個(gè)周期內(nèi)氣動(dòng)力對(duì)第i單元所做的功：

式中，γi為第i個(gè)單元首距輪轂的距離，γi＋1為第i個(gè)單元尾距輪轂的距離。通常將模態(tài)阻尼表示為對(duì)數(shù)衰減的形式，即：

式中fn為第n階振動(dòng)頻率；同樣第i個(gè)單元的n階模態(tài)氣動(dòng)阻尼比為：

對(duì)數(shù)衰減的形式為：

2 算 例

2.1 程序的整個(gè)迭代過程如下

迭代的終止條件由設(shè)置的收斂條件決定。并且在每個(gè)循環(huán)結(jié)束時(shí)保存相應(yīng)變量。

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

本文研究大型水平軸風(fēng)力機(jī)柔性葉片，所以選用NREL 5MW［10］風(fēng)力機(jī)為計(jì)算實(shí)例，這臺(tái)風(fēng)力機(jī)的葉片參數(shù)如表1所示。

表1 葉片的主要數(shù)據(jù)Table 1 Main datas of the blade

2.2 計(jì)算結(jié)果及結(jié)論

2.2.1 不考慮反饋的氣動(dòng)阻尼的計(jì)算

在穩(wěn)態(tài)失速下氣動(dòng)阻尼的大小在很大程度上依賴于翼型升力系數(shù)隨攻角變化曲線的斜率［2］，而且氣動(dòng)負(fù)阻尼一般出現(xiàn)在低階模態(tài)，氣動(dòng)阻尼系數(shù)用翼型函數(shù)表達(dá)式可以表示為：

本文采用能量損失法通過式（34）計(jì)算了關(guān)鍵節(jié)點(diǎn)的氣動(dòng)阻尼沿不同風(fēng)速的分布如圖3（a）所示，與翼型函數(shù)表達(dá)式預(yù)測(cè)的結(jié)果一致。

圖3 關(guān)鍵節(jié)點(diǎn)只考慮一種變形影響的一階揮舞和擺振阻尼對(duì)數(shù)衰減率Fig.3 Logarithmic decrement of 1st flap－wise and edge－wise mode damping for the important section considering one deformation

2.2.2 各種振動(dòng)變形對(duì)氣動(dòng)阻尼的影響分析

本文的計(jì)算結(jié)果包括以下7種情況，見表2。

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

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

從圖4中可以看出：揮舞與扭轉(zhuǎn)方向振動(dòng)變形的耦合作用在低風(fēng)速時(shí)對(duì)氣動(dòng)阻尼的影響不大，當(dāng)風(fēng)速增大時(shí)會(huì)使揮舞方向的氣動(dòng)阻尼減小。擺振與扭轉(zhuǎn)變形的耦合作用同樣在低風(fēng)速時(shí)對(duì)擺振方向的氣動(dòng)阻尼影響不大，當(dāng)風(fēng)速增大時(shí)會(huì)使擺振方向的氣動(dòng)阻尼增大。擺振—扭轉(zhuǎn)耦合對(duì)擺振方向的氣動(dòng)阻尼的影響比揮舞—扭轉(zhuǎn)耦合對(duì)揮舞方向阻尼的影響明顯。

圖4 關(guān)鍵節(jié)點(diǎn)考慮彎曲扭轉(zhuǎn)耦合變形的一階揮舞和擺振阻尼對(duì)數(shù)衰減率Fig.4 Logarithmic decrement of 1st flap－wise and edge－wise mode damping for the important section considering deformation coupling

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

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

3 結(jié)論

本文基于能量損失法建立了大型風(fēng)力機(jī)柔性葉片氣動(dòng)阻尼的分析模型。通過把結(jié)構(gòu)的振動(dòng)變形反饋到氣動(dòng)模型中，分析了其在揮舞、擺振、扭轉(zhuǎn)振動(dòng)變形影響下的氣動(dòng)阻尼特性。從得出的結(jié)論可以知道柔性葉片的氣動(dòng)阻尼相對(duì)于剛性葉片有較大變化，各種振動(dòng)變形對(duì)各個(gè)截面氣動(dòng)阻尼大小產(chǎn)生了很大的影響。

［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］王小虎.水平軸風(fēng)力機(jī)氣動(dòng)阻尼特性分析與優(yōu)化設(shè)計(jì)［D］.汕頭大學(xué)，2010.

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

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

［8］劉雄，陳嚴(yán)，葉枝全.水平軸風(fēng)力機(jī)氣動(dòng)性能計(jì)算模型［J］.太陽能學(xué)報(bào)，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.