李治國,陳 猛,張雅靜,高志鷹,汪建文
(1.內(nèi)蒙古工業(yè)大學(xué) 機械工程學(xué)院,內(nèi)蒙古 呼和浩特 010051;2.內(nèi)蒙古工業(yè)大學(xué) 能源與動力工程學(xué)院,內(nèi)蒙古 呼和浩特 010051;3.內(nèi)蒙古建筑職業(yè)技術(shù)學(xué)院 交通與市政工程學(xué)院,內(nèi)蒙古 呼和浩特 010070;4.內(nèi)蒙古工業(yè)大學(xué) 風(fēng)能太陽能利用技術(shù)教育部重點實驗室,內(nèi)蒙古 呼和浩特 010051)
風(fēng)力機葉片氣動彈性問題大多是揮舞、擺振以及扭轉(zhuǎn)之間的耦合,可形成破壞力極強的經(jīng)典顫振,是葉片設(shè)計時必須考慮的重點問題[1]。當(dāng)前風(fēng)力機逐漸向大型化發(fā)展,葉片越發(fā)細(xì)長,抗彎曲能力逐漸減弱,使葉片容易彎折破壞[2]。由于葉片扭轉(zhuǎn)剛度較大并且氣動力產(chǎn)生的扭矩較小,所以揮舞變形與擺振變形是葉片容易彎折破壞的主要因素。研究柔性葉片氣彈穩(wěn)定性對保護葉片,防止葉片的彎折破壞尤為重要。
文獻[3]~[5]基于ONERA氣動力模型對翼型的動態(tài)失速進行了研究,計算出了失速翼型的顫振臨界速度,得到了氣動力模型的適用范圍和在深失速等具有明顯多級失速現(xiàn)象下的計算精度。文獻[6]~[8]利用ONERA非定常氣動力模型,對翼型進行非線性氣動彈性分析,研究非定常氣動力對翼型的影響并建立了非定常氣彈動力學(xué)模型。常林[9]模擬了極端條件下葉片的揮舞和擺振方向受迫振動。Chaviaropoulos P K[10]基于常值變槳角和典型截面,研究了揮舞-扭轉(zhuǎn)、揮舞-擺振的顫振問題。楊浩南[11]采用修正動量葉素理論(BEM)和幾何精確梁理論(GEBT)耦合的葉片氣彈研究方法,研究了不同來流條件下葉片和塔筒的氣彈變形及其對風(fēng)力機氣動性能和載荷的影響。黃鑫祥[12]基于風(fēng)力機整機剛?cè)狁詈夏P停岢隽艘环N葉片動態(tài)氣彈扭轉(zhuǎn)變形分析的新方法,并分析葉片氣彈扭轉(zhuǎn)變形對風(fēng)力機氣動功率及氣彈穩(wěn)定性的影響。
研究人員已從不同角度對葉片氣彈穩(wěn)定性問題進行了研究,但對葉片變槳前后的氣彈穩(wěn)定性研究相對較少。為了探究大型風(fēng)力機柔性葉片在揮舞-擺振耦合作用下變槳前后的氣彈穩(wěn)定性,本文基于ONERA非線性氣動模型建立包括二維翼型非線性氣動升阻力方程及其揮舞-擺振耦合的簡化運動方程的氣彈模型。利用該模型計算美國國家能源實驗室(NREL)5 MW風(fēng)力機葉片中段的DU35-A17翼型變槳前后揮舞、擺振變形量,并與開源軟件FAST計算得到的結(jié)果進行對比分析,驗證本文建立的氣彈穩(wěn)定模型的正確性。
ONERA模型以微分方程組的形式給出了氣動載荷的表達式[13]。旋轉(zhuǎn)葉片揮舞-擺振耦合的簡化運動方程為[14]
圖1所示為風(fēng)力機葉片截面結(jié)構(gòu)與坐標(biāo)系示意圖。圖中,α為攻角,θ為扭角,θ0為槳距角,φ為相對風(fēng)速角,U為風(fēng)速,V為相對風(fēng)速,Ω為葉片旋轉(zhuǎn)角速度,r為旋轉(zhuǎn)半徑。
圖1 風(fēng)力機葉片截面坐標(biāo)示意圖Fig.1 Cross section coordinate diagram of wind turbine
本文采用ONERA非線性氣動力模型[15]。
1.2.1 非線性升力系數(shù)
通過式(2)得到葉片截面的非線性升力系數(shù)。
式中:CLa為線性氣動力部分對應(yīng)的與升力有關(guān)的氣動力系數(shù);CLb為非線性氣動力部分對應(yīng)的與升力有關(guān)的氣動力系數(shù);(')為關(guān)于時間或攻角的一階導(dǎo)數(shù);(″)為關(guān)于時間或攻角的二階導(dǎo)數(shù);CLy為與升力有關(guān)的空氣動力系數(shù);ΔCL為靜態(tài)空氣動力系數(shù)曲線線性部分的延長線與曲線非線性部分的差;SZ1,SZ2,SZ3,α0z,λ1,λ2均為與線性氣動力部分有關(guān)的參數(shù);r1z,r2z,r3z均為與雷諾數(shù)有關(guān)的參數(shù)。
ONERA模型中的線性部分是對經(jīng)典的Theodorsen氣動力的模擬,因此,ONERA模型中與線性氣動力部分有關(guān)的參數(shù)SZ1,SZ2,SZ3,α0z,λ1,λ2便可以依據(jù)兩者之間的關(guān)系確定[16]。參數(shù)值分別為SZ1=π,SZ2=π/2,SZ3=0,α0z=5.9,λ1=0.15,λ2=0.55。
為了方便計算,通常將靜態(tài)氣動力曲線簡化為折線的形式(圖2)。
圖2 靜態(tài)氣動力曲線簡化圖Fig.2 Simplified diagram of static aerodynamic curve
由圖2可以近似得到ΔCL的表達式。
1.2.2 非線性阻力系數(shù)
通過下式可以計算出葉片截面的阻力系數(shù)。
式中:CDa,CDb為阻力系數(shù)的組成部分;ΔCD為線性部分阻力系數(shù)與非線性部分阻力系數(shù)的差值;r1D,r2D,r3D均為與阻力系數(shù)有關(guān)的參數(shù)。
式中各參數(shù)的值為
最后通過將上述公式聯(lián)立解微分方程組即可求出葉片的揮舞、擺振變形量。本文選用四階龍格-庫塔法計算微分方程組。
式中:設(shè)X1為z;X2為z';X3為y;X4為y';X5為CLy;X6為CLb;X7為CLb';X8為CDb;X9為CDb'。
ONERA氣動力模型計算流程如圖3所示。
圖3 ONERA氣動力模型計算流程圖Fig.3 Calculation flow chart of ONERA aerodynamic model
本文選用NREL 5 MW風(fēng)力機葉片中的Du35-A17翼型進行計算,此翼型弦長為4.652 m,扭角為11.48°,截面的旋轉(zhuǎn)半徑r為15.85 m,額定風(fēng)速為11.4 m/s。風(fēng)力機風(fēng)速與槳距角變化關(guān)系見表1。
表1 不同風(fēng)速下風(fēng)力機槳距角Table 1 Pitch angle of the wind turbine under different wind speeds
為了驗證本文氣彈模型的準(zhǔn)確性,對NREL 5 MW風(fēng)力機在額定工況下的葉尖變形量進行計算并與NREL的FAST軟件進行對比,兩種方法計算結(jié)果如圖4所示。
由圖4可知:FAST計算的揮舞變形量為5.945 8 m,擺振變形量為-1.108 3×10-1~-1.528 7×10-1m,均值為-6.305 7×10-1m;ONERA模型計算的舞變形量為5.923 0 m,較FAST計算結(jié)果小0.383 5%,擺振變形量為-1.044 6×10-1~-2.802 6×10-1m,均值為-6.624 3×10-1m,較FAST計算結(jié)果大5.0526%。ONERA模型計算結(jié)果穩(wěn)定后與FAST計算結(jié)果基本吻合,可以認(rèn)為本文模型能夠用于分析葉片氣彈穩(wěn)定性。
圖4 額定工況下兩種方式計算的葉尖變形量Fig.4 Calculation of tip deflection in two ways under rated working conditions
利用ONERA模型對葉片運動進行分析,分別計算了當(dāng)風(fēng)速為5,8 m/s和10 m/s時葉片中段(r=15.85 m)的揮舞變形量與擺振變形量。此3種風(fēng)速皆在風(fēng)力機的設(shè)計工作范圍內(nèi),所以揮舞變形量與擺振變形量均是收斂的。圖5,6分別為未變槳時揮舞變形圖和擺振變形圖。
圖5 未變槳時揮舞變形圖Fig.5 Flapwise deflection without pitching
圖6 未變槳時擺振變形圖Fig.6 Edgewise deflection without pitching
由圖5,6可知:揮舞變形量最終會收斂為一個固定值,擺振變形量則穩(wěn)定在一個范圍內(nèi)作規(guī)律的變化;當(dāng)風(fēng)速為5 m/s時,葉片的揮舞變形量最終穩(wěn)定在2.834×10-2m,擺振變形量穩(wěn)定在-2.788×10-3~-2.312×10-3m;當(dāng)風(fēng)速為8 m/s時,葉片的揮舞變形量最終會穩(wěn)定在4.72×10-2m,擺振變形量穩(wěn)定在-7.089×10-3~-6.190×10-3m;當(dāng)風(fēng)速為10 m/s時,葉片的揮舞變形量最終會穩(wěn)定在5.108×10-2m,擺振變形量穩(wěn)定在-8.185×10-3~-7.157×10-3m。
隨著風(fēng)速的增大,葉片的揮舞與擺振方向的變形明顯增大。這是由于風(fēng)力機未變槳,在葉片z向、y向抗彎剛度相同的情況下,隨著風(fēng)速的增加,葉片表面所受風(fēng)載荷增加,使得葉片的變形量增大。相較于風(fēng)速為5 m/s時,風(fēng)速為10 m/s時的揮舞變形量增加近一倍,擺振變形量增加近兩倍。這是由于風(fēng)速增加一倍,葉片所受的風(fēng)載荷相應(yīng)的增加一倍,揮舞變形量主要受風(fēng)載荷的影響,所以增加近一倍,而擺振變形量不僅受風(fēng)載荷的影響,還受葉片旋轉(zhuǎn)的影響,風(fēng)輪的轉(zhuǎn)速會隨著風(fēng)速加快而加快,所以擺振變形量增加近兩倍。
通過查找NREL 5 MW風(fēng)力機手冊可知,風(fēng)力機的額定風(fēng)速為11.4 m/s,額定轉(zhuǎn)速為12.1 r/min,當(dāng)風(fēng)速為13 m/s時,葉片槳距角為6.60°,當(dāng)風(fēng)速為15 m/s時,葉片槳距角為10.45°。利用ONERA模型對葉片運動進行分析,分別計算當(dāng)風(fēng)速為11.4,13 m/s和15 m/s時葉片的揮舞變形量與擺振變形量。由于葉片的變槳,此3種風(fēng)速仍在風(fēng)力機的設(shè)計工作范圍內(nèi),所以揮舞變形量與擺振變形量均是收斂的。圖7,8分別為葉片變槳后揮舞變形圖和擺振變形圖。
圖7 變槳后揮舞變形圖Fig.7 Flapwise deflection after pitching
圖8 變槳后擺振變形圖Fig.8 Edgewise deflection after pitching
由圖7,8可知:揮舞變形量最終會收斂為一個固定值,擺振變形量最終會穩(wěn)定在一個范圍內(nèi)作規(guī)律的變化;當(dāng)風(fēng)速為11.4 m/s時,葉片的揮舞變形量最終穩(wěn)定在5.702×10-2m,擺振變形量穩(wěn)定在-1.006×10-2~-8.800×10-3m;當(dāng)風(fēng)速為13 m/s時,葉片的揮舞變形量最終會穩(wěn)定在5.035×10-2m,擺振變形量穩(wěn)定在-7.195×10-3~-6.198×10-3m;當(dāng)風(fēng)速為15 m/s時,葉片的揮舞變形量最終會穩(wěn)定在5.537×10-2m,擺振變形量穩(wěn)定在-7.905×10-3~-6.763×10-3m。當(dāng)風(fēng)力機變槳之后,葉片的揮舞、擺振變形量比額定風(fēng)速時的揮舞擺振變形量有所減小,這是由于隨著風(fēng)速的增加,風(fēng)力機要維持在額定功率下運動需要變槳,增大槳距角,槳距角的增大會改變?nèi)~片的抗彎剛度,葉片抵抗變形的能力增大,使得揮舞變形量與擺振變形量相對于變槳之前有所下降。而隨著風(fēng)速從13~15 m/s的進一步增大,葉片繼續(xù)變槳,槳距角增大,葉片所受風(fēng)載荷的影響超過抗彎剛度增大的影響,葉片揮舞變形量與擺振變形量會繼續(xù)隨風(fēng)速的增大而增大。
表2為葉片揮舞、擺振變形量計算結(jié)果。
表2 揮舞、擺振變形量計算結(jié)果Table 2 Calculation results of flapwise and edgewise deflection
本文基于ONERA非線性氣動模型,建立了包括二維翼型非線性氣動升阻力方程及其揮舞-擺振耦合運動方程的柔性葉片氣彈模型,計算并獲得葉片中段變槳前后的揮舞、擺振變形量變化曲線,得到了風(fēng)力機葉片在未變槳與變槳后,隨著風(fēng)速的增大,揮舞、擺振變形量的變化規(guī)律。
①在額定工況時,葉片通過向z軸正方向、y軸負(fù)方向的彎曲以抵消風(fēng)載荷對葉片的影響。揮舞變形量逐漸增大到一定值后,最終會穩(wěn)定在5.702×10-2m,擺振位移逐漸增大,最終會穩(wěn)定在-1.006×10-2~-8.800×10-3m。
②風(fēng)力機未變槳時,由于葉片抗彎剛度相同,揮舞、擺振變形量會隨風(fēng)速增大而增大,其中在風(fēng)載荷增加一倍的情況下,揮舞變形量會隨之增加近一倍達到2.274×10-2m,擺振變形量由于風(fēng)載荷與葉片旋轉(zhuǎn)的共同作用增加約兩倍達到5.121×10-3m。
③葉片變槳后,由于葉片抗彎剛度的增加,揮舞、擺振變形量比額定工況下的變形量有所減少。隨著來流風(fēng)速增加,風(fēng)載荷的影響超過抗彎剛度增大的影響,葉片揮舞變形量與擺振變形量會繼續(xù)隨風(fēng)速的增大而增大。不同風(fēng)速下的揮舞、擺振變形量總是收斂于某一個固定值,故葉片是氣彈穩(wěn)定的。