徐 磊，李德源

(廣東工業(yè)大學機電工程學院，廣州 510006)

引 言

單機大型化是現(xiàn)代風力機發(fā)展的必然方向，但柔性構件(如葉片等)振動對于機組氣動力的影響及與之相對應的氣動載荷變化對機組氣動彈性的反饋成為不能忽略的問題。在葉片設計過程中，需要準確地知道葉片各截面在空氣動力、重力與慣性力作用下的內力，因此，尋求包含非線性變形效應的風力機葉片各截面內力分析的方法成為了大型風力機柔性葉片強度與剛度校核與優(yōu)化設計的重要基礎［1-2］。

對于風力機葉片這種具有空間運動和彈性變形的構件，Molenaar［3］及 Holierhoek［4］等人引入“超級單元”(Superelement)對其進行離散，以反映其彈性變形。Molenaar［3］采用超級單元(無扭轉自由度)對風力機的柔性構件(如風輪葉片和塔架)進行離散，各剛體之間通過力元(彈簧、阻尼器)與鉸連接，這樣風力機系統(tǒng)可以用一離散的剛、柔混合多體系統(tǒng)(HMBS)來表示。Holierhoek［4］使用HMBS方法建立了 NM80風力機的力學模型，利用哈密頓原理導出了動力學方程，分析了超級單元個數(shù)及幾種不同的計算超級單元彈簧系數(shù)的方法對柔性件的固有頻率和阻尼系數(shù)計算精度的影響。

對于葉片空氣動力載荷的計算，目前常用的理論有葉素動量理論(Blade Element-Momentum Theory，BEM)、廣義動態(tài)尾流理論(Generalized Dynamic Wake Theory，GDW)以及計算流體動力學(Computational Fluid Dynamics，CFD)方法等［5］。由于 BEM 理論求解快速，且可結合其他非定常氣動模型，計算結果合理，應用較廣。GDW理論除了計算旋轉平面到尾流的軸向誘導速度外，還考慮了尾跡效應，氣動模型較為復雜［6］。而CFD方法則使用數(shù)值方法求解非線性微分方程組(Navier-Stokes方程)［7］，雖然計算精度較高，但計算量較大。

本文采用HMBS方法建立葉片多體系統(tǒng)的力學模型，基于虛位移原理導出了該系統(tǒng)的動力學方程的一般形式，編制了通用的動力學仿真程序。在此基礎上聯(lián)合BEM氣動模型得到氣彈耦合模型，在對動力學方程進行數(shù)值積分每一時間步中，氣動模型能實時地根據(jù)葉片各剛體的速度、角速度和姿態(tài)坐標(方向余弦陣)計算出作用在葉片各剛體上的氣動載荷，并將其加載到瞬時變形的葉片中，實現(xiàn)了葉片的氣彈耦合時域響應分析。

算例以美國可再生能源實驗室(NREL)發(fā)布的5 MW近海風力機葉片為研究對象［8］，計算了風力機葉片在運行工況下?lián)]舞、擺振與扭轉力矩的時域響應。所作研究對大型風力機葉片設計、優(yōu)化與保障風力機組的安全穩(wěn)定運行，具有重要的應用價值及意義。

1 柔性葉片的多體系統(tǒng)動力學建模

應用計算多體系統(tǒng)動力學中的Roberson-Wittenburg(R-W)方法［9］來建立柔性葉片的多體動力模型。為了準確地描述葉片在氣彈耦合分析中表現(xiàn)出的柔性行為，在多體系統(tǒng)模型中引入所謂的“超級單元”對柔性葉片進行動力學建模［4］，并將變形較小的輪轂和機艙等構件處理為剛體［10］，可在盡量少的自由度下得到與實際相符的結果。每個超級單元包含4個體，單元之間通過相鄰體聯(lián)接，而單元內部體通過帶彈簧和阻尼器的萬向節(jié)或轉動鉸相聯(lián)接，使其能反映葉片的橫向彎曲變形與扭轉變形。這樣葉片被處理成由剛體、鉸、力元和外力組成的多體系統(tǒng)。

1.1 超級單元模型和葉片的離散

將柔性構件離散為若干個超級單元，如圖1所示。其中B1和B2，B3和B4均由萬向節(jié)連接，而B2和B3則由旋轉鉸連接，由此來反映構件的彎曲與扭轉變形，彈簧與阻尼器約束剛體間的相對運動(如圖中的Cy1～Cy3)。因此，每個超級單元有4個旋轉和1個扭轉共5個自由度，圖1中兩個相互垂直平面內的彎曲由2個萬向節(jié)的廣義坐標 θ1，θ2，θ4，θ5來表示，扭轉則由轉動鉸的廣義坐標θ3來表示。

圖1 超級單元模型

以NREL 5 MW水平軸風力機61.5 m長葉片為算例，將其離散為4個超級單元來建立柔性葉片的拓撲模型。由于前后兩個超級單元中相鄰的兩剛體剛性聯(lián)接，可合并為一個剛體(如B4、B7和B10)，這樣葉片可分割成13個剛體，共21個自由度。在地面上建立慣性坐標系XYZ，在葉片的根部建立葉片坐標系X'Y'Z'(葉片的揮舞和擺振狀態(tài)的參考坐標系)，OO'的距離與輪轂質心到葉根的距離相同，以模擬葉片繞風輪主軸的轉動。

1.2 柔性葉片的多體系統(tǒng)動力學方程

根據(jù)R-W方法，按照規(guī)則標號對圖2所示的多體系統(tǒng)中各個體及鉸進行標號，體記作Bi(i=1，…，13)，其內接鉸記作Hi(i=1，…，13)。鉸H1的鉸點建立在慣性坐標系原點O上，其余內接鉸Hi的鉸點均建立在各剛體質心截面的弦長1/4處。以B1代表葉片上與輪轂固結的剛體，其運動被約束，對于其余體Bi，其內接鉸Hi的相對運動可通過固結在體Bi及其內接體上的坐標系的相對關系來描述。

圖2 剛體規(guī)則標號、慣性坐標系XYZ和葉片坐標系X'Y'Z'

取各鉸的轉動為廣義坐標，則從柔性葉片多體系統(tǒng)的拓撲構形可得其廣義坐標陣為:

式中，k為質點系質點數(shù)［rk為質點廣義坐標;δrk為其虛位移;為廣義加速度;Fk為外力主矢與主矩構成的廣義力。由此可導出葉片系統(tǒng)以q(t)為變量的二階常微分方程組:

其中，Z稱為廣義質量陣;z為廣義力陣，其具體表達形式可參閱文獻［9］。該方程組中廣義質量和廣義力陣均包含了以廣義坐標q為變量的方向余弦陣，而且這些變量隨時間變化，因此該方程組呈非線性。結合相應的初始條件，對上式進行數(shù)值積分求解，即可求得葉片中各個體的運動量q(t)、(t)及(t)。

1.3 模型的相關參數(shù)與模型驗證

柔性葉片離散為通過彈簧盒阻尼器連接的多個剛體的過程中，關鍵問題在于彈簧的剛度系數(shù)如何確定才能使模型能盡量符合柔性葉片的力學特性，確定原則如下:(1)在靜載荷下，離散模型的彈性變形應和柔性葉片變形相等［(2)超級單元應與相同尺寸的剛性梁有相同的質量和慣性性質［(3)超級單元的固有頻率應盡量接近連續(xù)梁的固有頻率。

NREL發(fā)布的5 MW近海風力機為上風向變速變槳距風力機，葉片截面剛度分布、質量分布與翼型等參數(shù)及各截面翼型氣動特性等參數(shù)可參考文獻［8］。據(jù)彈性梁的彎曲與扭轉理論［4］，可得各超級單元中彈簧的剛度系數(shù)(表1)。

表1 NREL－5MW風力機葉片彈簧剛度系數(shù)/(Nm/rad)

針對本文葉片多體模型，應用多體動力學軟件ADAMS/Vibration對該葉片的固有頻率進行了計算，結果與NREL提供的基于風力機氣彈分析軟件FAST的分析結果進行對比(表2)，驗證了本模型的有效性。

表2 計算結果對比

2 葉片氣動載荷計算與氣彈耦合分析

運用葉素-動量理論(BEM)計算作用在葉片各剛體上的氣動力，并采用普朗特葉尖修正［5］。假設某剛體Bi(i=1，…，13)各截面處的氣動載荷相同，并始終作用在Bi質心截面的氣動中心處(葉片質心截面弦上距前緣1/4處，如圖3)。則由葉素理論，風輪平面內某一環(huán)形區(qū)域所受的推力FT及轉矩Q計算公式為:

圖3 剛體Bi質心截面上的氣動參數(shù)及氣動中心連體基

由動量理論，風輪平面內某一環(huán)形區(qū)域所受的推力FT及轉矩Q計算方法為:

其中，r為氣動中心距離旋轉中心的距離;B為風輪葉片數(shù);ρ為空氣密度;c為Bi的質心截面處葉素的弦長;U∞為來流風速;ω為風輪旋轉角速度;W為作用在剛體Bi質心截面氣動中心切向速度ωr與來流風速U∞的合速度;α為剛體Bi質心截面處的攻角;L為剛體Bi沿葉片徑向的長度;CL、CD和CM分別是升力系數(shù)、阻力系數(shù)和力矩系數(shù)，它們?yōu)楣ソ铅恋暮瘮?shù);a為軸向誘導因子;a'為切向誘導因子;F為葉尖損失因子。

當葉片變形較大時，考慮葉片變形與揮舞與擺振速度對氣動力的影響，入流角φ及攻角α的計算公式為:

其中，φ為空氣入流角;θ為葉素的扭轉角;Ve-op為氣動中心揮舞速度;Ve-ip為氣動中心擺振速度。

結合以上公式，以及二維翼型升阻特性數(shù)據(jù)，便可通過迭代的方法計算出所需要的物理參數(shù)。具體迭代步驟可參考文獻［11］。

氣彈耦合仿真是在多體系統(tǒng)動力模塊的基礎上增加氣動載荷計算模塊來實現(xiàn)的，耦合仿真流程如圖4所示。仿真過程中如果當前時刻t小于仿真時間tend，則多體動力學模塊將會在每一個時間步中調用氣動模塊計算時變的氣動載荷。這些氣動載荷作為外力輸入到多體系統(tǒng)動力學模塊計算出下一時間步中改變了的葉片各剛體的姿態(tài)坐標(方向余弦陣)、速度和角速度，這些參數(shù)再輸入到氣動模塊，從而可計算出下一時間步的氣動載荷并再次耦合到葉片的結構動力學計算去，如此循環(huán)直到時間終點。

圖4 氣彈耦合分析流程圖

建立約束方程限制葉片轉速

其中q1為葉片相對于機艙的轉動自由度，ωe為額定轉速12.1 r/min。對于變轉速控制，此處亦可采用非線性約束方程進行約束。仍采用違約穩(wěn)定法計算葉片的動力學響應。文中取α =β=10，仿真時間為50 s，積分精度為10-5。

3 算例分析

3.1 模型及其參數(shù)的選取

基于以上的建模基礎，對文獻［8］提供的5 MW風力機葉片在隨機風速下的氣彈響應進行了數(shù)值分析。計算過程在隨機風速及額定轉速(12.1 r/min)下進行。

考慮到此風輪直徑(126 m)較大，葉片上各剛體(葉素)在不同的高度風速存在風剪，引入風速隨高度變化的公式加以修正:

式中，V、V0為距地面高度分別為H及H0處的平均風速，本文中，取V0=11.4 m/s，為輪轂高度H0=90 m處的風速;α為切變系數(shù)或粗糙度指數(shù)，是一個經(jīng)驗指數(shù)，其值為0.125 ～0.6，文中取 0.16。

3.2 葉片的內力分析

葉片變形導致相鄰兩個剛體之間產(chǎn)生相對轉動，因此剛體之間截面上的內力可通過剛體間的力元(彈簧、阻尼器)進行計算，彈簧力和阻力分別與剛體之間的相對轉動角度和轉動速度成正比。

圖5為葉片各剛體間的揮舞力矩圖。從中可以看出，葉片揮舞力矩為正(與風速方向相同)，其中葉根處揮舞力矩最大(B1與B2之間)。揮舞力矩由空氣動力載荷、重力載荷與慣性力所導致，由于風力機在運行過程中，始終受到隨機來流的壓迫作用;同時，各剛體的加速度也隨時間不斷變化，造成了揮舞力矩的波動。

圖5 葉片各剛體揮舞力矩

圖6為各剛體間的擺振力矩。從中可以看到，在風力機運行過程中擺振力矩一般為負(與風輪旋轉方向相同)，這是由于葉片始終受到升力的作用。而在這一過程中，葉片變形產(chǎn)生的結構彈性力與氣動力相互影響，導致葉片的擺振力矩作周期性的波動。

圖7為葉片超級單元(共4個)內的扭轉力矩。從中可以看到，在風力機運行過程中扭轉力矩一般為負，這是由于葉片各剛體氣動中心偏離扭轉軸，而剛體始終受到來流風所引起的空氣阻力的作用。

圖8為葉尖的運動軌跡。如圖8所示，由于葉片的預彎作用，故葉尖的初始點不在零點。由圖8可以看出，當風力機工作在工況下時，葉尖在氣動力及結構力的相互作用下，始終在一定范圍內不斷顫動。

4 結束語

本文應用超級單元法將柔性葉片離散為有限個體，并應用R-W建模方法，建立柔性葉片的非線性代數(shù)-微分氣彈耦合方程，結合BEM計算葉片的空氣動力和適當?shù)臄?shù)值求解方法，用不多的自由度，實現(xiàn)了大型風力機葉片的非線性氣彈耦合時域響應分析。算例對5 MW的風力機在紊流風速下的非線性氣彈響應進行了數(shù)值模擬。

圖6 葉片各剛體擺振力矩

圖7 葉片各超級單元內的扭轉力矩

圖8 葉尖在垂直于葉片徑向的平面內的運動軌跡

柔性葉片非線性多體動力學氣彈耦合方程的構建與數(shù)值求解可通過編制的軟件由計算機完成，大大節(jié)省建模與求解方程的時間以及可能引起的差錯。通過建立不同的約束條件，可以實現(xiàn)變速或定轉速風力機的時域分析。應用該方法，可以分析葉片上各位置處的截面內力，進而進行葉片的強度、剛度分析和氣彈穩(wěn)定性分析。通過時域響應的頻域分析，也可以實現(xiàn)葉片動力特性如固有頻率、結構阻尼與氣彈阻尼等與氣彈穩(wěn)定性相關的分析，分析的精度可通過增加超級單元的個數(shù)進行調整。在葉片設計或優(yōu)化階段，僅需通過輸入文件，調整輸入?yún)?shù)就能得到不同的優(yōu)化結果。

基于時域響應的大型風力機葉片氣彈耦合分析，對于大型風力機葉片的優(yōu)化設計和安全穩(wěn)定運行具有重要的意義，也是進行葉片氣彈穩(wěn)定性分析的重要基礎。本文的工作對于葉片強度與剛度校核與優(yōu)化設計有重要的應用價值。

［1］ Hansen M O L，S?rensen J N，Voutsinas S，et al.State of the art in wind turbine aerodynamics and aeroelasticity［J］.Progress in Aerospace Sciences，2006，42(4):285-330.

［2］ Kallesoe B S.Effect of steady deflections on the aeroelastic stability of a turbine blade［J］.Wind Energy，2011，14:209-224.

［3］ Molenaar D P.Cost-effective design and operation of variable speed wind turbines［D］.Delft:Delft University of Technology，2003.

［4］ Holierhoek J G.Aeroelasticity of Large Wind Turbines［D］.Delft:Delft University of Technology，2008.

［5］ Hansen M O L.Aerodynamics of wind turbines［M］.2nd ed.UK:Earthscan，2008.

［6］ Kim S，Sclavounos P D.Fully coupled response simulations of theme offshore structure in water depths of up to 10000 feet［C］//Proceedings of 11th International Offshore and Polar Engineering Conference(ISOPE)，Stavanger，Norway，June 17-22，2001:457-466.

［7］ Meng F Z.Aero-elastic Stability analysis for largescale wind turbines［D］.Delft:Delft University of Technology，2011.

［8］ Jonkman J，Butterfield S，Musial W，et al.Definition of a 5-MW reference wind turbine for offshore system development［R］.Springfield:U.S.National Renewable Energy Laboratory，2009.

［9］洪嘉振.計算多體系統(tǒng)動力學［M］.北京:高等教育出版社，1999.

［10］ Gebhardt C G，Preidikman S，J?rgensen M H，et al.Non-linear aeroelastic behavior of large horizontalaxis wind turbines:A multibody system approach［J］.International Journal of Hydrogen Energy 2012，37(19):14719-14724.

［11］ 李 春，葉 舟，高 偉，等.現(xiàn)代陸海風力機計算與仿真［M］.上海:上海科學技術出版社，2012.