杜 標，李雪斌，王 龍，李 亮，孫倫業(yè)

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

0 引言

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

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

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

1 控制方程及湍流模型

1.1 控制方程

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

連續(xù)性方程：

N-S方程：

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

1.2 S-A湍流模型

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

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

2 參數(shù)設置及網(wǎng)格劃分

2.1 數(shù)值模擬參數(shù)

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

翼型振蕩規(guī)律的函數(shù)如下：

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

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

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

2.2 葉片網(wǎng)格劃分

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

圖1 S809翼型網(wǎng)格

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

3 計算結果與分析

3.1 算例驗證

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

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

圖2 翼型算例驗證

3.2 氣動力系數(shù)遲滯閉環(huán)

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

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

圖3 升力系數(shù)遲滯閉環(huán)圖

3.3 流場漩渦的發(fā)展

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

4 結論

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

圖4 不同時刻流場圖

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

【】【】

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

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

[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]陳旭,郝輝,田杰,杜朝輝.水平軸風力機翼型動態(tài)失速特性的數(shù)值研究[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.