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

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

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

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

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

1.1 控制方程

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

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

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

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

連續(xù)性方程:

動量方程:

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

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

1.2 SSTk-ω湍流模型

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

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

2 二維翼型氣動性能分析

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

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

2.1.1 坐標系定義

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

圖1 坐標系的定義Fig.1 Definition of coordinates

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

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

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

2.2 二維翼型氣動性能算例分析

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

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

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

整個計算域為進流段和去流段都為10倍弦長，而側(cè)面為6倍弦長，具體為10L×6L×10L，其中L為弦長.邊界條件為:

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

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

3)外場邊界:外場邊界的設(shè)置與進口邊界條件一致;

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

2.1.4 數(shù)值計算方法

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

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

3 三維葉輪氣動性能分析

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

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

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

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

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

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

3.1.2 邊界條件和初始條件

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

3.1.3 數(shù)值計算方法

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

3.2 三維葉輪氣動性能算例分析

3.2.1 對象描述

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

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

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

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

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

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

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

2)葉片迎風(fēng)面，沿著翼型曲線前緣至后緣，壓力變化平緩，而背風(fēng)面壓力變化迅速.產(chǎn)生這種分布的原因是翼型截面曲率導(dǎo)致的，當(dāng)來流流經(jīng)翼型表面時，翼型幾何特性引起了翼型表面來流風(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軟件計算結(jié)果比較可以發(fā)現(xiàn):兩者計算結(jié)果較為接近，趨勢一致，但理論值略小.這是由于本文CFD計算采用粘流RANS方程，而GHBladed軟件基于葉素理論，因此在旋轉(zhuǎn)方向的阻力分量的理論計算較軟件大，得到轉(zhuǎn)矩比軟件小，因此CFD計算得到的葉片吸收功率比GHBladed計算值要小.圖9為葉尖速比－推力系數(shù)特性，理論計算結(jié)果與軟件計算結(jié)果相當(dāng)，理論值略小，這是由于推力分量中，升力貢獻占主要部分，阻力分量的差異導(dǎo)致軟件計算結(jié)果的偏大.

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

為了減小尾流的誘導(dǎo)損失，風(fēng)葉片環(huán)量的設(shè)計分布一般采用葉根和葉梢卸載，將載荷盡力均分到葉片中間區(qū)域［8］.圖10為2 MW風(fēng)機的葉片環(huán)量無因次化后的徑向分布，它反映了葉片徑向載荷的變化趨勢，葉片載荷最大位置在0.4R處.從圖中可以發(fā)現(xiàn)，2 MW風(fēng)機葉片的環(huán)量分布還是比較合理的.為了將其徑向環(huán)量更均勻地分布到葉片中間區(qū)域，作者調(diào)整了葉片中間區(qū)域的弦長和扭曲角，使葉片環(huán)量在中間分布得更均勻，計算結(jié)果表明優(yōu)化方案在設(shè)計尖速比附近的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ā)電機組葉片氣動性能分析方法，通過算例分析與比較，得到了以下結(jié)論:

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

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

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

［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)力機翼型的氣動模型及數(shù)值計算［J］.蘭州理工大學(xué)學(xué)報，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.