盛振國，李陳峰，任慧龍，劉小龍

(1.哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001;2.中國船級(jí)社 產(chǎn)品設(shè)計(jì)評(píng)估中心，北京 100006;3.中國船舶科學(xué)研究中心，江蘇 無錫 214082)

隨著能源和環(huán)境問題日益突出，儲(chǔ)量豐富、無污染、可再生的風(fēng)能逐漸受到人們的重視［1］.風(fēng)力發(fā)電機(jī)是將風(fēng)能轉(zhuǎn)換為電能的機(jī)械裝置，葉片是其重要組成部分.葉片以及葉輪的氣動(dòng)性能直接決定了風(fēng)力發(fā)電機(jī)組的效率.目前，國內(nèi)外有關(guān)大型風(fēng)力發(fā)電機(jī)組葉片氣動(dòng)性能研究有風(fēng)洞試驗(yàn)和數(shù)值模擬兩種方法［2-3］.風(fēng)洞試驗(yàn)數(shù)據(jù)可靠，但成本高、周期長，且單個(gè)獲得的流場(chǎng)信息有限.數(shù)值模擬方面，主要借鑒螺旋槳理論，建立了經(jīng)典的葉素動(dòng)量理論［1］(blade element-momentum，BEM).目前大多風(fēng)力機(jī)氣動(dòng)性能計(jì)算軟件是基于葉素理論開發(fā)的，葉素理論實(shí)質(zhì)上是把葉片當(dāng)作標(biāo)準(zhǔn)的二維問題來處理，難以真實(shí)地反映翼型的三維旋轉(zhuǎn)效應(yīng)和動(dòng)態(tài)失速等影響.隨著計(jì)算機(jī)技術(shù)的發(fā)展以及三維湍流技術(shù)的提高，計(jì)算流體力學(xué)方法(CFD)在研究復(fù)雜流場(chǎng)特性中起著越來越重要的作用，但在風(fēng)機(jī)氣動(dòng)性能研究中的應(yīng)用仍處于起步階段［4］.Sorensen等［2-3］使用不可壓縮RANS方法預(yù)報(bào)了風(fēng)機(jī)葉片的性能，與試驗(yàn)結(jié)果吻合較好.Madsen等［5］研究了偏航對(duì)氣動(dòng)負(fù)荷的影響.與BEM相比，CFD方法可以考慮三維旋轉(zhuǎn)效應(yīng)引起的失速延遲現(xiàn)象和動(dòng)態(tài)失速的影響，直接獲得翼型的三維氣動(dòng)特性和風(fēng)輪周圍詳細(xì)的流場(chǎng)特性，尤其對(duì)于新翼型的設(shè)計(jì)，不需要經(jīng)驗(yàn)值，即可得到功率特性.

根據(jù)大型風(fēng)力發(fā)電機(jī)組葉片的特點(diǎn)，采用CFD技術(shù)，基于RANS方法耦合SSTk-ω湍流模型，建立了風(fēng)機(jī)葉片氣動(dòng)性能分析方法，給出了二維翼型和三維葉輪氣動(dòng)性能分析的網(wǎng)格生成、邊界條件與數(shù)值求解方法等.采用本文方法對(duì)NACA64-618翼型的二維氣動(dòng)性能和2MW風(fēng)機(jī)葉片三維氣動(dòng)性能進(jìn)行了分析，與試驗(yàn)值及GHBladed軟件計(jì)算結(jié)果的比較證明了有關(guān)方法的準(zhǔn)確性，在此基礎(chǔ)上對(duì)2MW風(fēng)機(jī)翼型進(jìn)行了優(yōu)化，完善了其氣動(dòng)性能.

1 風(fēng)機(jī)氣動(dòng)性能的CFD數(shù)值模擬技術(shù)

1.1 控制方程

流體連續(xù)性方程和RANS方程如下:

式中:ui、uj為速度時(shí)均量;為速度脈動(dòng)量;ρ為密度;μ為流體粘性系數(shù);p為壓力.其中，對(duì)于二維問題，i=1，2;三維問題，i=1，2，3.

由于風(fēng)機(jī)葉輪為三維旋轉(zhuǎn)對(duì)稱結(jié)構(gòu)，因此將控制方程轉(zhuǎn)化到旋轉(zhuǎn)坐標(biāo)系下.對(duì)于靜止坐標(biāo)系下的描述速度場(chǎng)的絕對(duì)速度V與旋轉(zhuǎn)坐標(biāo)系下描述速度場(chǎng)的相對(duì)速度Vr之間的關(guān)系如下:

式中:Ω為指角速度向量(即旋轉(zhuǎn)坐標(biāo)系的角速度);r是旋轉(zhuǎn)坐標(biāo)系中的位置向量.

連續(xù)性方程:

動(dòng)量方程:

式中:τ是應(yīng)力張量，它包含粘應(yīng)力和湍流應(yīng)力2部分.對(duì)于湍流應(yīng)力項(xiàng)采用渦粘性進(jìn)行描述，即:

式中:μt為渦粘性系數(shù)，計(jì)算采用Spalart-Allmaras湍流模式.

1.2 SSTk-ω湍流模型

對(duì)于翼型氣動(dòng)性能的 CFD分析，SSTk-ω模型［6］是常用的湍流模型之一，它混合了k-ω模型和k-ε模型，使得該湍流模型同時(shí)具有了k-ε模型計(jì)算近壁面區(qū)域粘性流動(dòng)的可靠性和模型計(jì)算遠(yuǎn)場(chǎng)自由流動(dòng)的精確性.

式中:Γk和Γω為擴(kuò)散系數(shù)，μt為渦粘性系數(shù)，Gk和Gω為湍流產(chǎn)生項(xiàng)，Yk和Yω為湍流耗散項(xiàng)，α*為低雷諾數(shù)修正系數(shù)，σk和σω分別是k和ω對(duì)應(yīng)的湍流普朗特?cái)?shù)，Dω為擴(kuò)散項(xiàng).

2 二維翼型氣動(dòng)性能分析

風(fēng)力機(jī)葉片都是三維的，但是數(shù)值計(jì)算中三維計(jì)算所需的網(wǎng)格數(shù)較多、占用的計(jì)算機(jī)資源較多、計(jì)算周期較長.為了計(jì)算快捷，工程中很多計(jì)算都以二維翼型為研究對(duì)象，將空間流動(dòng)簡化成了平面流動(dòng).

2.1 二維翼型氣動(dòng)性能分析參數(shù)設(shè)置

2.1.1 坐標(biāo)系定義

坐標(biāo)原點(diǎn)角度的定義如圖1所示.所有坐標(biāo)以弦長為特征長度進(jìn)行無量綱化，因此，翼型前端為x=－0.5，翼型后端為x=0.5.

圖1 坐標(biāo)系的定義Fig.1 Definition of coordinates

2.1.2 網(wǎng)格生成與控制

在數(shù)值計(jì)算中，計(jì)算網(wǎng)格的質(zhì)量直接影響到計(jì)算結(jié)果的精度和收斂性.對(duì)于二維翼型分析，采用的網(wǎng)格種類為多塊結(jié)構(gòu)化網(wǎng)格，網(wǎng)格數(shù)量控制在6× 104左右.在靠近翼型處網(wǎng)格適當(dāng)加密，以便能夠較精確的模擬壁面附近的流動(dòng).整個(gè)網(wǎng)格分布及翼型尾部網(wǎng)格放大見圖2.在后緣處由于考慮強(qiáng)度和穩(wěn)定性方面問題時(shí)，經(jīng)常處理為隨邊厚度不為0.法;離散得到的代數(shù)方程使用Gauss-Seidel迭代求解.

圖2 網(wǎng)格劃分示意Fig.2 Computational grid domain of airfoil

2.2 二維翼型氣動(dòng)性能算例分析

NACA64-618翼型是目前較為成熟的一種翼型，本文對(duì)該翼型在攻角范圍－180°～180°的氣動(dòng)性能進(jìn)行了分析，并與試驗(yàn)結(jié)果［7］進(jìn)行了比較，如圖3所示.相關(guān)計(jì)算參數(shù)為:變槳轉(zhuǎn)矩中心位于0.25倍弦長處;雷諾數(shù)為3.0×106;計(jì)算條件取標(biāo)準(zhǔn)空氣密度1.225 g/m3，空氣運(yùn)動(dòng)粘性系數(shù)取1.789 4×10－5.

圖3 NACA64-618翼型在－180°～180°的氣動(dòng)特性曲線Fig.3 Aerodynamic performance of airfoil NACA64-618 at attack angle－180°～180°

2.1.3 求解區(qū)域和邊界條件

整個(gè)計(jì)算域?yàn)檫M(jìn)流段和去流段都為10倍弦長，而側(cè)面為6倍弦長，具體為10L×6L×10L，其中L為弦長.邊界條件為:

1)進(jìn)口邊界條件:其速度等于來流速度，即V=V∞;

2)出口邊界條件:壓力出口，假定壓力等于大氣壓;

3)外場(chǎng)邊界:外場(chǎng)邊界的設(shè)置與進(jìn)口邊界條件一致;

4)物面條件:葉片表面設(shè)定為無滑移邊界條件，即V=0.

2.1.4 數(shù)值計(jì)算方法

二維翼型氣動(dòng)性能分析通過直接求解二維粘性不可壓RANS方程，微分方程的離散使用有限體積法，其中對(duì)流項(xiàng)采用二階迎風(fēng)差分格式，擴(kuò)散項(xiàng)采用中心差分格式;壓力和速度耦合采用的SIMPLE方

由于試驗(yàn)［7］只給出了小攻角范圍內(nèi)的有關(guān)升力系數(shù)、阻力系數(shù)和俯仰力矩系數(shù).從圖3可以發(fā)現(xiàn)，當(dāng)攻角較小時(shí)，翼型的升力系數(shù)隨著攻角的增大而增大;攻角為14°時(shí)翼型的升力系數(shù)達(dá)到最大值，隨著攻角繼續(xù)增加，翼型的升力系數(shù)開始下降，而阻力系數(shù)則急劇上升，翼型的氣動(dòng)性能開始惡化.NACA64-618翼型在雷諾數(shù)為3.0×106時(shí)的失速攻角應(yīng)該在14°附近，理論結(jié)果與試驗(yàn)結(jié)果吻合.分析大攻角下理論計(jì)算結(jié)果可以發(fā)現(xiàn)，此時(shí)的氣動(dòng)性能呈非定常特性.實(shí)際應(yīng)用時(shí)，需要將有關(guān)氣動(dòng)性能在攻角度范圍－180°～180°進(jìn)行光順處理，再導(dǎo)入GHBladed或Focus等風(fēng)機(jī)分析軟件應(yīng)用.阻力系數(shù)和俯仰力矩系數(shù)基本對(duì)稱，存在的差異主要由于翼型有拱度，不是完全對(duì)稱剖面.

3 三維葉輪氣動(dòng)性能分析

二維方法可以較快捷的計(jì)算翼型的氣動(dòng)性能，但是該方法無法考慮葉片的三維效應(yīng)，尤其當(dāng)攻角達(dá)到或超過失速攻角時(shí)，葉片表面失速旋渦具有極強(qiáng)的三維性，采用三維方法可以更真實(shí)地模擬流場(chǎng)，獲得風(fēng)機(jī)的氣動(dòng)性能.

3.1 三維葉輪氣動(dòng)性能分析參數(shù)設(shè)置

3.1.1 計(jì)算區(qū)域和網(wǎng)格劃分

整體計(jì)算域圓柱體區(qū)域，見圖4.進(jìn)口在風(fēng)機(jī)前3.1D，出口在風(fēng)機(jī)后5.168D，圓柱半徑3.1D，其中為風(fēng)輪直徑.長度方向?yàn)閄方向，向下游為正，Y方向?yàn)槿~片參考線方向，半徑朝外方向?yàn)檎琙方向滿足右手法則.坐標(biāo)原點(diǎn)在風(fēng)輪的旋轉(zhuǎn)中心，風(fēng)輪截面垂直于X軸，迎著來流.

計(jì)算區(qū)域整體上被分為2個(gè)子區(qū)，包圍葉輪的圓餅型旋轉(zhuǎn)區(qū)域和其他部分的靜止區(qū)域.在靜止區(qū)域主要以結(jié)構(gòu)化網(wǎng)格進(jìn)行剖分.對(duì)于旋轉(zhuǎn)區(qū)域網(wǎng)格形式以非結(jié)構(gòu)網(wǎng)格為主，在葉片根部、梢部及導(dǎo)隨邊進(jìn)行了網(wǎng)格加密.

圖4 計(jì)算區(qū)域外圍網(wǎng)格示意Fig.4 Computational grid domain of wind turbine blade

3.1.2 邊界條件和初始條件

計(jì)算區(qū)域的進(jìn)口為速度進(jìn)口邊界條件，速度值設(shè)為來流風(fēng)速.在計(jì)算域的上邊界也同樣設(shè)為進(jìn)口速度邊界條件.在計(jì)算域的出口為壓力出口邊界條件，壓力值設(shè)為環(huán)境壓力.計(jì)算區(qū)域的兩個(gè)側(cè)面設(shè)為周期邊界條件.葉片和輪轂表面設(shè)為不可滑移物面邊界.這里不考慮地面的邊界層效應(yīng)，所以地面設(shè)為滑移物面邊界.計(jì)算初始值全域使用統(tǒng)一的來流風(fēng)速.

3.1.3 數(shù)值計(jì)算方法

控制方程采用有限體積法進(jìn)行離散，離散過程為在網(wǎng)格單元上對(duì)控制方程實(shí)施積分從而得到離散的代數(shù)方程組.這種離散方法可以保證控制方程的守恒性，具有較高的離散精度，可以方便地處理復(fù)雜幾何問題.在離散的過程中對(duì)流項(xiàng)使用二階迎風(fēng)格式，擴(kuò)散項(xiàng)采用二階中心格式，時(shí)間項(xiàng)采用一階隱式格式，對(duì)于壓強(qiáng)和速度的耦合采用SIMPLE算法.離散代數(shù)方程采用逐點(diǎn)Gauss-Seidel迭代法求解，并且采用代數(shù)多重網(wǎng)格方法加快求解的收斂速度.

3.2 三維葉輪氣動(dòng)性能算例分析

3.2.1 對(duì)象描述

2MW風(fēng)機(jī)葉輪的幾何參數(shù)、運(yùn)行參數(shù)等相關(guān)說明見表1.葉根到葉梢的翼型剖面輪廓線和三維幾何造型見圖5.

表1 2MW風(fēng)機(jī)特性描述Table1 Description of 2MW wind turbine characteristics

圖5 葉片的幾何建模Fig.5 Geometrical modeling of wind turbine blade

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

基于三維方法，建立了2MW風(fēng)機(jī)葉輪的氣動(dòng)性能分析模型，對(duì)額定風(fēng)速下的流場(chǎng)及葉片氣動(dòng)性能進(jìn)行了分析.

圖6為葉片的壓力分布云圖，可以發(fā)現(xiàn):

1)葉片迎風(fēng)面所受風(fēng)壓為正壓(壓力面)，而背風(fēng)面基本處于負(fù)壓(吸力面).由伯努利效應(yīng)得知:流速越快，壓力越低;流速越慢，壓力越高.因而葉片背風(fēng)面風(fēng)速大于迎風(fēng)面風(fēng)速.葉片迎風(fēng)面壓強(qiáng)高于大氣壓產(chǎn)生壓力，背風(fēng)面壓強(qiáng)低于大氣壓產(chǎn)生吸力，由此對(duì)翼型產(chǎn)生升力，這也是葉片能夠旋轉(zhuǎn)的原因.

2)葉片迎風(fēng)面，沿著翼型曲線前緣至后緣，壓力變化平緩，而背風(fēng)面壓力變化迅速.產(chǎn)生這種分布的原因是翼型截面曲率導(dǎo)致的，當(dāng)來流流經(jīng)翼型表面時(shí)，翼型幾何特性引起了翼型表面來流風(fēng)速的變化，因而造成了壓力分布的相應(yīng)變化.

圖6 葉片壓力分布云圖Fig.6 Pressure distribution on turbine blade

圖7為葉表面的速度矢量圖，可以觀察到翼型繞流現(xiàn)象，且背風(fēng)面風(fēng)速大于迎風(fēng)面風(fēng)速.

圖7 葉片速度矢量圖Fig.7 velocity vector distributions on turbine blade

圖8為葉尖速比－風(fēng)能利用系數(shù)特性，與GHBladed軟件計(jì)算結(jié)果比較可以發(fā)現(xiàn):兩者計(jì)算結(jié)果較為接近，趨勢(shì)一致，但理論值略小.這是由于本文CFD計(jì)算采用粘流RANS方程，而GHBladed軟件基于葉素理論，因此在旋轉(zhuǎn)方向的阻力分量的理論計(jì)算較軟件大，得到轉(zhuǎn)矩比軟件小，因此CFD計(jì)算得到的葉片吸收功率比GHBladed計(jì)算值要小.圖9為葉尖速比－推力系數(shù)特性，理論計(jì)算結(jié)果與軟件計(jì)算結(jié)果相當(dāng)，理論值略小，這是由于推力分量中，升力貢獻(xiàn)占主要部分，阻力分量的差異導(dǎo)致軟件計(jì)算結(jié)果的偏大.

圖8 葉尖速比-風(fēng)能利用系數(shù)特性Fig.8 Characteristics of TSR-Cp

為了減小尾流的誘導(dǎo)損失，風(fēng)葉片環(huán)量的設(shè)計(jì)分布一般采用葉根和葉梢卸載，將載荷盡力均分到葉片中間區(qū)域［8］.圖10為2 MW風(fēng)機(jī)的葉片環(huán)量無因次化后的徑向分布，它反映了葉片徑向載荷的變化趨勢(shì)，葉片載荷最大位置在0.4R處.從圖中可以發(fā)現(xiàn)，2 MW風(fēng)機(jī)葉片的環(huán)量分布還是比較合理的.為了將其徑向環(huán)量更均勻地分布到葉片中間區(qū)域，作者調(diào)整了葉片中間區(qū)域的弦長和扭曲角，使葉片環(huán)量在中間分布得更均勻，計(jì)算結(jié)果表明優(yōu)化方案在設(shè)計(jì)尖速比附近的Cp比原型提高了3%左右.

圖9 葉尖速比-推力系數(shù)特性圖Fig.9 Characteristics of TSR～CT

圖10 葉片環(huán)量無因次化后的徑向分布Fig.10 Distribution of circulation along radii of wind turbine blade

4 結(jié)論

本文基于CFD技術(shù)，建立了大型風(fēng)力發(fā)電機(jī)組葉片氣動(dòng)性能分析方法，通過算例分析與比較，得到了以下結(jié)論:

1)小攻角下的二維氣動(dòng)性能分析與試驗(yàn)結(jié)果吻合良好，大攻角下的氣動(dòng)性能計(jì)算結(jié)果趨勢(shì)合理，證明本文建立的二維翼型氣動(dòng)性能分析方法是可行的，可為風(fēng)機(jī)分析軟件提供重要的氣動(dòng)輸入數(shù)據(jù);

2)2MW風(fēng)機(jī)葉輪的三維氣動(dòng)性能計(jì)算結(jié)果與GHBladed軟件吻合良好，證明本文建立的三維翼型氣動(dòng)性能分析方法也是可行的.壓力分布與速度矢量分布顯示三維方法可以更真實(shí)地模擬流場(chǎng)，考慮葉片的三維效應(yīng).對(duì)2MW風(fēng)機(jī)葉片翼型的優(yōu)化，進(jìn)一步體現(xiàn)了三維方法的優(yōu)勢(shì).

因此，基于CFD技術(shù)，采用RANS方程耦合SST湍流模型，預(yù)報(bào)二維和三維翼型的氣動(dòng)性能是可行的.本文研究成果對(duì)于大型風(fēng)力發(fā)電機(jī)組葉片氣動(dòng)性能的設(shè)計(jì)與優(yōu)化具有一定的參考價(jià)值.

［1］劉萬餛，張志英，李銀鳳.風(fēng)能與風(fēng)力發(fā)電技術(shù)［M］.北京:化學(xué)工業(yè)出版社，2006:1-24.

［2］SORENSEN N，MICHELSEN J.Aerodynamic predictions for the unsteady aerodynamics experiment phase II rotor at the national renewable energy laboratory［EB/OL］.［1999-01-11］.http://rotorcraft.arc.nasa.gov/publications/files/ S?rensen_AIAA1999.pdf.

［3］SORENSEN N，MICHELSEN J，SCHRECK S.Navier-Stokes predictions of the NREL phase VI rotor in the NASA Ames 80-by-120 wind tunnel［C］//ASME 2002 Wind Energy Symposium.Reno，USA，2002:94-105.

［4］李仁年，李銀然，王秀勇，等.風(fēng)力機(jī)翼型的氣動(dòng)模型及數(shù)值計(jì)算［J］.蘭州理工大學(xué)學(xué)報(bào)，2010，36(3):65-68.

LI Rennian，LI Yinran，WANG Xiuyong，et al.Aerodynamic model of airfoil for wind turbine and its numeric computation［J］.Journal of Lanzhou University of Technology，2010，36 (3):65-68.

［5］MADSEN H，SORENSEN N，SCHRECK S.Yaw aerodynamics analyzed with three codes in comparison with experiment［C］//ASME 2003 Wind Energy Symposium.Reno，USA，2003:94-103.

［6］MENTER F R.Multiscale model for turbulent flows［C］// AIAA 24th Fluid Dynamics Conference.American Institute of Aeronautics and Astronautics.Orlando，USA，1993:1-21.

［7］ABBOTT I H，VON DOENHOFF A E.Theory of wing sections［M］.New York:Dover Publications，INC.，1959: 428-430.

［8］TONY B.風(fēng)能技術(shù)［M］.北京:科學(xué)出版社，2007:74-76.