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

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

0 引言

風(fēng)電作為一種環(huán)保可再生資源受到世界各國的重視，而在風(fēng)力機(jī)組工作中，其成本的75%～90%主要來源于制造和維護(hù)風(fēng)力機(jī)。風(fēng)力機(jī)需要維護(hù)主要是因?yàn)樵谑褂弥惺艿綇?fù)雜大氣環(huán)境的影響，其中動(dòng)態(tài)失速問題尤為嚴(yán)重。動(dòng)態(tài)失速指一個(gè)振蕩周期內(nèi)的壓力面在超過其臨界迎角時(shí)繞流流場發(fā)生失速和非定常分離的現(xiàn)象[1]，與靜態(tài)失速相比較，葉片所受動(dòng)態(tài)載荷更大，其問題更貼近風(fēng)力機(jī)實(shí)際工況。

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

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

1 控制方程及湍流模型

1.1 控制方程

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

連續(xù)性方程：

N-S方程：

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

1.2 S-A湍流模型

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

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

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

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

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

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

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

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

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

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

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

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

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

3 計(jì)算結(jié)果與分析

3.1 算例驗(yàn)證

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

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

圖2 翼型算例驗(yàn)證

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

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

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

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

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

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

4 結(jié)論

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

圖4 不同時(shí)刻流場圖

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

【】【】

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

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

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

[2]文曉慶,柳陽威,方樂,陸利蓬.提高k-ωSST模型對(duì)翼型失速特性的模擬能力[J].北京航空航天大學(xué)學(xué)報(bào),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]陳旭,郝輝,田杰,杜朝輝.水平軸風(fēng)力機(jī)翼型動(dòng)態(tài)失速特性的數(shù)值研究[J].太陽能學(xué)報(bào),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]俞國華.水平軸風(fēng)力機(jī)葉片失速問題研究[D].上海交通大學(xué),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]周光垌.流體力學(xué),上冊(cè)[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.