李治國，陳 猛，張雅靜，高志鷹，汪建文

（1.內(nèi)蒙古工業(yè)大學(xué) 機(jī)械工程學(xué)院，內(nèi)蒙古 呼和浩特 010051；2.內(nèi)蒙古工業(yè)大學(xué) 能源與動(dòng)力工程學(xué)院，內(nèi)蒙古 呼和浩特 010051；3.內(nèi)蒙古建筑職業(yè)技術(shù)學(xué)院 交通與市政工程學(xué)院，內(nèi)蒙古 呼和浩特 010070；4.內(nèi)蒙古工業(yè)大學(xué) 風(fēng)能太陽能利用技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，內(nèi)蒙古 呼和浩特 010051）

0 引言

風(fēng)力機(jī)葉片氣動(dòng)彈性問題大多是揮舞、擺振以及扭轉(zhuǎn)之間的耦合，可形成破壞力極強(qiáng)的經(jīng)典顫振，是葉片設(shè)計(jì)時(shí)必須考慮的重點(diǎn)問題[1]。當(dāng)前風(fēng)力機(jī)逐漸向大型化發(fā)展，葉片越發(fā)細(xì)長，抗彎曲能力逐漸減弱，使葉片容易彎折破壞[2]。由于葉片扭轉(zhuǎn)剛度較大并且氣動(dòng)力產(chǎn)生的扭矩較小，所以揮舞變形與擺振變形是葉片容易彎折破壞的主要因素。研究柔性葉片氣彈穩(wěn)定性對(duì)保護(hù)葉片，防止葉片的彎折破壞尤為重要。

文獻(xiàn)[3]～[5]基于ONERA氣動(dòng)力模型對(duì)翼型的動(dòng)態(tài)失速進(jìn)行了研究，計(jì)算出了失速翼型的顫振臨界速度，得到了氣動(dòng)力模型的適用范圍和在深失速等具有明顯多級(jí)失速現(xiàn)象下的計(jì)算精度。文獻(xiàn)[6]～[8]利用ONERA非定常氣動(dòng)力模型，對(duì)翼型進(jìn)行非線性氣動(dòng)彈性分析，研究非定常氣動(dòng)力對(duì)翼型的影響并建立了非定常氣彈動(dòng)力學(xué)模型。常林[9]模擬了極端條件下葉片的揮舞和擺振方向受迫振動(dòng)。Chaviaropoulos P K[10]基于常值變槳角和典型截面，研究了揮舞-扭轉(zhuǎn)、揮舞-擺振的顫振問題。楊浩南[11]采用修正動(dòng)量葉素理論（BEM）和幾何精確梁理論（GEBT）耦合的葉片氣彈研究方法，研究了不同來流條件下葉片和塔筒的氣彈變形及其對(duì)風(fēng)力機(jī)氣動(dòng)性能和載荷的影響。黃鑫祥[12]基于風(fēng)力機(jī)整機(jī)剛?cè)狁詈夏Ｐ?，提出了一種葉片動(dòng)態(tài)氣彈扭轉(zhuǎn)變形分析的新方法，并分析葉片氣彈扭轉(zhuǎn)變形對(duì)風(fēng)力機(jī)氣動(dòng)功率及氣彈穩(wěn)定性的影響。

研究人員已從不同角度對(duì)葉片氣彈穩(wěn)定性問題進(jìn)行了研究，但對(duì)葉片變槳前后的氣彈穩(wěn)定性研究相對(duì)較少。為了探究大型風(fēng)力機(jī)柔性葉片在揮舞-擺振耦合作用下變槳前后的氣彈穩(wěn)定性，本文基于ONERA非線性氣動(dòng)模型建立包括二維翼型非線性氣動(dòng)升阻力方程及其揮舞-擺振耦合的簡化運(yùn)動(dòng)方程的氣彈模型。利用該模型計(jì)算美國國家能源實(shí)驗(yàn)室（NREL）5 MW風(fēng)力機(jī)葉片中段的DU35-A17翼型變槳前后揮舞、擺振變形量，并與開源軟件FAST計(jì)算得到的結(jié)果進(jìn)行對(duì)比分析，驗(yàn)證本文建立的氣彈穩(wěn)定模型的正確性。

1 葉片氣彈方程

1.1 柔性葉片揮舞-擺振耦合運(yùn)動(dòng)方程

ONERA模型以微分方程組的形式給出了氣動(dòng)載荷的表達(dá)式[13]。旋轉(zhuǎn)葉片揮舞-擺振耦合的簡化運(yùn)動(dòng)方程為[14]

圖1所示為風(fēng)力機(jī)葉片截面結(jié)構(gòu)與坐標(biāo)系示意圖。圖中，α為攻角，θ為扭角，θ0為槳距角，φ為相對(duì)風(fēng)速角，U為風(fēng)速，V為相對(duì)風(fēng)速，Ω為葉片旋轉(zhuǎn)角速度，r為旋轉(zhuǎn)半徑。

圖1 風(fēng)力機(jī)葉片截面坐標(biāo)示意圖Fig.1 Cross section coordinate diagram of wind turbine

1.2 ONERA氣動(dòng)力模型

本文采用ONERA非線性氣動(dòng)力模型[15]。

1.2.1 非線性升力系數(shù)

通過式（2）得到葉片截面的非線性升力系數(shù)。

式中：CLa為線性氣動(dòng)力部分對(duì)應(yīng)的與升力有關(guān)的氣動(dòng)力系數(shù)；CLb為非線性氣動(dòng)力部分對(duì)應(yīng)的與升力有關(guān)的氣動(dòng)力系數(shù)；（＇）為關(guān)于時(shí)間或攻角的一階導(dǎo)數(shù)；（″）為關(guān)于時(shí)間或攻角的二階導(dǎo)數(shù)；CLy為與升力有關(guān)的空氣動(dòng)力系數(shù)；ΔCL為靜態(tài)空氣動(dòng)力系數(shù)曲線線性部分的延長線與曲線非線性部分的差；SZ1，SZ2，SZ3，α0z，λ1，λ2均為與線性氣動(dòng)力部分有關(guān)的參數(shù)；r1z，r2z，r3z均為與雷諾數(shù)有關(guān)的參數(shù)。

ONERA模型中的線性部分是對(duì)經(jīng)典的Theodorsen氣動(dòng)力的模擬，因此，ONERA模型中與線性氣動(dòng)力部分有關(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。

為了方便計(jì)算，通常將靜態(tài)氣動(dòng)力曲線簡化為折線的形式（圖2）。

圖2 靜態(tài)氣動(dòng)力曲線簡化圖Fig.2 Simplified diagram of static aerodynamic curve

由圖2可以近似得到ΔCL的表達(dá)式。

1.2.2 非線性阻力系數(shù)

通過下式可以計(jì)算出葉片截面的阻力系數(shù)。

式中：CDa，CDb為阻力系數(shù)的組成部分；ΔCD為線性部分阻力系數(shù)與非線性部分阻力系數(shù)的差值；r1D，r2D，r3D均為與阻力系數(shù)有關(guān)的參數(shù)。

式中各參數(shù)的值為

最后通過將上述公式聯(lián)立解微分方程組即可求出葉片的揮舞、擺振變形量。本文選用四階龍格-庫塔法計(jì)算微分方程組。

式中：設(shè)X1為z；X2為z＇；X3為y；X4為y＇；X5為CLy；X6為CLb；X7為CLb＇；X8為CDb；X9為CDb＇。

ONERA氣動(dòng)力模型計(jì)算流程如圖3所示。

圖3 ONERA氣動(dòng)力模型計(jì)算流程圖Fig.3 Calculation flow chart of ONERA aerodynamic model

2 算例分析

本文選用NREL 5 MW風(fēng)力機(jī)葉片中的Du35-A17翼型進(jìn)行計(jì)算，此翼型弦長為4.652 m，扭角為11.48°，截面的旋轉(zhuǎn)半徑r為15.85 m，額定風(fēng)速為11.4 m/s。風(fēng)力機(jī)風(fēng)速與槳距角變化關(guān)系見表1。

表1 不同風(fēng)速下風(fēng)力機(jī)槳距角Table 1 Pitch angle of the wind turbine under different wind speeds

2.1 氣彈模型驗(yàn)證

為了驗(yàn)證本文氣彈模型的準(zhǔn)確性，對(duì)NREL 5 MW風(fēng)力機(jī)在額定工況下的葉尖變形量進(jìn)行計(jì)算并與NREL的FAST軟件進(jìn)行對(duì)比，兩種方法計(jì)算結(jié)果如圖4所示。

由圖4可知：FAST計(jì)算的揮舞變形量為5.945 8 m，擺振變形量為-1.108 3×10-1～-1.528 7×10-1m，均值為-6.305 7×10-1m；ONERA模型計(jì)算的舞變形量為5.923 0 m，較FAST計(jì)算結(jié)果小0.383 5%，擺振變形量為-1.044 6×10-1～-2.802 6×10-1m，均值為-6.624 3×10-1m，較FAST計(jì)算結(jié)果大5.0526%。ONERA模型計(jì)算結(jié)果穩(wěn)定后與FAST計(jì)算結(jié)果基本吻合，可以認(rèn)為本文模型能夠用于分析葉片氣彈穩(wěn)定性。

圖4 額定工況下兩種方式計(jì)算的葉尖變形量Fig.4 Calculation of tip deflection in two ways under rated working conditions

2.2 未變槳時(shí)葉片揮舞擺振變形量

利用ONERA模型對(duì)葉片運(yùn)動(dòng)進(jìn)行分析，分別計(jì)算了當(dāng)風(fēng)速為5，8 m/s和10 m/s時(shí)葉片中段（r=15.85 m）的揮舞變形量與擺振變形量。此3種風(fēng)速皆在風(fēng)力機(jī)的設(shè)計(jì)工作范圍內(nèi)，所以揮舞變形量與擺振變形量均是收斂的。圖5，6分別為未變槳時(shí)揮舞變形圖和擺振變形圖。

圖5 未變槳時(shí)揮舞變形圖Fig.5 Flapwise deflection without pitching

圖6 未變槳時(shí)擺振變形圖Fig.6 Edgewise deflection without pitching

由圖5，6可知：揮舞變形量最終會(huì)收斂為一個(gè)固定值，擺振變形量則穩(wěn)定在一個(gè)范圍內(nèi)作規(guī)律的變化；當(dāng)風(fēng)速為5 m/s時(shí)，葉片的揮舞變形量最終穩(wěn)定在2.834×10-2m，擺振變形量穩(wěn)定在-2.788×10-3～-2.312×10-3m；當(dāng)風(fēng)速為8 m/s時(shí)，葉片的揮舞變形量最終會(huì)穩(wěn)定在4.72×10-2m，擺振變形量穩(wěn)定在-7.089×10-3～-6.190×10-3m；當(dāng)風(fēng)速為10 m/s時(shí)，葉片的揮舞變形量最終會(huì)穩(wěn)定在5.108×10-2m，擺振變形量穩(wěn)定在-8.185×10-3～-7.157×10-3m。

隨著風(fēng)速的增大，葉片的揮舞與擺振方向的變形明顯增大。這是由于風(fēng)力機(jī)未變槳，在葉片z向、y向抗彎剛度相同的情況下，隨著風(fēng)速的增加，葉片表面所受風(fēng)載荷增加，使得葉片的變形量增大。相較于風(fēng)速為5 m/s時(shí)，風(fēng)速為10 m/s時(shí)的揮舞變形量增加近一倍，擺振變形量增加近兩倍。這是由于風(fēng)速增加一倍，葉片所受的風(fēng)載荷相應(yīng)的增加一倍，揮舞變形量主要受風(fēng)載荷的影響，所以增加近一倍，而擺振變形量不僅受風(fēng)載荷的影響，還受葉片旋轉(zhuǎn)的影響，風(fēng)輪的轉(zhuǎn)速會(huì)隨著風(fēng)速加快而加快，所以擺振變形量增加近兩倍。

2.3 變槳后葉片揮舞擺振變形量

通過查找NREL 5 MW風(fēng)力機(jī)手冊可知，風(fēng)力機(jī)的額定風(fēng)速為11.4 m/s，額定轉(zhuǎn)速為12.1 r/min，當(dāng)風(fēng)速為13 m/s時(shí)，葉片槳距角為6.60°，當(dāng)風(fēng)速為15 m/s時(shí)，葉片槳距角為10.45°。利用ONERA模型對(duì)葉片運(yùn)動(dòng)進(jìn)行分析，分別計(jì)算當(dāng)風(fēng)速為11.4，13 m/s和15 m/s時(shí)葉片的揮舞變形量與擺振變形量。由于葉片的變槳，此3種風(fēng)速仍在風(fēng)力機(jī)的設(shè)計(jì)工作范圍內(nèi)，所以揮舞變形量與擺振變形量均是收斂的。圖7，8分別為葉片變槳后揮舞變形圖和擺振變形圖。

圖7 變槳后揮舞變形圖Fig.7 Flapwise deflection after pitching

圖8 變槳后擺振變形圖Fig.8 Edgewise deflection after pitching

由圖7，8可知：揮舞變形量最終會(huì)收斂為一個(gè)固定值，擺振變形量最終會(huì)穩(wěn)定在一個(gè)范圍內(nèi)作規(guī)律的變化；當(dāng)風(fēng)速為11.4 m/s時(shí)，葉片的揮舞變形量最終穩(wěn)定在5.702×10-2m，擺振變形量穩(wěn)定在-1.006×10-2～-8.800×10-3m；當(dāng)風(fēng)速為13 m/s時(shí)，葉片的揮舞變形量最終會(huì)穩(wěn)定在5.035×10-2m，擺振變形量穩(wěn)定在-7.195×10-3～-6.198×10-3m；當(dāng)風(fēng)速為15 m/s時(shí)，葉片的揮舞變形量最終會(huì)穩(wěn)定在5.537×10-2m，擺振變形量穩(wěn)定在-7.905×10-3～-6.763×10-3m。當(dāng)風(fēng)力機(jī)變槳之后，葉片的揮舞、擺振變形量比額定風(fēng)速時(shí)的揮舞擺振變形量有所減小，這是由于隨著風(fēng)速的增加，風(fēng)力機(jī)要維持在額定功率下運(yùn)動(dòng)需要變槳，增大槳距角，槳距角的增大會(huì)改變?nèi)~片的抗彎剛度，葉片抵抗變形的能力增大，使得揮舞變形量與擺振變形量相對(duì)于變槳之前有所下降。而隨著風(fēng)速從13～15 m/s的進(jìn)一步增大，葉片繼續(xù)變槳，槳距角增大，葉片所受風(fēng)載荷的影響超過抗彎剛度增大的影響，葉片揮舞變形量與擺振變形量會(huì)繼續(xù)隨風(fēng)速的增大而增大。

表2為葉片揮舞、擺振變形量計(jì)算結(jié)果。

表2 揮舞、擺振變形量計(jì)算結(jié)果Table 2 Calculation results of flapwise and edgewise deflection

3 結(jié)論

本文基于ONERA非線性氣動(dòng)模型，建立了包括二維翼型非線性氣動(dòng)升阻力方程及其揮舞-擺振耦合運(yùn)動(dòng)方程的柔性葉片氣彈模型，計(jì)算并獲得葉片中段變槳前后的揮舞、擺振變形量變化曲線，得到了風(fēng)力機(jī)葉片在未變槳與變槳后，隨著風(fēng)速的增大，揮舞、擺振變形量的變化規(guī)律。

①在額定工況時(shí)，葉片通過向z軸正方向、y軸負(fù)方向的彎曲以抵消風(fēng)載荷對(duì)葉片的影響。揮舞變形量逐漸增大到一定值后，最終會(huì)穩(wěn)定在5.702×10-2m，擺振位移逐漸增大，最終會(huì)穩(wěn)定在-1.006×10-2～-8.800×10-3m。

②風(fēng)力機(jī)未變槳時(shí)，由于葉片抗彎剛度相同，揮舞、擺振變形量會(huì)隨風(fēng)速增大而增大，其中在風(fēng)載荷增加一倍的情況下，揮舞變形量會(huì)隨之增加近一倍達(dá)到2.274×10-2m，擺振變形量由于風(fēng)載荷與葉片旋轉(zhuǎn)的共同作用增加約兩倍達(dá)到5.121×10-3m。

③葉片變槳后，由于葉片抗彎剛度的增加，揮舞、擺振變形量比額定工況下的變形量有所減少。隨著來流風(fēng)速增加，風(fēng)載荷的影響超過抗彎剛度增大的影響，葉片揮舞變形量與擺振變形量會(huì)繼續(xù)隨風(fēng)速的增大而增大。不同風(fēng)速下的揮舞、擺振變形量總是收斂于某一個(gè)固定值，故葉片是氣彈穩(wěn)定的。