覆冰對風(fēng)力機(jī)專用翼型氣動(dòng)性能影響的數(shù)值研究

張聘亭, 楊 濤, 盧緒祥

(1.華中科技大學(xué)能源與動(dòng)力工程學(xué)院430074;2.長沙理工大學(xué)能源與動(dòng)力工程學(xué)院410114)

風(fēng)能是一種清潔可再生能源,面對越來越嚴(yán)重的能源危機(jī)和環(huán)境問題,利用綠色風(fēng)能是普遍的趨勢和潮流.截至2010年年底,中國風(fēng)力發(fā)電的裝機(jī)容量已經(jīng)達(dá)到了44 730 MW,占全球裝機(jī)容量的22.4%[1].我國風(fēng)能資源主要集中在北部(東北、華北、西北)地區(qū)、東南沿海及其島嶼以及內(nèi)陸個(gè)別地區(qū),不同地域的環(huán)境差別很大,大部分地區(qū)氣候變化急劇,溫差幅度較大,都有不同程度的冰凍.處于臨界狀態(tài)的雨、雪、霧或露遇到低溫的設(shè)備和金屬結(jié)構(gòu)表面時(shí)會在設(shè)備以及結(jié)構(gòu)的表面結(jié)冰,因此運(yùn)行在寒冷氣候環(huán)境中的風(fēng)力機(jī)經(jīng)常面臨嚴(yán)重的結(jié)冰問題[2-3].風(fēng)力機(jī)葉片的覆冰問題在20世紀(jì)90年代就引起了國外學(xué)者的關(guān)注:Bose[4-5]研究了小型水平軸風(fēng)力機(jī)覆冰的過程和形態(tài),并觀測了覆冰過程中風(fēng)力機(jī)輸出功率的變化,對壓力面和吸力面的覆冰狀態(tài)做了分析;隨著CFD技術(shù)的發(fā)展,國內(nèi)外的學(xué)者開始利用CFD技術(shù)對風(fēng)力機(jī)葉片的覆冰問題進(jìn)行研究,Fortin等[3]和Imamura等[6]利用CFD技術(shù)對風(fēng)力機(jī)旋轉(zhuǎn)過程中的覆冰過程進(jìn)行了模擬,Fortin等[3]還對風(fēng)力機(jī)葉片融冰過程進(jìn)行了模擬,并對覆冰后風(fēng)力機(jī)葉片的性能進(jìn)行了分析;李聲茂等[7]對寒冷氣候條件下直線翼垂直軸風(fēng)力機(jī)葉片表面結(jié)冰進(jìn)行了數(shù)值模擬.

風(fēng)力機(jī)葉片是風(fēng)力機(jī)獲取風(fēng)能的關(guān)鍵部位,其氣動(dòng)性能直接關(guān)系到風(fēng)力機(jī)的運(yùn)行效率,甚至關(guān)系到風(fēng)力機(jī)運(yùn)行的安全性.風(fēng)力機(jī)翼型的升阻力曲線直接關(guān)系到風(fēng)力機(jī)的性能,其斜率也影響著風(fēng)力機(jī)運(yùn)行的穩(wěn)定性[8],覆冰所導(dǎo)致的翼型失速被推遲也會致使風(fēng)力機(jī)過載[2].因此,研究風(fēng)力機(jī)翼型在覆冰情況下的氣動(dòng)性能改變,對風(fēng)力機(jī)覆冰狀態(tài)下運(yùn)行決策起著重要的作用.

筆者采用通用CFD軟件Fluent 6.3,利用S-A湍流模型,對風(fēng)力機(jī)專用翼型S809不同覆冰狀態(tài)下的靜態(tài)失速特性進(jìn)行了研究,得到了靜態(tài)失速特性下的升阻力曲線,并利用 Fluent的自定義函數(shù)UDF以及動(dòng)網(wǎng)格對覆冰翼型的動(dòng)態(tài)失速特性進(jìn)行了研究.通過和非覆冰翼型的比較,得到了覆冰對翼型性能和動(dòng)態(tài)失速特性的影響.

1 數(shù)值模型

1.1. S-A湍流模型

Spalart-Allmaras湍流模型由Spalart和 Allmaras于1992年提出[9].該模型是一種相對簡單的通過求解工作變量ν～的輸運(yùn)方程得到湍流運(yùn)動(dòng)黏度νt的單方程湍流模型,其具體形式如下:

S-A模型在計(jì)算受制于壓力梯度的邊界層流動(dòng)中取得了很好的結(jié)果,因此,S-A模型被普遍應(yīng)用于透平機(jī)械的數(shù)值計(jì)算中.文獻(xiàn)[10]～文獻(xiàn)[12]分別利用S-A模型對翼型動(dòng)態(tài)失速特性進(jìn)行了仿真,得到了比較正確的結(jié)果.

1.2 數(shù)值方法

S809翼型是美國國家可再生能源實(shí)驗(yàn)室(NREL)設(shè)計(jì)的水平軸風(fēng)力機(jī)專用翼型之一.根據(jù)Bose[3]對風(fēng)力機(jī)葉片覆冰形態(tài)的研究,筆者建立了S809翼型以及S809不同覆冰形態(tài)翼型的數(shù)值模型.本文的數(shù)值模型在Gambit2.0中建立,整體計(jì)算域如圖1所示,總長度為45c,高度為30c,c為翼型的弦長.圖2為S809翼型及其周圍網(wǎng)格的局部放大圖,翼型周圍采用結(jié)構(gòu)化C型網(wǎng)格.在正式計(jì)算之前,對計(jì)算網(wǎng)格進(jìn)行了無關(guān)性分析,翼型表面分布了300點(diǎn),邊界層第一層高度為翼型弦長的 1×10-4,總網(wǎng)格數(shù)為18萬.根據(jù)文獻(xiàn)[3]對不同溫度及環(huán)境條件下覆冰形態(tài)數(shù)值仿真的結(jié)果建立了翼型可能的覆冰模型,并給出了周圍的網(wǎng)格,見圖3.S809翼型未覆洋情況定義為R0.

圖1 整體計(jì)算域Fig.1 The region of integ ral computation

圖2 S809翼型及其周圍的網(wǎng)格局部放大圖Fig.2 S809 airfoil and local enlargement of the surrounding g rid

圖3 翼型不同的覆冰形態(tài)及其周圍的網(wǎng)格局部放大圖Fig.3 Different icing patterns on airfoil and local enlargement of the surrounding grid

數(shù)值計(jì)算中對連續(xù)方程、二維動(dòng)量方程、能量方程、S-A單方程湍流模型四個(gè)方程進(jìn)行耦合求解,方程采用隱式格式,當(dāng)流場中所有計(jì)算點(diǎn)上的最大殘差小于1×10-6時(shí),認(rèn)為計(jì)算收斂.

2 靜態(tài)翼型數(shù)據(jù)的求解和討論

靜態(tài)翼型數(shù)值計(jì)算的邊界條件:來流馬赫數(shù)為0.3,雷諾數(shù)為1×106.計(jì)算不同攻角下S809翼型的升力系數(shù),并分析不同覆冰情況對升力系數(shù)的影響.計(jì)算得到S809翼型的升力系數(shù)曲線如圖4所示.

圖4 S809翼型不同覆冰情況下升力系數(shù)Fig.4 T he lift coefficient under different icing conditions of S809 airfoil

在翼型覆冰朝吸力面方向生長的情況下,如圖3(a)和圖3(b)的R1和R2覆冰狀態(tài)下,覆冰降低了翼型的彎曲,造成了翼型升力系數(shù)的下降.當(dāng)攻角較小時(shí)(＜6°),翼型的升力系數(shù)改變并不明顯,但在較大的攻角范圍內(nèi),無論是在線性區(qū)域還是失速區(qū)域,R1和R2的覆冰情況都帶來了升力系數(shù)的較大下降.當(dāng)攻角為0°時(shí),在覆冰形態(tài) R1的情況下,覆冰層后就已經(jīng)出現(xiàn)了分離渦;且隨著攻角的增加,分離渦不斷地向尾緣方向生長(見圖5),造成了升力系數(shù)的較大下降.對于覆冰形態(tài)R2,覆冰層后分離渦的生長慢于R1,但同樣也造成了升力系數(shù)的下降,見圖6.

圖7為攻角同為8°時(shí),不同的覆冰形態(tài)對葉片周圍分離渦的影響.從圖7可以看出:和R0覆冰形態(tài)相比,R1和R2的覆冰形態(tài)不但造成了前緣的分離渦,并且也加速了尾部分離渦的生成,造成了翼型失速的提前,引起升力系數(shù)的下降;對于R3覆冰形態(tài),覆冰沿弦向生長,當(dāng)攻角較小時(shí),翼型的升力系數(shù)改變并不明顯,但隨著攻角的不斷增加,R3覆冰形態(tài)使失速提前了,造成了升力系數(shù)的下降;對于R4覆冰形態(tài),覆冰主要朝壓力面方向生長,覆冰增加了翼型彎度,在一定的程度下反而提高了升力系數(shù),并且R4覆冰形態(tài)推遲了翼型尾緣分離渦的形成,一定程度上推遲了失速,增加了升力系數(shù),可能造成風(fēng)力機(jī)運(yùn)行的過載.

圖5 R1覆冰形態(tài)下覆冰層后的局部流場Fig.5 Local flow field following the icing layer of R1 pattern

圖6 R2覆冰形態(tài)下覆冰層后的局部流場Fig.6 Local flow field following the icing layer of R2 pattern

圖7 攻角為8°時(shí),S809翼型不同覆冰情況下翼型周圍的流場Fig.7 The surrounding flow field under different icing conditions of S809 airfoil at an attack angle of 8

3 動(dòng)態(tài)翼型數(shù)據(jù)的求解和討論

動(dòng)態(tài)翼型的數(shù)值計(jì)算邊界條件同靜態(tài)翼型數(shù)據(jù):來流馬赫數(shù)為0.3,雷諾數(shù)為1×106.根據(jù)文獻(xiàn)[12],翼型繞 x/c=0.25(x表示該點(diǎn)到翼型前緣的距離)處作正弦周期振蕩,正弦變化規(guī)律為:

式中 :α0=20°;Δα=11°;取衰減頻率為 k=0.4,則振蕩頻率ω=83.132 5;時(shí)間步長取0.001 26 s.

經(jīng)過4個(gè)周期的計(jì)算后得到了基本一致的滯回曲線,則認(rèn)為流場趨于穩(wěn)定.分別計(jì)算S809翼型及其覆冰情況下的動(dòng)態(tài)升力系數(shù)曲線,如圖8所示.

圖8 S809翼型不同覆冰狀態(tài)下的動(dòng)態(tài)升力系數(shù)曲線Fig.8 The dy namic lift coefficient curve under different icing conditions of S809 airfoil

從翼型的升力系數(shù)曲線來看,其動(dòng)態(tài)失速特性和靜態(tài)失速特性有著明顯的區(qū)別.在靜態(tài)流場中,翼型的升力系數(shù)和攻角是一一對應(yīng)的.然而在動(dòng)態(tài)流場中,當(dāng)翼型突然從低于失速攻角增加到高于失速攻角時(shí),則可以看到翼型升力系數(shù)相對于靜態(tài)升力系數(shù)的過調(diào),在翼型的振蕩過程中,翼型的升力系數(shù)曲線形成了一條滯回曲線,這種動(dòng)態(tài)失速特性也是失速顫振氣動(dòng)負(fù)阻尼的來源.

結(jié)合圖8和圖9分析覆冰對動(dòng)態(tài)失速特性的影響.對于S809翼型,和未覆冰狀態(tài)(見圖9(a))相比,R1和R2覆冰狀態(tài)提前了翼型失速過程的發(fā)生,而R4覆冰狀態(tài)(圖9(d))則推遲了翼型失速過程的發(fā)生.這和靜態(tài)翼型數(shù)據(jù)的分析結(jié)論是一致的,也可以看出覆冰對翼型失速的影響.對于R4覆冰狀態(tài),其升力系數(shù)曲線的斜率要大于 R0無覆冰狀態(tài)(見圖8),這從氣動(dòng)彈性穩(wěn)定的角度考慮是有益的.

圖9 S809翼型不同覆冰狀態(tài)下翼型周圍流場圖Fig.9 The surrounding flow field under different icing conditions of S809 airfoil

4 結(jié) 論

(1)向吸力面生長的覆冰形態(tài)會造成覆冰層后的分離渦,并且隨著攻角的增加分離渦會向后緣生長,造成升力系數(shù)的較大下降,同時(shí)也會造成阻力系數(shù)的增加;

(2)翼型的前緣覆冰會對翼型尾緣分離渦的生成產(chǎn)生影響,對于向吸力面生長的覆冰,會促進(jìn)尾緣分離渦的生成,而向壓力面生長的覆冰則能抑制尾緣分離渦的生長,從而造成對翼型升力系數(shù)的影響;沿弦向生長的覆冰情況,則對尾緣分離渦的生成影響較小;

(3)在動(dòng)態(tài)失速情況下,翼型周圍流場比較復(fù)雜,覆冰的形態(tài)對翼型升力系統(tǒng)造成影響規(guī)律也比較復(fù)雜,但其分離渦生長的基本規(guī)律同靜態(tài)流場時(shí);

(4)從翼型覆冰情況下的靜態(tài)流場特性分析,覆冰破壞了翼型的流線,從而直接影響了翼型的氣動(dòng)性能.當(dāng)翼型升力系數(shù)降低時(shí),風(fēng)力機(jī)的出力將受到影響,經(jīng)濟(jì)性降低;當(dāng)翼型失速被推遲,又容易造成風(fēng)力機(jī)的過載,影響安全性;

(5)從翼型覆冰情況下的動(dòng)態(tài)流場特性分析,覆冰改變了翼型動(dòng)態(tài)升力系數(shù)曲線的斜率;而升力系數(shù)曲線的斜率將影響風(fēng)力機(jī)的氣動(dòng)彈性穩(wěn)定性;

(6)由于風(fēng)力機(jī)覆冰情況復(fù)雜,不同的覆冰厚度、不同的覆冰形態(tài)等都會對覆冰翼型的氣動(dòng)性能造成影響.本文的工作僅是對該問題從數(shù)值模擬的角度開展前期的研究.為了了解風(fēng)力機(jī)在覆冰情況下的運(yùn)行狀況,提高風(fēng)力機(jī)在覆冰情況下運(yùn)行的安全性和穩(wěn)定性,還需要長期的探索.

[1]李俊峰,蔡豐波,唐文倩,等.風(fēng)光無限:中國風(fēng)電發(fā)展報(bào)告 2011[M].北京:中國環(huán)境科學(xué)出版社,2011.

[2]DALILI N,EDRISY A,CARRIVEAU R.A review of surface engineering issues critical to wind turbine performance[J].Renewable and Sustainable Energy Reviews,2009,13(2):428-438.

[3]FORTIN G,PERRON J.Wind turbine icing and deicing[C]//47th AIAA Aerospace Sciences Meeting Including The NewHorizons Forumand Aerospace Exposition.Orlando,USA:AIAA,2009.

[4]BOSE N.Icing on a small horizontal-axis wind turbine-part 1:glaze ice profiles[J].Journal of Wind Engineering and Industrial Aerodynamics,1992,45(1):75-85.

[5]BOSE N.Icing on a small horizontal-axis wind turbine-part 2:three dimensional ice and wet snow formations[J].Journal of Wind Engineering and Industrial Aerodynamics,1992,45(1):87-96.

[6]IMAM URA Y,TODA K,YAMAMOTO M.Numerical simulation of ice accretion on wind turbine blade[C]//Abstracts of the Papers Presented at the Minisymposia Sessions of WCCM VI in Conjunction with APCOM.Beijing:[s.n.],2004.

[7]李聲茂,李巖.垂直軸風(fēng)力機(jī)葉片翼型結(jié)冰的數(shù)值計(jì)算[J].動(dòng)力工程學(xué)報(bào),2011,31(3):214-219.LI Shengmao,LI Yan.Numerical simulation on icing of a blade aerofoil for vertical axis wind turbines[J].Journal of Chinese Society ofPower Engineering,2011,31(3):214-219.

[8]HANSEN M.Aerodynamics of wind turbines[M].[S.1.]:Earthscan,2008.

[9]SPALA RT P R,ALLM ARAS S R.A one-equation turbulence model for aerodynamic flows[J].La Recherche Aerospatiale,1994,1(1):5-21.

[10]李杰.風(fēng)力機(jī)葉片翼型氣動(dòng)特性數(shù)值研究[D].北京:華北電力大學(xué)能源動(dòng)力與機(jī)械工程學(xué)院,2009.

[11]趙旭,肖俊,席德科.水平軸風(fēng)力機(jī)翼型設(shè)計(jì)與動(dòng)態(tài)失速數(shù)值模擬[J].太陽能學(xué)報(bào),2009,30(3):348-354.ZHAO Xu,XIAO Jun,XI Deke.The design of airfoils and the simulation of dynamic stall of horizontal axis wind turbines[J].Acta Energiae Solaris Sinica,2009,30(3):348-354.

[12]陳旭,郝輝,田杰,等.水平軸風(fēng)力機(jī)翼型動(dòng)態(tài)失速特性的數(shù)值研究[J].太陽能學(xué)報(bào),2003,24(6):735-740.CHEN Xu,HAO Hui,TIAN Jie,et al.Investigation on airfoil dynamic stall of horizontal axis wind turbine[J].Acta Energiae Solaris Sinica,2003,24(6):735-740.