大型風(fēng)力機氣動彈性響應(yīng)計算研究

呂計男，劉子強，趙 玲，冉景洪

（中國航天空氣動力技術(shù)研究院，北京 100074）

0 引 言

隨著風(fēng)力發(fā)電事業(yè)的蓬勃發(fā)展，風(fēng)力機的功率已經(jīng)由初期的千瓦級小型風(fēng)力機上升為現(xiàn)在的兆瓦級大型風(fēng)力機，風(fēng)力機葉片長度和重量都呈現(xiàn)快速增加的趨勢。風(fēng)力機葉片的氣動彈性分析需要氣動力載荷及結(jié)構(gòu)數(shù)據(jù)。目前，氣動力載荷計算主要以動量葉素理論（BEM）及計算空氣動力學(xué)（CFD）方法為主。BEM方法簡單快捷，應(yīng)用方便，理論成熟度高，一直是國內(nèi)外風(fēng)力機氣動載荷及氣動彈性問題研究的重要手段。CFD方法可以提供精確的流場描述，但對計算機硬件的要求很高，計算效率較低，特別是對氣動彈性分析所需的非定常氣動力計算尚存在很大難度。

本文選用南京航空航天大學(xué)設(shè)計的NH1500風(fēng)輪葉片作為研究對象。該葉片長度40.5m，三維構(gòu)型由南航優(yōu)化后NH02系列翼型疊加組成。針對大型風(fēng)力機葉片氣動彈性響應(yīng)計算，本課題組開發(fā)了一套基于BEM方法的氣動彈性分析程序，該程序可以考慮偏航工況、計及葉尖及輪轂損失、考慮翼型厚度和寬度、風(fēng)剪切的影響，并可以進(jìn)行失速修正。本文針對NH1500葉片，比較了本程序計算的氣動載荷結(jié)果與商用軟件GH Bladed的結(jié)果，驗證了程序可靠性，并給出風(fēng)力機整機結(jié)構(gòu)在氣動載荷作用下的響應(yīng)結(jié)果。

1 葉片動量葉素理論氣動力建模

1.1 理論基礎(chǔ)

動量理論［1］主要用來估算風(fēng)力機的理想輸出功率，該理論描述的是風(fēng)力機葉輪整體與空氣氣流的能量交換情況。葉素理論［1］將風(fēng)力機葉片簡化為有限個葉素疊加的形式，葉片三維氣動特性由葉素的氣動性能沿徑向積分得到，該理論以葉素表面氣體流動的狀態(tài)分析葉片的受力和能量交換?？紤]到采用動量理論與葉素理論描述的葉片氣動力與力矩在物理本質(zhì)上的同一性，可以建立平衡方程得出動量葉素理論中軸向及切向速度誘導(dǎo)因子的顯式表達(dá)式［2］。

在理論的基礎(chǔ)上，可以通過葉尖損失系數(shù)、輪轂損失系數(shù)考慮由于葉尖及根部氣流沿葉片產(chǎn)生的二次流動而引起的力矩的減?。?，3］。

葉柵理論（Cascade Theory）研究葉片厚度和寬度的影響導(dǎo)致的葉素攻角的改變［1］。當(dāng)周向速度誘導(dǎo)因子大于0.5時，風(fēng)輪工作在湍流尾流狀態(tài)，這時需要對動量理論進(jìn)行修正。本文的計算模型采用的是 Glauert提出的修正表達(dá)式［4－5］：

1.2 風(fēng)力機坐標(biāo)系統(tǒng)建模

根據(jù)動量葉素理論得到的風(fēng)力機葉片所受的氣動力與力矩是相當(dāng)于葉素局部坐標(biāo)系統(tǒng)的。由于風(fēng)力機氣動性能最優(yōu)設(shè)計的要求，一般情況下，風(fēng)輪轉(zhuǎn)軸傾角，槳葉錐角不為零，這樣使得氣動力／力矩需要完成從葉素局部坐標(biāo)到全局坐標(biāo)的坐標(biāo)變換。此外，葉片旋轉(zhuǎn)過程中，工作方位角對坐標(biāo)系統(tǒng)轉(zhuǎn)換關(guān)系的影響也需要考慮。文獻(xiàn)［2］中對該部分內(nèi)容進(jìn)行了詳盡的討論。

大型風(fēng)力機的葉片很長，槳葉不同位置的風(fēng)速相差很大，一般情況下需要考慮風(fēng)剪切的影響。本文采用的是指數(shù)模型［4，6］：

其中，H代表葉素高度；HR為參考高度；UR為參考高度處的風(fēng)速；ε為經(jīng)驗風(fēng)剪切指數(shù)，本文計算模型取值為0.1667。

1.3 計算流程

利用動量葉素理論獲得葉片氣動力主要分為速度誘導(dǎo)因子迭代求解及葉片氣動性能計算兩部分內(nèi)容，文獻(xiàn)［7］中詳盡了敘述了速度誘導(dǎo)因子迭代的求解過程，這里由于篇幅限制，不再累述。程序的總流程如圖1所示。

圖1 風(fēng)力機結(jié)構(gòu)響應(yīng)求解流程圖Fig.1 Structure response solution scheme

1.4 風(fēng)力機葉片氣動力模型結(jié)果分析

南京航空航天大學(xué)針對NH1500葉片采用GH Bladed軟件進(jìn)行了葉片氣動性能的校核。NH1500葉片為南航自主設(shè)計，長度40.5m，截面翼型沿葉片展向分別為 NH02＿40、NH02＿35、NH02＿30、NH02＿25、NH02＿21、NH02＿18、NH02＿15，翼型形狀及氣動性能數(shù)據(jù)為南航提供。針對該葉片，通過編寫程序，完成了葉片氣動力的建模工作，對比結(jié)果顯示自編程序的計算結(jié)果與南航提供的采用GH Bladed軟件計算的結(jié)果吻合的很好，驗證了自編程序的正確性，計算結(jié)果如圖2所示。

圖2 功率系數(shù)計算結(jié)果對比圖Fig.2 Power coefficient contrast

圖2中縱坐標(biāo)表示功率系數(shù)，橫坐標(biāo)表示來流風(fēng)速。CpCAAA為本課題組自編程序得到的功率系數(shù)數(shù)據(jù)，CpNH為南航采用GH Bladed軟件計算得到的功率系數(shù)數(shù)據(jù)。

2 風(fēng)力機全機有限元建模

2.1 葉片結(jié)構(gòu)有限元建模

風(fēng)力機整體結(jié)構(gòu)在氣動力作用下的響應(yīng)問題是風(fēng)力機氣動彈性研究的一個重要方向。結(jié)構(gòu)設(shè)計需要根據(jù)響應(yīng)結(jié)果獲得彎矩載荷等重要數(shù)據(jù)。

風(fēng)力機葉片外形沿葉片展向變化且內(nèi)部結(jié)構(gòu)復(fù)雜，進(jìn)行結(jié)構(gòu)三維有限元建模難度很大，使得分布?xì)鈩恿ψ饔孟碌慕Y(jié)構(gòu)響應(yīng)計算的代價極大。對響應(yīng)問題，局部變形對整體響應(yīng)結(jié)果的影響通常不予考慮。借鑒航空氣動彈性問題的做法，將大展弦比機翼結(jié)構(gòu)簡化為梁模型處理，這里將大型風(fēng)力機葉片也做簡化處理。利用南航提供的結(jié)構(gòu)剛度及質(zhì)量分布數(shù)據(jù)通過模態(tài)分析可以獲得葉片的結(jié)構(gòu)動力學(xué)特性。葉片簡單梁模型與三維有限元模型的結(jié)構(gòu)動力學(xué)特性對比如表1所示。

對比結(jié)果顯示，除了第四階固有振動的頻率誤差略大之外，其它各階頻率誤差都在有限元模型建模允許誤差范圍之內(nèi)。

2.2 塔架結(jié)構(gòu)有限元建模

進(jìn)行風(fēng)力機全機的響應(yīng)計算，柔性塔架的影響需要考慮。目前，南航設(shè)計的大型風(fēng)力機只提供了葉片的結(jié)構(gòu)屬性，無塔架屬性可供參考。本部分內(nèi)容主要根據(jù)塔架剛度設(shè)計準(zhǔn)則［8］，對塔架構(gòu)成部件及邊界條件進(jìn)行簡化，建立了塔架的有限元模型。

表1 葉片簡單梁模型與三維有限元模型的結(jié)構(gòu)動力學(xué)特性對比圖Table1 Structural dynamic characteristic of beam and 3DFE model

塔架結(jié)構(gòu)主要包括塔體及塔頂集中質(zhì)量。根據(jù)剛度設(shè)計準(zhǔn)則，塔架結(jié)構(gòu)的一階固有頻率以及其上10%，和三階固有頻率以及其下10%，需要避開風(fēng)力機工作頻率的1倍頻和3倍頻。南航設(shè)計的風(fēng)力機額定轉(zhuǎn)速17.2RPM／min，即工作頻率為0.287Hz，三倍頻為0.861Hz。根據(jù)剛度設(shè)計準(zhǔn)則，塔架固有振動頻率需要避開的頻率范圍為：0.287Hz～0.316Hz，0.775Hz～0.861Hz。

本文關(guān)于塔架的有限元建模沒有考慮法蘭上的螺栓、塔架內(nèi)部的附屬結(jié)構(gòu)等構(gòu)件；在塔架頂端創(chuàng)建一個質(zhì)量單元節(jié)點，模擬塔頂質(zhì)量，質(zhì)量點的位置設(shè)置在機艙、輪轂和葉片的合重心位置處；此外，塔架基礎(chǔ)按固支邊界條件處理。

按照以上條件，塔架的結(jié)構(gòu)動力學(xué)特性如表2所示。

表2 塔架的結(jié)構(gòu)動力學(xué)特性Table2 Structural dynamic characteristic of tower

2.3 風(fēng)力機全機建模

通過商用有限元軟件建立風(fēng)力機全機模型，葉片結(jié)構(gòu)采用梁單元模擬，塔頂質(zhì)量采用點質(zhì)量單元模擬。葉片工作方位角每轉(zhuǎn)過10°動量葉素氣動力求解模塊求解一次氣動力并保存結(jié)果到數(shù)據(jù)文件。結(jié)構(gòu)分析軟件讀入氣動力數(shù)據(jù)，氣動載荷以數(shù)據(jù)表形式施加到結(jié)構(gòu)單元上，含邊界條件的有限元模型如圖3所示。

圖3 風(fēng)力機全機有限元模型示意圖Fig.3 FE model of total wind turbine

3 計算結(jié)果與討論

采用軟件的瞬態(tài)動力學(xué)分析模塊，結(jié)構(gòu)動力學(xué)方程求解采用直接積分法中的Newmark方法，結(jié)構(gòu)阻尼采用Rayleigh阻尼形式。Rayleigh阻尼中的α，β為不依賴于頻率的常數(shù)，且阻尼對穩(wěn)態(tài)振動的振幅沒有影響。

考慮柔性塔架作用時，參考高度處來流風(fēng)速為12m／s時，葉片葉尖節(jié)點振動位移隨時間的變化如圖4－圖5所示。由于初始條件設(shè)置為零位移、零速度，振動的初始階段會產(chǎn)生跳躍，是強迫振動的過渡階段，隨著時間的進(jìn)行，振動逐漸進(jìn)入簡諧激勵下強迫振動的穩(wěn)態(tài)階段。

3.1 揮舞方向振動

葉片X方向及Z（揮舞）方向振動的時程如圖4所示。X方向位移最大值27.27cm，最小值23.91cm，平衡位置為25.59cm，振幅1.68cm。揮舞方向位移最大值20.07cm，最小值17.77cm，平衡位置為18.92cm，振幅1.15cm。

圖4 葉片葉尖X及Z方向振動位移時程圖Fig.4 Response of blade in Xand Zdirection

3.2 擺振方向振動

Y（擺振）方向位移最大值4.51m，最小值4.14m，平衡位置為4.33m，振幅18.85cm，時程如圖5所示。

圖5 葉片葉尖Y方向振動位移時程圖Fig.5 Response of blade in Y－direction

3.3 塔架振動

塔架的振動產(chǎn)生結(jié)構(gòu)內(nèi)部載荷，塔架的振動情況將影響結(jié)構(gòu)的彎矩載荷及疲勞設(shè)計。圖6給出塔架頂端節(jié)點的時域振動情況。不考慮塔架對葉片的影響，如葉片經(jīng)過塔架時加速逆風(fēng)效應(yīng)時，塔架頂端迎風(fēng)方向位移最大值7.37cm，最小值6.97cm，平衡位置為7.17cm，振幅0.2cm，時程如圖6所示。

圖6 塔架頂端振動位移時程圖Fig.6 Response of tower in time domain

3.4 幾何非線性

大型風(fēng)力機葉片長度非常大，受力時會發(fā)生大位移小應(yīng)變的幾何非線性效應(yīng)，這時平衡方程和幾何關(guān)系都是非線性的。早期的非線性有限元分析基本上是線性分析的擴(kuò)展，只能針對個別具體問題進(jìn)行分析。近年來基于非線性連續(xù)介質(zhì)力學(xué)原理的有限元分析有了很大的發(fā)展，商用軟件也有相應(yīng)成熟的模塊求解該類問題。本部分利用商用有限元軟件分別進(jìn)行了葉片線性響應(yīng)分析與幾何非線性響應(yīng)分析，對比結(jié)果如圖7和表3所示。

結(jié)果顯示，對葉片葉尖節(jié)點，考慮幾何非線性時葉片的變形較小，呈現(xiàn)出比線性分析結(jié)果剛硬的性質(zhì)，這是葉片初位移效應(yīng)及初應(yīng)力效應(yīng)共同作用的結(jié)果。

圖7 幾何非線性對葉片振動的影響Fig.7 Influence of geometry nonlinearity

表3 幾何非線性對葉片振動的影響Table3 Influence of geometry nonlinearity

4 結(jié) 論

本文建立了一套適用于大型風(fēng)力機氣動彈性響應(yīng)的快速計算方法。該方法可以對葉片所受氣動力進(jìn)行快速預(yù)測，并快速的給出結(jié)構(gòu)在氣動力作用下的響應(yīng)情況，提供了一種可供風(fēng)力機設(shè)計人員借鑒的技術(shù)手段。

［1］ROBERT E WILSON，PETER B SLISSAMAN，STEL N WALKER.Aerodynamic performance of wind turbines［M］.Department of Mechanical Engineering，Oregon State University，1976.

［2］劉雄，陳嚴(yán)，葉枝全.水平軸風(fēng)力機氣動性能計算模型［J］.太陽能學(xué)報，2005，26（6）：792－800.

［3］PRANDTL L，TIETJENS O G.Applied hydro and aeromechanics［M］.Dover Publications，1957.

［4］DET NORSKE VERITAS and RIS NATIONAL LABORATORY.Guidelines for design of wind turbines［M］.Det Norske Veritas and Ris National Laboratory，2002.

［5］GLAUERT H.The analysis of experimental results in the windmill brake and vortex ring states of an airscrew［R］.Reports and Memoranda，1926.

［6］SPERA D A.Wind turbine technology［M］.New York ASME Press，1994.

［7］ROBERTS E.WILSON，PETER B.S.LISSAMAN.Applied aerodynamics of wind power machines［M］.Or－egon State University Corvallis，1974.

［8］趙立新.風(fēng)力發(fā)電機塔架的有限元分析與優(yōu)化設(shè)計［D］.吉林：吉林大學(xué)，2008.