分裂葉尖概念型風力機葉片的氣動設計與數值優化研究

張震宇，王同光，陳 立，許波峰，王 瓏，羅 源

（1．南京航空航天大學 江蘇省風力機設計高技術研究重點實驗室，江蘇 南京 210016；2．中國空氣動力研究與發展中心 空氣動力學國家重點實驗室，四川 綿陽 621000）

0 引 言

水平軸風力機的風能輸出受到來流風況的強烈影響，在偏航以及中等、大風速條件下由于葉片背風面的邊界層容易發生流動分離，并會導致扭矩輸出和能量利用效率的明顯降低，以及分離噪聲等一系列影響。風力機葉片的上述特性通過風洞實驗和計算流體力學數值模擬不斷被證實［1－2］。

另一方面，水平軸風力機的大展弦比布局一般會產生明顯的葉尖渦結構；在三維旋轉主流流動的作用下，螺旋狀發展的葉尖渦會對葉片產生比較強的誘導影響［3］。

如何準確捕捉風力機葉片上述流動結構是實現合理高效氣動設計的前提。本文采用數值方法，獲得了以NREL phase VI風力機葉片為原型的流場與葉尖渦結構的準確模擬結果，并與縮比模型的PIV實驗結果形成了良好的對比。在此基礎上，提出了基于多片小翼概念的分裂葉尖形式概念型水平軸風力機葉片設計方案，并針對其綜合氣動性能利用快速非支配排序遺傳算法（NSGA－II）進行了多目標優化研究并獲得了滿足給定約束的優化結果。

1 數值模擬

1．1 模型與數值方法

水平軸風力機葉片流動具有明顯的三維旋轉和強烈的湍流特征。基于工程允許的計算能力，除了大渦模擬（LES）與脫落渦模擬（DES）方法外，數值求解雷諾平均Navier－Stokes（RANS）方程仍然是目前工程上能夠兼顧湍流流動計算量和計算準確度的一種有效方法。本文采用該方法實現了風力機葉片基本流場的數值模擬。

由于暫不考慮側滑等復雜風況，通過在上述區域中采用旋轉坐標系求解三維不可壓RANS方程，原始的非定常葉片繞流可以轉化為新坐標系下的準定常流動進行求解。其控制方程形式在常規笛卡爾坐標系下的形式基礎上應補充附加源項以考慮在非慣性系下的離心力與哥氏力影響。

根據葉片繞流的周向對稱性，計算區域采用一半流場并根據低速不可壓流動特性分別施加速度進口、壓力出口，無滑移葉片物面以及周期邊界條件，如圖1所示。

圖1 風力機葉片流場計算區域Fig．1 Domain for the half flow field of a wind turbine blade

計算區域大小沿流向（圖1中－z方向）為（－2R，8R），沿展向為5R（R為葉片長度）。總網格單元數超過6百萬。為了在葉片下游區域計算網格局部加密以捕捉尾渦結構，總網格單元數接近2千2百萬。

上述控制方程的離散形式通過有限體積法進行求解。其對流項采用二階迎風格式，擴散項為中心差分。壓力與速度的耦合采用SIMPLE算法實現［4］。RANS方程的封閉采用 Menter的k－ωSST模型［5］。

1．2 結果與討論

利用上述方法對于原始NREL phase VI風力機葉片進行了計算，并與NREL盲比實驗以及文獻［2］中S?rensen等人的計算結果進行了對比，見圖2。

圖2 低速軸扭矩結果對比Fig．2 Comparison of low－speed shaft torque

利用目前方法，尾流加密（22mil）和不加密（6mil）兩套網格算得的低速軸扭矩結果在主要風速范圍內與實驗存在良好一致性。除10m／s之外主要計算結果與實驗相對誤差均在6%之內，說明對于該風力機葉片的整體性能把握較好。10m／s風速下出現較大偏差的主要原因是此風況下葉片吸力面邊界層流動未完全分離，而目前所用RANS算法內尚未引入對于流動轉捩位置的預估功能，因此導致該風速下對吸力面分離流動的影響把握不準。

圖3 尾流區渦量強度分布（周向角60°，風速12m／s，尖速比4．91）Fig．3 Vorticity in the wake（azimuth angle of 60°，wind speed of 12m／s，and tip speed ratio of 4．91）

在此基礎上本文利用尾流區加密網格準確捕捉了近尾跡區葉尖渦的發展過程。圖3中為與風洞PIV實驗結果的對比。該實驗設置的詳細介紹見文獻［6］。由圖3可見，在給定風速下，計算得到的尾渦渦核位置與渦核強度與實驗非常一致。

2 分裂葉尖布局設計

2．1 翼尖小翼概念

長期以來關于固定翼飛行器的研究表明，翼尖小翼能夠有效的起到提高升力，抑制誘導阻力的效果。在此基礎上基于多個小翼的分裂葉尖形式證明在機翼增升減阻方面具有更為優良的效果［7］。小翼設計在風力機葉片的設計中在近年來也得到了應用［8］。本文提出了基于多片小翼概念的分裂葉尖形式概念型水平軸風力機葉片設計方案，以研究該布局對于風力機葉片氣動性能的控制效果。其基本形式如圖4所示。

2．2 分裂葉尖數值模擬

作為初步設計，在原始葉片85%展長位置改為由三段相互獨立的小翼組成的分裂葉尖布局。各小翼均無扭轉，梢根比均為0．75，后掠角均為45°，從前緣向后緣上反角依次為40°，20°，0°。

圖4 帶三個小翼的NREL phase VI葉片模型Fig．4 NREL phase VI model with three winglets

由于采用分裂葉尖形式的氣動布局，多個小翼將原來集中于葉尖附近的環量進行了重新分配。圖5中所示下游近尾跡區渦核附近的渦量極值隨尾跡周向角的分布情況描述了這種重新分配環量的效果。分裂葉尖的布局（見圖5中winglet1結果）使得近尾跡區渦核附近的渦量極值比原始葉片（見圖5中basic結果）大為降低，從而有效減弱下游集中渦的誘導作用。

圖5 尾流區渦量極值沿周向角的分布Fig．5 Azimuthal distribution of vorticity maximum in the wake

另一方面，各小翼的厚度均接近6%，且帶有大后掠角，這種設計使得在中等和大風速下葉片吸力面在前緣附近獲得更大的吸力峰（如圖6中winglet1黑線所示），并體現為相應風速下軸扭矩輸出的明顯增長，在10m／s－25m／s范圍內，增長幅度在7%－40%之間，見圖7中winglet1結果。

圖6 分裂葉尖布局80%展長位置壓力分布，25m／sFig．6 Pressure distribution at 80%span，25m／s

圖7 分裂葉尖布局低速軸扭矩Fig．7 Low－speed shaft torque of the blade with multiple winglets

3 數值優化

在上述工作基礎上，對葉片關鍵幾何參數利用快速非支配排序遺傳算法（NSGA－II）進行了多目標數值優化。為問題簡化考慮，優化風速定為20m／s，小翼根稍比均為0．75。

該優化問題中需對以下兩目標函數最小化：

其中Cp為風能利用系數，CT為軸向推力系數。由于CFD方法所需計算量過大，優化部分的上述目標函數采用自由渦尾跡方法計算得出，具體算法見文獻［10］。

優化變量及約束條件為：

小翼上反角（向吸力面方向為正）：

小翼扭角：

3．1 遺傳算法

快速非支配排序遺傳算法（Fast and Elitist Nondominated Sorting Genetic Algorithm，NSGAII）是Deb等人于2002年在對NSGA算法進行了大幅改進的基礎上提出的一種高效分類排序算法［9］。該算法引入了Pareto支配關系和個體分布評價函數，從而可以將復雜多目標優化問題轉化成可比較個體優劣的單目標優化問題，能夠用于進行任意目標個數的數值優化工作。該算法根據Pareto支配關系對種群內個體首先進行非支配分層以確定各個個體的層級；再用擁擠距離評價函數控制各個層級個體分布；然后用個體的層級個體之間距離做選擇操作，并結合傳統遺傳算法的其它操作過程，最終得到Pareto最優解集。

3．2 結果與討論

優化種群大小選取為36，迭代次數300次。得到的Pareto最優解集分布如圖8所示。利用自由渦尾跡方法算得的基本型NREL PHASE VI葉片的風能利用系數Cp與軸向推力系數CT分別為0．02286和0．1635。比較發現，目前得到的優化解風能利用系數分布均明顯高于該參考值，相對提高27%～31%不等；同時軸向推力系數CT與NREL PHASE VI原型葉片相比變化不大。

圖8 Pareto最優解集分布Fig．8 Pareto solution sets

由NSGA－II算法得到的上述最優解集在給定目標函數空間中給出了一組滿足相應約束條件的最優解，由于該解集任一個最優解并不更優于其他解，因此在最終的優化方案取決于研究者根據問題性質的具體分析與著眼點。由圖8得知，在目前的優化范圍內，隨著風能利用系數的提高，一般伴隨有軸向推力的增長。選取三個典型狀態如表1所示，并分別對應于軸向推力最小、折中、能量轉化率最大三個考慮因素。

4 結 論

借助于準確的數值模擬方法得知，分裂葉尖形式氣動布局能夠有效地分散風力機葉尖渦量的分布，從而減小其誘導影響；通過合理布置小翼，能夠使葉片端部區域在中等和大風速下獲得額外的前緣吸力，實現更高的氣動效率，并為進一步在此風速范圍的氣動控制策略提供可行的技術準備。

基于NSGA－II遺傳算法的優化方法能夠在給定的性能約束范圍內給出準確的全局最優解分布。對于上述問題，在不明顯增大軸向推力的前提下，以風能利用系數形式給出的氣動效率得到明顯改善，比原始設計提高達30%左右。以上討論表明，NSGA－II遺傳算法與自由渦尾跡方法結合，是一種行之有效的風力機葉片氣動性能優化方法。

［1］ 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］．NREL／TP－500－29955，NREL，2001．

［2］ S?RENSEN N N，MICHELSEN J，SCHRECK S．Navier－Stokes predictions of the NREL phase VI rotor in the NASA Ames 80ft×120ft wind tunnel［J］．Wind Energy，2002，5（2－3）：151－169．

［3］ VERMEER L J，SORENSEN J N，CRESPO A．Wind turbine wake aerodynamics［J］．Progress in Aerospace Sciences，2003，39（6－7）：467－510．

［4］ VERSTEEG H K，MALALASEKERA W．An introduction to computational fluid dynamics－the finite volume method［M］．2nd Ed．，Pearson Education Limited，2007

［5］ MENTER FR．Two－equation eddy－viscosity turbulence models for engineering applications［J］．AIAA Journal，1994，32（8）：1598－1605．

［6］ 肖京平，武杰，陳立，史喆羽，風力機葉尖渦尾跡結構PIV測量實驗研究［J］．應用數學和力學，2011，32（6）：683－692．

［7］ 陳明巖，齊孟卜，翼尖帆片的增升減阻研究［J］．航空學報，1994，15（6）：641－646．

［8］ JOHANSEN J，SORENSEN N N．Aerodynamic investigation of winglets on wind turbine blades using CFD［R］．Ris?－R－1543（EN），2006．

［9］ DEB K，PRATA P A，AGARWAL S，et al．A fast and elitist multi－objective genetic algorithm：NSGA－Ⅱ［J］．IEEE Transactions on Evolutionary Computation，2002，6（2）：182－197．

［10］許波峰，王同光．基于自由渦尾跡法和面元法全耦合風力機氣動特性計算［J］．南京航空航天大學學報，2011，43（5）：592－597．