基于ANSYS的10 MW風(fēng)力機葉片彎扭耦合特性研究

VasilisRiziotisGiannisSeferainTohidBagherpour

(1.華中科技大學(xué)中歐清潔與可再生能源學(xué)院，武漢 430074； 2.華中科技大學(xué)能源與動力工程學(xué)院，武漢 430074； 3.School of Mechanical Engineering, National Technology University of Athens, 9 HeroonPolytechneiou Street, GR15780 Athens, Greece)

0 引 言

隨著風(fēng)力機發(fā)電容量的不斷增大，風(fēng)力機葉片趨于大型化的同時，承受的載荷也越來越大。目前主流葉片都采用復(fù)合材料制成。在葉片中引入材料的彎扭耦合(BTC)效應(yīng)成為一種流行的葉片減載荷方式[1]。彎扭耦合是指當(dāng)葉片在風(fēng)載荷作用下，彎曲的同時發(fā)生扭轉(zhuǎn)變形的現(xiàn)象。這種現(xiàn)象會導(dǎo)致攻角的改變，減少葉片承受的載荷。該方法是通過改變?nèi)~片梁帽上或者表皮上各向異性材料的鋪層角度來實現(xiàn)鋪層水平上法向應(yīng)力和剪切應(yīng)力的耦合，從而引起葉片的彎曲和扭轉(zhuǎn)耦合。

對葉片彎扭耦合的模擬是基于不均勻各向異性梁理論和有限元理論。Blasques等[2]結(jié)合Saint-Venant理論開發(fā)了二維梁截面特性分析工具BECAS。Bagherpour等[3]采用BECAS結(jié)合hGAST(一維各向異性梁分析工具)模擬了葉片的彎扭耦合特性。Fedorov和Berggreen[4]提出彎扭耦合系數(shù)計算方法。Hang Meng[5]等使用BECAS和一維各向異性梁模擬在彎扭耦合影響下葉片承受的載荷用于疲勞壽命分析。上述文獻都是將3D葉片分為了2D翼型截面和1D梁結(jié)合的形式來分析彎扭耦合特性。這種簡化的方式將模型置于不同的計算環(huán)境，且網(wǎng)格劃分稀疏，增大了結(jié)果的不確定性。此外當(dāng)需要通過改變材料鋪層角度在葉片中引入彎扭耦合特性時，計算過程繁瑣。文中通過建立葉片整體三維有限元模型分析彎扭耦合特性?？捎行П苊庥捎谟嬎悱h(huán)境不同帶來的誤差，且網(wǎng)格劃分密集，增加了結(jié)果的可靠性。另外直接通過改變參數(shù)化編程語言中參數(shù)的方式在葉片中引入不同的鋪層角度和鋪層位置，極大的提高了計算效率。國內(nèi)在風(fēng)力機葉片彎扭耦合領(lǐng)域的文獻較少，隨著海上風(fēng)電的發(fā)展，大容量風(fēng)力機組已成必然趨勢。劉宇航等[6]和安立強[7]等在ANSYS中建立了三維葉片模型，分析了葉片在氣動載荷作用下的變形情況。但由于采用的5 MW 葉片數(shù)據(jù)，翼型截面結(jié)構(gòu)簡單，結(jié)果不理想。文中采用了最新的DTU 10 MW風(fēng)力機葉片數(shù)據(jù)，分析葉片在不同鋪層角度下的氣動彈性響應(yīng)和該葉片能夠取得的最大扭轉(zhuǎn)變形，彎扭耦合效應(yīng)明顯。為其他研究人員提供參考。

1 有限元模型

1.1 箱形梁有限元模型

葉片具有長寬比，長厚比大的特點，常用箱形梁做初步模擬。文中首先在ANSYS中建立了箱形梁有限元模型，與已有的實驗數(shù)據(jù)和其他數(shù)值模擬結(jié)果作對比分析。

采用了Chadra[8]實驗梁數(shù)據(jù)。其為薄壁復(fù)合材料梁，長0.762 m，恒定的箱型截面，尺寸為13.6*24.2 mm(0.537*0.953 inch ),壁厚7.6 mm(0.03 inch)，對稱鋪層，上下為[45]6，左右為[±45]3,材料為AS4/3501-6單向石墨/環(huán)氧樹脂，經(jīng)坐標系轉(zhuǎn)換后完整的材料特性見表1。

表1箱型梁材料特性

材料名稱彈性模量/GPa泊松比ExEyEzGxyGxzGyzVxyVxzVyz密度/kg·m-3AS4/3501-61420.810.8166.133.770.30.420.421580

邊界條件為一端采用固定約束，另一端加載垂直集中力4.448 N。在ANSYS中選擇適用于復(fù)合材料建模的shell181單元，它是每個節(jié)點具有6個自由度的4節(jié)點單元，適用于模擬分層的復(fù)合殼或夾層結(jié)構(gòu)。此外需要注意的是，鋪層的偏置情況需選擇向內(nèi)偏置。在箱型梁各截面建立中心節(jié)點，并通過MPC多點約束方法耦合參考點與截面上各從節(jié)點來分析梁的變形情況。

圖1展示了ANSYS結(jié)果與Chadra的實驗解，Stablein和Hansen[9]，Smith和Chopra[10]的數(shù)值解以及Bagherpour使用一維梁模型和二維截面結(jié)合模擬得到的數(shù)值結(jié)果的對比。在扭轉(zhuǎn)角度上ANSYS很好的捕捉了變化趨勢；在彎曲方向上，ANSYS表現(xiàn)出良好的預(yù)測能力，相比于其他軟件，ANSYS的結(jié)果更接近實驗解。

圖1 梁模型的實驗解和數(shù)值解

1.2 DTU 10MW 實驗風(fēng)機

該風(fēng)力機設(shè)計風(fēng)況等級為IEC IA級，為變速變槳距，三葉片風(fēng)力機，風(fēng)輪直徑為178.3 m，輪轂中心距地面119 m。切入風(fēng)速為4 m/s，切出風(fēng)速為25 m/s，額定風(fēng)速為11.4 m/s，在額定轉(zhuǎn)速9.6rpm下達到10 MW額定輸出功率。葉片總長度為89.166 m，從葉根到葉尖共分為101個不同的翼型截面，3個腹板支撐結(jié)構(gòu)。如圖2所示每個翼型截面分為11個部分。

圖2 翼型截面示意

葉片由單軸向布，雙軸向布，三軸向布和巴沙木四種材料組成，材料特性如表2所示。材料的鋪層順序和厚度以及葉片的預(yù)彎和預(yù)扭隨葉片展向分布可參見文獻[11]。關(guān)于DTU10 MW風(fēng)力機的更多參數(shù)可參見文獻[12]。翼型截面坐標系為X軸沿弦線方向由尾緣指向前緣；Y軸垂直于X軸由吸力面指向壓力面，Z軸由葉片根部指向葉尖。載荷計算也將遵循該坐標系。

1.3 建立葉片有限元模型

通過設(shè)置材料屬性-輸入關(guān)鍵點-生成翼型輪廓-生成面-定義鋪層-網(wǎng)格劃分的步驟建立葉片殼有限元模型。利用參數(shù)化設(shè)計語言APDL的循環(huán)特性實現(xiàn)快速建模。葉片模型共由19569個單元和19007個節(jié)點組成，其中包括101個在不同翼型截面位置處建立的MPC主節(jié)點。

在葉根加固定約束后，前八階振型如圖3所示，葉片振型從低階到高階逐漸復(fù)雜，一階振型為葉片一階揮舞振動，二階振型為葉片一階擺振振動，三階振型為二階揮舞振動，四階振型為二階擺振振動，五階振型為三階揮舞振動，六階振型為一階扭轉(zhuǎn)振動，七階振型為高階擺振振動，八階振型為高階揮舞振動。

圖3 葉片前八階振型

對應(yīng)的固有頻率與其他模擬空氣彈性變形工具的結(jié)果[13]作對比，見表3。其中hGAST，NEREA，Cp-Lambda，HAWC2采用的是梁理論，NISA是三維有限元模型?？梢夾NSYS的模擬結(jié)果與其他軟件結(jié)果具有良好的一致性。

表2葉片材料特性

材料類型楊氏模量(Gpa)剪切模量(Gpa)泊松比ExEyEzGxGyGzVyxVyzVxz密度(kg/m3)單軸向布14.9313.4241.635.005.005.000.030.260.241915.50雙軸向布13.9212.1013.9211.504.534.530.030.270.531845.00三軸向布14.7012.1021.809.414.534.530.030.270.471845.00巴沙木0.052.730.050.020.150.150.010.010.50110.00

表3葉片固有頻率對比

ModeANSYShGASTNEREACp-LambdaHAWC2NISA11st flap0.601 430.620.620.620.610.6421st edge0.934 580.940.940.940.930.9632nd flap1.722 91.761.741.761.741.8542nd edge2.807 82.82.792.82.772.8653rd flap3.518 63.593.523.63.573.7661st torsion5.486 15.45.36-6.66.0173rd edge5.587 95.735.615.745.75.8284th flap6.067 26.096.036.116.11-

2 葉片的彎扭耦合特性

通過改變梁帽位置單軸向布的鋪層角度在葉片結(jié)構(gòu)中引入彎扭耦合。鋪層角度在葉片上下表面對稱分布。計算葉片在11.4 m/s額定風(fēng)速下，即風(fēng)力機正常發(fā)電工況下葉片的彎曲和扭轉(zhuǎn)變形情況，其中風(fēng)為穩(wěn)態(tài)風(fēng)。

通過葉素動量理論(BEM)和葉片變形量迭代得到葉片受到的氣動載荷。過程如下：首先在ANSYS中輸入一個初始載荷，之后將對應(yīng)的葉片變形輸入BEM程序得到新的載荷，重復(fù)該過程直到葉片變形穩(wěn)定，得到該葉片變形下對應(yīng)的載荷。鋪層角度為0度和10度條件下葉片受到的氣動載荷設(shè)為載荷1和載荷2，如圖4所示。氣動載荷在葉根區(qū)域較小，沿著葉展方向逐漸增大，接近葉尖位置時急速下降。載荷1和載荷2均為均布載荷，經(jīng)過式1轉(zhuǎn)化為集中力后通過MPC節(jié)點耦合的方式施加到葉片殼有限元模型的101個截面上。

設(shè)截面i距葉根距離為X，受到的集中力為F，則：

(1)

其中，(0

圖4 額定風(fēng)況下氣動載荷

圖5展示了葉片在0度鋪層下的變形情況及10度鋪層下加載載荷1和載荷2得到的彎扭變形。葉片的彎扭變形連續(xù)，沿著葉展方向不斷增大，在葉尖達到最大值。在額定風(fēng)速下，0度鋪層葉片葉尖彎曲變形為8.89m，扭轉(zhuǎn)變形為1.06度。10度鋪層葉片在載荷2作用下，葉尖彎曲變形為9.68m，扭轉(zhuǎn)變形為4.16度。10度鋪層葉片在兩種情況下葉尖最大扭轉(zhuǎn)角度誤差為3.6%。由于計算載荷的迭代過程復(fù)雜，接下來的計算將直接采用在0度鋪層條件下的載荷。

圖5 葉片在0度鋪層和10度鋪層下變形

由于材料鋪層角度和耦合區(qū)域的不同會影響葉片的彎扭耦合性能。這里設(shè)定葉尖扭轉(zhuǎn)角度在2度到6度以0.5度的間隔變化，通過改變鋪層角度和起始位置的不同組合以得到預(yù)期的葉尖扭轉(zhuǎn)角度。選定鋪層起始位置為0%，20%，35%或50%，結(jié)果如圖6所示。從圖中可以看到耦合區(qū)域越靠近葉尖，達到指定的葉尖扭轉(zhuǎn)角度需要的鋪層角度越大；彎扭耦合系數(shù)在鋪層角度為18.5度之后下降；可達到的最大葉尖扭轉(zhuǎn)變形為6度。這表明如果繼續(xù)增大鋪層角度，不會造成葉片扭轉(zhuǎn)變形的繼續(xù)增大，從而不會造成葉片載荷的降低。為設(shè)計具有優(yōu)良彎扭耦合性能的葉片提供參考。

圖6 彎扭耦合特性圖

3 結(jié)束語

(1)在ANSYS里選用shell 181單元建立箱形梁模型，改變梁表面材料鋪層角度，梁在自由端集中力作用下的變形與已有的實驗結(jié)果和數(shù)值結(jié)果基本一致。在彎曲變形上ANSYS模擬結(jié)果更接近實驗結(jié)果。

(2)建立了10 MW葉片殼有限元模型，通過改變梁帽位置材料鋪層角度在葉片中引入彎扭耦合效應(yīng)。由于采用參數(shù)化設(shè)計語言，在快速實現(xiàn)模型建立之外，還極大的方便了模型參數(shù)的修改，例如鋪層角度和耦合位置的改變，極大的提高了計算效率。

(3)通過繪制彎扭耦合特性圖發(fā)現(xiàn)葉片的扭轉(zhuǎn)變形不會隨著鋪層角度的增大而無限增大，為設(shè)計具有優(yōu)良氣動彈性裁剪性能的葉片提供參考。

(4)文中模擬了葉片在彎扭耦合作用下的變形情況，沒有涉及葉片載荷量的計算，接下來可在這一方面繼續(xù)研究。此外，由于載荷的降低，在葉片疲勞失效范圍內(nèi)可以考慮減少葉片材料，達到經(jīng)濟性的目的。