基于ONERA模型的柔性葉片氣彈穩定性研究

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

（1.內蒙古工業大學 機械工程學院，內蒙古 呼和浩特 010051；2.內蒙古工業大學 能源與動力工程學院，內蒙古 呼和浩特 010051；3.內蒙古建筑職業技術學院 交通與市政工程學院，內蒙古 呼和浩特 010070；4.內蒙古工業大學 風能太陽能利用技術教育部重點實驗室，內蒙古 呼和浩特 010051）

0 引言

風力機葉片氣動彈性問題大多是揮舞、擺振以及扭轉之間的耦合，可形成破壞力極強的經典顫振，是葉片設計時必須考慮的重點問題[1]。當前風力機逐漸向大型化發展，葉片越發細長，抗彎曲能力逐漸減弱，使葉片容易彎折破壞[2]。由于葉片扭轉剛度較大并且氣動力產生的扭矩較小，所以揮舞變形與擺振變形是葉片容易彎折破壞的主要因素。研究柔性葉片氣彈穩定性對保護葉片，防止葉片的彎折破壞尤為重要。

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

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

1 葉片氣彈方程

1.1 柔性葉片揮舞-擺振耦合運動方程

ONERA模型以微分方程組的形式給出了氣動載荷的表達式[13]。旋轉葉片揮舞-擺振耦合的簡化運動方程為[14]

圖1所示為風力機葉片截面結構與坐標系示意圖。圖中，α為攻角，θ為扭角，θ0為槳距角，φ為相對風速角，U為風速，V為相對風速，Ω為葉片旋轉角速度，r為旋轉半徑。

圖1 風力機葉片截面坐標示意圖Fig.1 Cross section coordinate diagram of wind turbine

1.2 ONERA氣動力模型

本文采用ONERA非線性氣動力模型[15]。

1.2.1 非線性升力系數

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

式中：CLa為線性氣動力部分對應的與升力有關的氣動力系數；CLb為非線性氣動力部分對應的與升力有關的氣動力系數；（＇）為關于時間或攻角的一階導數；（″）為關于時間或攻角的二階導數；CLy為與升力有關的空氣動力系數；ΔCL為靜態空氣動力系數曲線線性部分的延長線與曲線非線性部分的差；SZ1，SZ2，SZ3，α0z，λ1，λ2均為與線性氣動力部分有關的參數；r1z，r2z，r3z均為與雷諾數有關的參數。

ONERA模型中的線性部分是對經典的Theodorsen氣動力的模擬，因此，ONERA模型中與線性氣動力部分有關的參數SZ1，SZ2，SZ3，α0z，λ1，λ2便可以依據兩者之間的關系確定[16]。參數值分別為SZ1=π，SZ2=π/2，SZ3=0，α0z=5.9，λ1=0.15，λ2=0.55。

為了方便計算，通常將靜態氣動力曲線簡化為折線的形式（圖2）。

圖2 靜態氣動力曲線簡化圖Fig.2 Simplified diagram of static aerodynamic curve

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

1.2.2 非線性阻力系數

通過下式可以計算出葉片截面的阻力系數。

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

式中各參數的值為

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

式中：設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

2 算例分析

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

表1 不同風速下風力機槳距角Table 1 Pitch angle of the wind turbine under different wind speeds

2.1 氣彈模型驗證

為了驗證本文氣彈模型的準確性，對NREL 5 MW風力機在額定工況下的葉尖變形量進行計算并與NREL的FAST軟件進行對比，兩種方法計算結果如圖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計算結果小0.383 5%，擺振變形量為-1.044 6×10-1～-2.802 6×10-1m，均值為-6.624 3×10-1m，較FAST計算結果大5.0526%。ONERA模型計算結果穩定后與FAST計算結果基本吻合，可以認為本文模型能夠用于分析葉片氣彈穩定性。

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

2.2 未變槳時葉片揮舞擺振變形量

利用ONERA模型對葉片運動進行分析，分別計算了當風速為5，8 m/s和10 m/s時葉片中段（r=15.85 m）的揮舞變形量與擺振變形量。此3種風速皆在風力機的設計工作范圍內，所以揮舞變形量與擺振變形量均是收斂的。圖5，6分別為未變槳時揮舞變形圖和擺振變形圖。

圖5 未變槳時揮舞變形圖Fig.5 Flapwise deflection without pitching

圖6 未變槳時擺振變形圖Fig.6 Edgewise deflection without pitching

由圖5，6可知：揮舞變形量最終會收斂為一個固定值，擺振變形量則穩定在一個范圍內作規律的變化；當風速為5 m/s時，葉片的揮舞變形量最終穩定在2.834×10-2m，擺振變形量穩定在-2.788×10-3～-2.312×10-3m；當風速為8 m/s時，葉片的揮舞變形量最終會穩定在4.72×10-2m，擺振變形量穩定在-7.089×10-3～-6.190×10-3m；當風速為10 m/s時，葉片的揮舞變形量最終會穩定在5.108×10-2m，擺振變形量穩定在-8.185×10-3～-7.157×10-3m。

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

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

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

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

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

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

表2為葉片揮舞、擺振變形量計算結果。

表2 揮舞、擺振變形量計算結果Table 2 Calculation results of flapwise and edgewise deflection

3 結論

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

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

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

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