不同風速下風力機動態特性研究

杜 標，李雪斌，王 龍，李 亮，孫倫業

（1.安徽理工大學 機械工程學院，淮南 232001；2.安徽理工大學 力學與光電物理學院，淮南 232001）

0 引言

風電作為一種環保可再生資源受到世界各國的重視，而在風力機組工作中，其成本的75%～90%主要來源于制造和維護風力機。風力機需要維護主要是因為在使用中受到復雜大氣環境的影響，其中動態失速問題尤為嚴重。動態失速指一個振蕩周期內的壓力面在超過其臨界迎角時繞流流場發生失速和非定常分離的現象[1]，與靜態失速相比較，葉片所受動態載荷更大，其問題更貼近風力機實際工況。

因為獲取翼型在不同工況的完整實驗數據進行的風洞試驗費用較高，利用計算流體力學（CFD）數值模擬[2]研究翼型氣動特性的方法得到廣泛使用。在眾多風力機翼型中，S809翼型是美國國家可再生能源實驗室（NREL）風力機動力學實驗[3]所采用的葉片翼型，其風洞實驗數據非常充分。S.L.Yang等人[4]在90年代最先選取S809翼型進行湍流數值模擬，在初始附著流動階段，升力系數與實驗值吻合程度很高，只有在大攻角下與實驗值不太一致。陳旭等人[5]利用S-A湍流模型模擬S809翼型俯仰運動時的動態繞流流場，并將動態與靜態繞流流場進行對比，發現在動態工況下翼型的升力和阻力系數發生了顯著的變化。McCroskey等人[6]進行動態失速翼型的氣動力實驗，對比分析翼型在不同折合頻率和馬赫數下的升阻力系數。以上研究不僅有利于深入了解翼型動態失速過程中的流動現象和改進工程中廣泛采用的經驗解析模型[7]，而且便于風力機葉片的設計制造與維護。

本文計算分析了S809翼型在不同風速下動態特性，在S809翼型1/4弦長處繞其作正弦周期振蕩，振蕩的幅度為10°，折合頻率為0.026，選擇分別在來流速度為15m/s、20m/s、25m/s、30m/s、35m/s和40m/s的工況下，對翼型進行二維數值模擬計算，并且采用算例與俄亥俄州立大學（OSU）風洞試驗的數據[8]對比驗證，研究風速對動態失速的影響以及不同時刻的流場分布。

1 控制方程及湍流模型

1.1 控制方程

本文對于動態失速中做俯仰運動的S809翼型利用CFD方法模擬計算，控制方程采用在笛卡爾坐標系下，二維連續性方程和不可壓縮流體的N-S方程[9]：

連續性方程：

N-S方程：

式中，ρ為密度，kg/m3；ui、uj（i ，j = x，y，z）為時均速度分量，m/s；μ為流體動力黏性系數，u為脈動平均速度，m/s；kg/(m·s)；P為流體靜壓，kg·m/s；為雷諾應力項，kg/(m·s2)；fi為體積力，kg/（m·s）2；Fi為附加源項。

1.2 S-A湍流模型

渦粘模式在工程湍流問題中得到廣泛應用，其雷諾應力為：

對于翼型的動態繞流流場，本文選取了Spalart-Allmaras(S-A)湍流模型[10]。其表達式如下：

2 參數設置及網格劃分

2.1 數值模擬參數

本文研究的對象為繞1/4弦長處作正弦運動的S809翼型。參數設置如下：折合頻率為k=0.026，初始攻角α0=8°，震蕩幅角α1=10°，來流速度分別為15m/s、20m/s、25m/s、30m/s、35m/s和40m/s共六種工況。

翼型振蕩規律的函數如下：

式中，α(t)為瞬時攻角，α0為平均攻角，α1為振蕩幅度，f為振蕩頻率，t為時間。

在振蕩翼型中，折合頻率是描述問題的非定常程度的一個重要參數。定義折合頻率的公式如下：

折合頻率數值大小不僅反映了振蕩運動對翼型主流運動影響的程度，而且體現其主流運動和繞轉軸方向振蕩運動的互相作用。

2.2 葉片網格劃分

本文計算S809翼型的振蕩運動選用了動網格方法[12]，圖1分別為翼型計算域網格和放大后翼型附近網格。整個計算域均采用三角形網格，共計約26400個網格單元，計算域分為內部動區域和外部靜止區域。從外部靜止區域到內部交界處網格相近以此減小插值誤差。翼型處的網格密度較大，可在圖1(b)放大圖呈現。翼型計算域在進口邊界設置來流速度，出口邊界設置為靜壓力，在翼型表面邊界滿足無滑移壁面條件。

圖1 S809翼型網格

模擬計算設置時間步長為0.001秒，且需要檢查時間步長的無關性，為了得到周期性穩定重復的流場，每個工況經過6個以上連續周期的計算，并提取最穩定的數據進行后處理以及與OSU實驗數據對比。

3 計算結果與分析

3.1 算例驗證

圖2是S809翼型在風速32m/s條件下靜態失速和動態失速升力系數遲滯圖算例驗證，相關計算參數如下：雷諾數Re=106，折合頻率k=0.050，振蕩頻率f=1.2Hz，初始攻角α0=8°，振蕩幅角α1=5.5°。如圖2(a)所示，靜態失速中計算值在小攻角下與實驗值從圖中可以看出升力系數計算值與OSU實驗值在8°攻角前十分接近，在攻角超過8°以后計算值略大于實驗值，實驗值趨于平穩，說明此時出現流動分離。

由圖2(b)可知，利用S-A湍流模型模擬計算的氣動力遲滯閉環曲線與OSU實驗值具有良好的一致性，兩個遲滯閉環曲線均為“8”字形，CFD計算值在剛開始的上仰過程中升力系數值與OSU實驗值相持平，當攻角達到10°以后，計算值略大于實驗值并達到峰值，因為在翼型表面出現漩渦，并且漩渦會附著在翼型表面上。在下俯階段，大攻角下的計算值大于實驗值，此時存在流動分離現象，湍流模型對于該現象模擬的精確性直接影響計算結果。此算例驗證說明本文所采用的CFD方法模擬風力機翼型動態失速是可行的。

圖2 翼型算例驗證

3.2 氣動力系數遲滯閉環

圖3(a)～(f)分別是來流速度為15m/s、20m/s、25m/s、30m/s、35m/s和40m/s共六個工況的升力系數遲滯閉環圖，模擬計算均選擇振蕩幅角為10°。由圖3六種工況的升力系數遲滯閉環規律可得，當攻角到達13°～15°范圍內會出現最大升力系數，而且當風速不斷增大時，達到升力系數的峰值攻角就越大；升力系數的峰值大小在1.3～1.5之間，同樣與風速呈線性增長關系；由圖3(f)可知，當風速達到40m/s時，升力系數峰值有最值1.5；升力系數遲滯回歸線所包圍的面積在風速為35m/s之內隨著風速的變大而增大，隨后風速增大而面積減小。

在此次工況計算模擬下，由于所取的折合頻率較低，流動有部分都處于附著狀態，當攻角逐漸增大會出現流動分離現象。當來流速度逐漸增大時，翼型上仰到大攻角以及后面的下俯過程中會發生振蕩波動，原因在于翼型表面上有分離渦和部分逆流，尾緣處分離渦會逐漸脫落。當翼型在高風速大攻角下運動時，翼型上的振蕩載荷逐漸增大，且超出其設計值導致結構受損，所以研究葉片動態失速升力系數遲滯規律對風力機性能維護和受損預防有重要意義。

圖3 升力系數遲滯閉環圖

3.3 流場漩渦的發展

為了更好的研究風力機翼型動態失速的流場特性以及流場漩渦的發展過程，本文選取了風速v=30m/s工況下同一振蕩周期中不同時刻的翼型流線圖。如圖4所示，分別是翼型在同一周期內上仰角度為10.11°、14.53°、17.30°、17.91°以及下俯角度為16.11°、14.85°、11.77°、7.44°的流場圖。在上仰角達到14.53°左右，開始出現尾緣渦，在此之前，流動處于完全附著狀態，升力系數隨攻角的增大而增大。隨著時間的流動，升力系數達到峰值以后，風力機翼型尾緣產生了反方向的漩渦，當攻角進一步增大后，漩渦開始向前緣延伸，升力開始逐漸減小，圖4(d)表明翼型即將達到最大值18°，此時漩渦開始變小且有脫落的趨勢。隨后翼型開始下俯階段，尾緣渦逐漸變小直到攻角為14°到后基本消散，隨著攻角的逐漸減小，流動分離點向翼型的后端移動直到完全貼附。

4 結論

本文使用CFD方法，并利用S-A湍流模型，以S809翼型1/4弦長點為旋轉中心，對作周期正弦振蕩的翼型流場進行模擬計算。選取來流速度15m/s至40m/s范圍中不同風速的6個典型工況，分別得到不同工況下的升力系數遲滯閉環圖和同一振蕩周期中的流場分布圖。并得出以下結論：

圖4 不同時刻流場圖

1）先用算例與OSU風洞實驗進行對比，研究升力系數遲滯規律，兩個遲滯閉環曲線均為“8”字形，數值大小十分接近，可得CFD計算值與OSU實驗值具大小趨勢比較一致，證明了CFD方法計算風力機動態失速的可行性。

【】【】

2）通過不同風速下的升力系數遲滯圖可得，隨著風速的逐漸增大，升力系數的峰值呈線性增長關系；達到升力系數峰值的攻角越大，升力系數的峰值大小也逐漸增加并且范圍在1.3～1.5之間；升力系數遲滯閉環包圍的面積先隨著風速的增加而變大，當風速大于35m/s后，遲滯閉環面積開始減小。

3）在風速為30m/s工況的同一振蕩周期不同時刻的流場中，在上仰過程中，升力增大的同時出現流動分離現象，導致尾緣端出現漩渦，達到峰值會出現逆向小渦，并且當攻角不斷增大時，漩渦開始逐漸脫落；在下俯階段，隨著攻角的減小，尾緣渦逐漸變小直至消散。

[1]周姣,楊科,張磊,等.基于大渦模擬的風力機翼型失速機制探討[J].工程熱物理學報,2013(11):2048-2051.

[2]文曉慶,柳陽威,方樂,陸利蓬.提高k-ωSST模型對翼型失速特性的模擬能力[J].北京航空航天大學學報,2013,39(8):1127-1132.

[3]HAND M M,SIMMS D A,FINGERSH L J,et al. Unsteady aerodynamics experiment phase VI: wind tunnel test configurations and available data campaigns[R].Colorado: National Renewable Energy Laboratory,2001.

[4]S.L.Yang.Incompressible Navier- Stokes Computation of the NREL Airfoils Using a Symmetric Total Variational Diminishing Scheme[J].Journal of Solar Energy Engineering,1994,116(4):174-182.

[5]陳旭,郝輝,田杰,杜朝輝.水平軸風力機翼型動態失速特性的數值研究[J].太陽能學報,2003,24(6):735-740.

[6]G.R.Srinivasan,W.J.Mccroskey.Navier-Stokes calculations of hovering rotor flowfields[J].Journal of Aircraft,1988,25(10):865-874.

[7]俞國華.水平軸風力機葉片失速問題研究[D].上海交通大學,2013.

[8]Ramsay R R，Hoffman M J，Gregorek G M.Effects of grit roughness and pitch oscillations on the S809 airfoil[R].Golden,Colorado:National Renewable Energy Laboratory，1999,1-165.

[9]周光垌.流體力學,上冊[M].高等教育出版社,1992.

[10]Zhang Y, Bai J Q, Xu J L,et al.A one-equation turbulence model for recirculating flows[J].Science China Physics Mechanics &Astronomy,2016,59(6):664701.

[11]Dennis DJ.Coherent structures in wall-bounded turbulence[J].Anais Da Academia Brasileira De Ciências,2015,87(ahead):513-537.

[12]Long Wang, Lun-ye Sun, Guang Wu, et al.Research on algorithm of blade vibration for general wind turbine[A],Proc.SPIE9903,Seventh International Symposium on Precision Mechanical Measurements[C],990325.