考慮轉捩影響的風力機翼型氣動特性計算研究

侯銀珠 ,宋文萍 ,張 坤

(西北工業大學翼型、葉柵空氣動力學國家級科技重點實驗室,陜西 西安 710072)

0 引 言

風力機依靠風輪葉片獲取風能,葉片的氣動性能好壞直接影響著風力機的總體效率。構成葉片氣動外形的翼型決定著葉片的性能,因此準確計算風力機翼型氣動性能是進行葉片性能分析和氣動設計的重要基礎。傳統的風力機通常選用NACA系列的航空翼型,這些翼型對前緣粗糙度非常敏感,前緣受到污染會導致翼型性能大幅度下降[1]。鑒于航空翼型不能很好地滿足風力機及其特殊運行環境的要求,80年代后期風電發達國家開始對風力機翼型進行研究,并成功開發出風力機葉片專用翼型,比如美國NREL的S系列、丹麥 RIS ?-A 系列、瑞典 FFA-W 系列和荷蘭DU系列[2],目前DU翼型已被直徑從29m到100m、最大功率從350KW到3.5MW的十多種不同型號的風電機組所采用[3]。

設計風力機翼型的技術基礎大多是從60年代后期發展起來的,隨著計算機硬件及軟件技術的發展,先進的CFD技術已廣泛應用于不同類型氣動外形的計算。針對低雷諾數、高升力系數狀態下的風力機運行條件,采用考慮粘性的NS控制方程分析葉片翼型的流場是非常必要的。由于風力機的工作環境和運行工況很復雜,來流方向不確定,風速大小多變及葉片表面存在的粗糙度,都會導致葉片上的流動更容易發生從層流向湍流的轉捩。深刻認識翼型表面邊界層內的流態,對于正確預估翼型升力和阻力、控制并減小流動分離以及翼型的優化設計有著重要意義。然而,由于對轉捩機理認識的不足,多數對風力機翼型的繞流數值模擬都采用全湍流計算,難以準確地預估翼型表面邊界層中轉捩、分離等復雜現象。

為了更可靠地預測風力機翼型的氣動性能,本文在求解RANS方程的過程中耦合Mark Drela[4]提出的簡化eN轉捩判斷方法進行轉捩點位置的自動判斷,在翼型繞流數值計算中考慮流動轉捩因素提高現有流場求解器的計算精度。

1 計算方法

在NS方程中耦合eN方法進行轉捩點判斷是一個迭代的過程,每一次轉捩點位置的迭代需要調用三組程序。

首先是繞風力機翼型流動的計算程序,其功能是給層流邊界層方程的求解提供物面壓力分布。在連續介質假設下,忽略徹體力和熱源的二維雷諾平均NS方程可表述為如下形式:

其次是二維不可壓縮層流邊界層方程求解程序,其功能是由RANS方程提供的物面壓力分布確定外邊界并求解出翼型表面層流邊界層的解。本文采用二維不可壓縮層流邊界層方程,通過Falkner-Skan坐標變換轉變原邊界層方程的形式以降低邊界層方程的求解難度,采用Keller提出的“BOX”格式離散求解邊界層方程。

再者就是Mark Drela提出的的eN數據庫轉捩判斷方法,這一方法將層流邊界層的動量厚度、位移厚度及形狀因子等參數與擾動的放大因子N聯系起來,并給出了他們之間具體的關系,因此,利用邊界層的解可以快速的得到邊界層流動發生轉捩的位置。eN數據庫的具體形式為:

動量厚度雷諾數Reθ和空間擾動放大因子N的關系為:

臨界動量厚度雷諾數Reθ可由形狀因子來確定:

那么空間擾動放大因子N可由下式得到:

本文利用上述程序利用如圖1所示流程實現在雷諾平均NS方程的迭代計算中耦合轉捩點的自動判斷,具體過程參見文獻[7]。

2 算例及結果分析

根據上述方法本文對馬赫數0.1,雷諾數1.0×106,攻角范圍為-15.8°～ 28.1°的 DU91-W2-250翼型繞流分別進行了全湍流和耦合轉捩判斷的數值模擬計算,翼型幾何形狀如圖2所示。計算網格為貼體的512×128 C型結構化網格,外邊界距離翼型表面為15倍弦長,其中翼型上下表面共分布384個點,對翼型前后緣和表面處的網格進行了加密處理,圖3所示為網格局部放大圖。

圖1 轉捩點位置的自動判斷過程Fig.1 Aotomatic transition prediction process

圖2 翼型幾何形狀Fig.2 Airfoil profile

圖3 計算網格放大圖Fig.3 Computation grid with a local zoom

圖4給出了采用全湍流和耦合轉捩判斷兩種方法計算得到的DU91-W2-250翼型升力、阻力特性曲線并與實驗值的比較結果,圖中實驗數據記為“Experimental”,全湍流計算結果記為“Cal-Fully Turbulent”,加轉捩判斷的計算結果記為“Cal-Transitional”。從圖中可以明顯看出,在風力機翼型的使用工況攻角變化范圍(-10°～10°)內,加轉捩判斷的計算結果與實驗數據基本吻合,而全湍流計算的升力系數在-5°～10°攻角范圍內比實驗值要低0.1左右,阻力系數在攻角-10°～10°范圍內也略高于實驗值,這說明在考慮流動轉捩因素后的雷諾平均NS方程計算對翼型的升阻力特性的預測精度都有了一定提高。同樣可以看到當攻角小于-10°或大于10°后,計算結果與實驗值變化趨勢基本一致,但計算的升力系數明顯偏高而阻力系數偏小,這是由于數值計算所采用的湍流模型與真實情況還存在差異,高估了升力,低估了阻力,難以準確模擬大分離流動。

圖4 翼型升、阻力系數隨攻角變化曲線Fig.4 Curves of lift and drag coefficient versus angle of attack

圖 5分別給出了翼型攻角為7.686°、9.742°時采用全湍流和耦合轉捩判斷計算的表面壓力系數分布與實驗數據的對比。從圖可以看出,計算結果與實驗測量結果基本一致,尤其是在攻角為7.686°、9.742°時,加轉捩判斷的壓力分布與實驗值符合很好,而全湍流計算的上翼面前半段壓力分布與前兩者相差較大,表明考慮轉捩影響后的雷諾平均NS方程求解能更準確地模擬翼型繞流流場。

圖5 翼型表面壓力分布計算與實驗結果比較Fig.5 Comparison of the pressure distribution between the calculated and the experimental data

3 小 結

本文以風力機翼型DU91-W2-250為研究對象,通過求解雷諾平均NS方程對其繞流進行數值模擬,計算采用SA湍流模型,并耦合eN數據庫方法進行流動轉捩的自動判斷。采用該方法計算了雷諾數1.0×106時,攻角在-15.8°～28.1°范圍內該翼型的升力和阻力特性曲線及表面壓力分布,并與全湍流計算結果和實驗數據比較,對比結果顯示:

(1)在正常工況攻角范圍內,耦合轉捩判斷的計算結果與實驗結果相比符合很好,與全湍流假設得到的計算結果相比與實驗值更吻合,說明要準確模擬風力機翼型使用工況下的氣動特性,必須考慮轉捩帶來的影響;

(2)隨著攻角不斷增大,轉捩點位置提前到翼型前緣,導致翼型的整個上表面幾乎全是湍流,此時再考慮流動轉捩因素進行轉捩判斷的計算與全湍流計算的效果相同。

[1]葉枝全,黃繼雄,陳嚴.適用于風力機的新翼型氣動性能的實驗研究[J].太陽能學報,2003,24(4):548-554.(YE Z Q,HUANG J X,CHEN Y.Experimental study on aerodynamic performances of new airfoils for windturbine[J].Acta Energiae Solaris Sinicas,2003,24(4):548-554.)

[2]BERTAGNOLIO F,SORENSEN N,JOHANSEN J,FUGLSANG P.Wind turbine airfoil catalogue[R].Ris? National Laboratory,Roskilde,Denmark,August 2001.

[3]TIMMER W A,ROOIJ R P J O M.Summary of the delft university wind turbine dedicated airfoils[R].AIAA-2003-0352.

[4]DRELA M,GILES M B.Viscous-inviscid analysis of transonic and low reynold number airfoils[J].AIAA Journal,1986,25(10):1347-1355.

[5]韓忠華.旋翼繞流的高效數值計算方法及主動流動控制研究[D].西安:西北工業大學,2007.(HAN Z H.Efficient method for simulation of viscous flows past helicopter rotors and active flow control[D].Xi'an:Northwestern Polytechnical University,2007.)

[6]SPALART P R,ALLMARAS S R.A one-equation turbulence model for aerodynamic flows[R].AIAA-92-0439.

[7]張坤,宋文萍.基于轉捩自動判斷的翼型阻力精確計算方法研究[A].第六屆西北地區計算物理學術會議論文集[C].2008:274-279.(ZHANG K,SONG W P.Accurate airfoils drag calculations based on the automatic transition prediction[A].In:6th computational physics conference papers of northwest in China[C].2008,pp.274-279.)