史先路，趙 凡，趙 凱，，趙軍華，3

(1.江南大學(xué) 機(jī)械工程學(xué)院，江蘇 無錫 214122；2.中國空氣動(dòng)力研究與發(fā)展中心，四川 綿陽 621000；3.江蘇省食品先進(jìn)制造裝備技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 無錫 214122）

0 引言

風(fēng)力機(jī)常年運(yùn)行在惡劣環(huán)境下，葉片表面易表現(xiàn)出流動(dòng)分離、失速等非定常氣動(dòng)特性，設(shè)計(jì)新一代高效、輕量化的風(fēng)力渦輪機(jī)，需要更準(zhǔn)確地預(yù)測氣動(dòng)載荷。

作為風(fēng)力機(jī)葉片設(shè)計(jì)及其氣動(dòng)性能評(píng)估的常用方法，葉素動(dòng)量理論(BEM）[1]由于其計(jì)算簡單，能快速預(yù)測葉片氣動(dòng)性能而備受青睞。文獻(xiàn)[2]利用BEM理論分析了風(fēng)力機(jī)葉片氣動(dòng)載荷的分布規(guī)律以及風(fēng)速大小對(duì)葉片氣動(dòng)載荷特性的影響。作為一種常用的數(shù)值模擬方法，計(jì)算流體動(dòng)力學(xué)(CFD）可用于解決所有常見的流動(dòng)問題。文獻(xiàn)[3]，[4]利用CFD方法在均勻來流條件下對(duì)風(fēng)力機(jī)葉片表面流動(dòng)分離情況進(jìn)行了研究，但未考慮大氣邊界層的影響。文獻(xiàn)[5]基于CFD計(jì)算方法研究了一種小型風(fēng)力機(jī)的氣動(dòng)特性，將計(jì)算結(jié)果與基于BEM計(jì)算方法得到的結(jié)果進(jìn)行了對(duì)比，結(jié)果表明，兩種計(jì)算方法得到的氣動(dòng)載荷特性相似。文獻(xiàn)[6]基于BEM與CFD計(jì)算方法研究了風(fēng)力機(jī)的功率和載荷特性，采用Spalart-Allmaras湍流模型針對(duì)海上變速變槳風(fēng)力機(jī)分析了不同雷諾數(shù)下兩種計(jì)算方法的結(jié)果差異，但未考慮失速狀態(tài)下葉片表面的流動(dòng)分離對(duì)計(jì)算結(jié)果的影響。

研究人員利用BEM及CFD數(shù)值模擬方法對(duì)風(fēng)力機(jī)葉片的氣動(dòng)特性進(jìn)行了較為深入的研究，但是，基于BEM與CFD兩種計(jì)算方法對(duì)風(fēng)力機(jī)葉片氣動(dòng)載荷的對(duì)比研究較少，且集中在兩葉片風(fēng)力機(jī)或者功率較小的小型風(fēng)力機(jī)。本文采用修正的BEM和不同來流條件下CFD數(shù)值模擬兩種方法研究大型風(fēng)力機(jī)三維旋轉(zhuǎn)效應(yīng)導(dǎo)致的葉尖損失、動(dòng)態(tài)失速等現(xiàn)象，分析對(duì)比兩種風(fēng)力機(jī)氣動(dòng)性能計(jì)算方法，為風(fēng)力機(jī)的設(shè)計(jì)、功率預(yù)測和氣動(dòng)載荷的研究以及BEM的進(jìn)一步修正提供理論支撐。

1 氣動(dòng)性能計(jì)算方法

1.1 修正后的葉素動(dòng)量理論

BEM通過對(duì)作用在葉片展向上每個(gè)微段上的力進(jìn)行積分得到作用于風(fēng)輪上的推力及轉(zhuǎn)矩。圖1為葉片剖面和氣流角受力關(guān)系圖。圖中：入流角φ為總速度W與風(fēng)輪轉(zhuǎn)動(dòng)平面的夾角，由局部攻角(α）和葉素槳距角(β）組成，變量V1為入流風(fēng)速，Ω為風(fēng)輪轉(zhuǎn)速，a為軸向誘導(dǎo)因子，b為周向誘導(dǎo)因子。根據(jù)作用在葉片截面上的升力(L）和阻力(D）的關(guān)系可以得到葉素上氣動(dòng)合力的軸向分量和切向分量。

圖1 葉片剖面和氣流角受力關(guān)系Fig.1 Force relation diagram of blade profile and flow angle

作用在葉片單元上的法向力和轉(zhuǎn)矩可以通過升力系數(shù)(CL）和阻力系數(shù)(CD）的關(guān)系得到，假設(shè)葉片單元的弦長(c）已知，葉片數(shù)為N，則法向力dFA和轉(zhuǎn)矩dT可以表示為

由于模型中假設(shè)每個(gè)葉片截面與相鄰截面之間沒有相互作用，使得該模型無法捕捉葉片上發(fā)生的三維效應(yīng)，尤其是在葉尖和葉根附近。此外，氣流繞風(fēng)輪葉片剖面流動(dòng)使得葉片表面氣流分離，進(jìn)而導(dǎo)致梢部和根部翼型表面的環(huán)量減少，相應(yīng)計(jì)算的轉(zhuǎn)矩也會(huì)減小。因此，利用BEM計(jì)算時(shí)需要考慮修正因子，通常使用文獻(xiàn)[7]的葉根、葉尖損失修正因子來校正剖面數(shù)據(jù)。

式中：rn為輪轂半徑；Ft為葉尖損失修正因子；Fr為葉根損失修正因子。

令葉尖葉根修正因子F=Ft?Fr，則修正后的法向力和轉(zhuǎn)矩可以表示為

文獻(xiàn)[8]在文獻(xiàn)[7]的修正因子的基礎(chǔ)上，引入一個(gè)新的修正項(xiàng)用于法向力系數(shù)Cn和切向力系數(shù)Ct，并將其代入經(jīng)典的葉素理論中得到推力與轉(zhuǎn)矩。

同 理，令F′=Ft′?F′，其 中g(shù)=exp[-c1(Nλ-c2）]，c1=0.125，c2=21，則修正后的法向力和轉(zhuǎn)矩可以表示為

1.2 CFD數(shù)值模擬

1.2.1建立計(jì)算模型及網(wǎng)格無關(guān)性驗(yàn)證

以某1.5 MW大型水平軸風(fēng)力機(jī)為研究對(duì)象，選取風(fēng)機(jī)葉片常用翼型S809，設(shè)計(jì)功率為1.5 MW，風(fēng) 速V為12 m/s，風(fēng) 輪 轉(zhuǎn) 速 為22.6 r/min，直徑D為71 m。設(shè)置內(nèi)域?yàn)樾D(zhuǎn)域，外域?yàn)殪o止域，劃分非結(jié)構(gòu)化四面體網(wǎng)格并對(duì)葉片附近以及內(nèi)流域區(qū)域進(jìn)行網(wǎng)格加密，當(dāng)網(wǎng)格數(shù)量大于960萬時(shí)，功率誤差變化很小，因此最終確定數(shù)值計(jì)算的總網(wǎng)格數(shù)為960.4萬。網(wǎng)格無關(guān)性驗(yàn)證見表1。

表1 網(wǎng)格無關(guān)性驗(yàn)證Table 1 Validation of grid independence

將模型導(dǎo)入求解器中進(jìn)行求解，定義湍流模型為SST k-ω，設(shè)置速度入口和壓力出口。壓力速度耦合采用SIMPLE算法，離散方式為二階迎風(fēng)格式。將速度入口邊界條件以UDF的形式加載到求解器中分別進(jìn)行均勻風(fēng)場和非定常流動(dòng)風(fēng)場的風(fēng)力機(jī)三維數(shù)值模擬，在葉尖和葉根處設(shè)置壓力監(jiān)測點(diǎn)，監(jiān)測葉片表面壓力隨時(shí)間的變化。壓力監(jiān)測點(diǎn)位置如圖2所示。

圖2 壓力監(jiān)測點(diǎn)位置Fig.2 Location of the pressure monitoring point

1.2.2湍流風(fēng)場的模擬方法

由于風(fēng)力機(jī)葉片在旋轉(zhuǎn)過程中會(huì)發(fā)生流動(dòng)分離、失速延遲等現(xiàn)象，而復(fù)雜的工作環(huán)境，例如大氣邊界層、湍流、風(fēng)剪切、風(fēng)速和風(fēng)向的瞬變性以及偏航等會(huì)使翼型及葉片的氣動(dòng)行為呈現(xiàn)強(qiáng)烈的非定常特性，目前大多數(shù)關(guān)于水平軸風(fēng)力機(jī)的研究沒有考慮到這種不穩(wěn)定性。

非定常流的模型復(fù)雜，模擬真實(shí)脈動(dòng)風(fēng)需要大量計(jì)算資源與時(shí)間，在理想情況下，可以將高頻湍流風(fēng)視為均勻來流、剪切來流以及湍流，其中湍流風(fēng)可以被理想化為隨時(shí)間變化的具有平均值、幅值和頻率的正弦曲線形式的風(fēng)速序列。此外，文獻(xiàn)[9]的研究結(jié)果表明，正弦函數(shù)形式的脈動(dòng)風(fēng)速序列非常適合替代湍流來研究非定常流動(dòng)行為，因此，本文使用式(11）表示某一局部時(shí)間點(diǎn)的瞬時(shí)風(fēng)速。

式 中：vˉ(z）為 距 地 面 高 度z處 的 平 均 風(fēng) 速；k為 湍流強(qiáng)度，為風(fēng)速的標(biāo)準(zhǔn)差與平均風(fēng)速之比；f為風(fēng)的脈動(dòng)頻率；t為流動(dòng)時(shí)間。

為了突出脈動(dòng)風(fēng)的波動(dòng)特征，對(duì)式(11）中的相應(yīng)參數(shù)進(jìn)行指定，來設(shè)置較高的頻率和振幅，其中f取1 Hz，k取0.1，z垂直于地面方向向上。

平均風(fēng)速隨高度的變化可以用指數(shù)率分布來描述。

式 中：γ為 切 變 系 數(shù) 或 粗 糙 度 指 數(shù)，取0.27；vˉ(z0）為距地面高度z0=10 m處的平均風(fēng)速，取6.4 m/s。

2 結(jié)果與討論

2.1 翼型氣動(dòng)性能分析

在進(jìn)行三維數(shù)值模擬前，先分析S809翼型的氣動(dòng)性能，以驗(yàn)證湍流模型的準(zhǔn)確性，俄亥俄州立大學(xué)(OSU）對(duì)翼型S809進(jìn)行了不同雷諾數(shù)條件下的風(fēng)洞試驗(yàn)[10]。本文為了改善失速狀態(tài)下數(shù)值模擬的準(zhǔn)確性，在均勻來流條件下使用SST k-ω湍流模型，選取封閉常數(shù) β*作為調(diào)整參數(shù)。β*會(huì)改變湍動(dòng)能k和湍流耗散率ω，在翼型數(shù)值模擬中校正湍流模型的封閉常數(shù)并用于葉片的數(shù)值模擬。獲得不同攻角(α）下翼型的CL和CD，并將仿真結(jié)果與OSU的二維靜態(tài)測試數(shù)據(jù)和BEM計(jì)算結(jié)果進(jìn)行對(duì)比(圖3）。

圖3 不同攻角下翼型升力系數(shù)對(duì)比Fig.3 Comparison of airfoil lift coefficients at different angles of attack

由圖3可知：在 α為0～8 °時(shí)，CL呈線性增長，在 α為8～16°時(shí)，CL趨于平緩，在 α超過16°后，CL急劇下降，調(diào)整封閉參數(shù)后的CFD模擬結(jié)果和實(shí)驗(yàn)數(shù)據(jù)幾乎完全吻合；在α>8°時(shí)，BEM計(jì)算結(jié)果與CFD模擬結(jié)果和實(shí)驗(yàn)結(jié)果差異較大，表明BEM不能很好地預(yù)測翼型在失速狀態(tài)下的氣動(dòng)性能，需要進(jìn)一步修正。

圖4為不同 α下的翼型流線圖。由圖4可知，隨著 α的增大，失速區(qū)域逐漸增大且分離點(diǎn)逐漸向前緣移動(dòng)，可將3個(gè)階段稱為附著流動(dòng)狀態(tài)(α為0～8°）、輕 失 速 狀 態(tài)(α為8～16°）和 深 度失 速 狀 態(tài)(α>16°）。

圖4 不同攻角下翼型流線圖Fig.4 Flow diagram of airfoil at different angles of attack

2.2 葉片表面極限流線圖

圖5為均勻來流條件下不同葉尖速比(TSR）下的極限流線。當(dāng)風(fēng)速分別為12，14，16.8 m/s和21 m/s時(shí)，對(duì) 應(yīng) 的TSR分 別 為7，6，5和4。在 離 心力和葉片展向壓力梯度的作用下，氣流在葉根處發(fā)生分離并沿葉片展向流動(dòng)，展向壓力梯度隨著攻角和葉尖速比的變化而變化。

圖5 葉片表面極限流線Fig.5 Limiting streamline of blade surface

由圖5可知：當(dāng)風(fēng)速為12 m/s時(shí)，所有流線均從前緣指向后緣，葉片壓力面均為附著流動(dòng)，沒有發(fā)生失速，除葉根小部分區(qū)域外，葉片吸力面基本未發(fā)生分離流動(dòng)；隨著風(fēng)速的增加，葉根部位分離點(diǎn)向前緣移動(dòng)，葉片吸力面分離區(qū)域逐漸增大，葉片展向有明顯的流動(dòng)分離線；當(dāng)風(fēng)速增大到21 m/s時(shí)，葉片處于深度失速狀態(tài)，吸力面流動(dòng)完全分離，沿葉片展向有速度分量，在靠近葉根部位產(chǎn)生一個(gè)較大的氣流漩渦。

2.3 葉片表面壓力頻譜特性分析

為了進(jìn)一步研究非定常來流對(duì)風(fēng)力機(jī)葉片表面壓力波動(dòng)的影響，對(duì)各葉尖和葉根部位截面監(jiān)測點(diǎn)所測壓力進(jìn)行FFT變換，得到壓力功率譜(圖6）。壓力功率譜總體上符合Kolmogorov-5/3功率定律，所有截面在頻率f=0.376 Hz時(shí)均出現(xiàn)明顯的峰值，而且該頻率與風(fēng)輪的旋轉(zhuǎn)頻率相等，風(fēng)輪的旋轉(zhuǎn)頻率fn=n/60(n為風(fēng)輪轉(zhuǎn)速），說明壓力脈動(dòng)的主頻為1倍風(fēng)輪轉(zhuǎn)動(dòng)頻率，即f=1fn。

圖6 翼型壓力的功率譜特性Fig.6 Power spectrum characteristics of airfoil pressure

由圖6(a）可知：當(dāng)f=1fn時(shí)，所有測點(diǎn)壓力功率譜均出現(xiàn)明顯的峰值，峰值從大到小依次是x/c等 于0.85，0.11，0.32和0.57，靠 近 后 緣 的x/c為0.85時(shí)的功率譜峰值最大，說明葉片后緣壓力含有的主頻能量最高；在高頻段，靠近尾緣的x/c=0.85呈現(xiàn)明顯的鋸齒狀分布特性，其余所有壓力功率譜均呈現(xiàn)紊亂狀態(tài)，表明葉尖部位翼型吸力面對(duì)失速旋渦的敏感程度較高，尤其是尾緣部位。

由圖6(b）可知：當(dāng)f=1fn時(shí)，所有測點(diǎn)壓力功率譜均出現(xiàn)明顯的峰值，峰值從大到小依次是x/c等于0.11，0.32，0.85和0.57，最接近前緣的x/c等于0.11的功率譜峰值最大，表明前緣對(duì)來流最敏感，后緣次之而翼型中部最弱；在頻率f為2～11 Hz時(shí)，功率譜峰值點(diǎn)不明顯，表明葉尖部位翼型壓力面對(duì)失速旋渦的敏感程度較低。

由于葉根部位翼型受轉(zhuǎn)頻影響較小，使得葉根在2～5倍轉(zhuǎn)頻的功率譜峰值不明顯[圖6(c），(d）]，但總體上符合Kolmogorov-5/3功率定律。

總體而言，壓力功率譜在高頻段均呈現(xiàn)相對(duì)平緩的變化趨勢，而且壓力面壓力功率譜的平緩變化趨勢更明顯，翼型越接近葉根，平緩變化的功率譜所占頻段越寬；由于靠近葉根翼型的攻角較大，吸力面受到順壓梯度作用較強(qiáng)，流動(dòng)分離嚴(yán)重，使得吸力面靠近葉根部位的功率譜的平緩變化區(qū)域較小。

2.4 BEM與CFD葉片載荷預(yù)測結(jié)果對(duì)比

本文使用QBlade進(jìn)行BEM計(jì)算來預(yù)測葉片氣動(dòng)載荷，考慮了失速延遲效應(yīng)的影響，分別利用兩種修正方法對(duì)BEM進(jìn)行修正并將修正后的結(jié)果與CFD計(jì)算結(jié)果進(jìn)行對(duì)比。將葉片弦長、扭角等參數(shù)導(dǎo)入到QBlade中完成葉片模型的建立并進(jìn)行氣動(dòng)性能分析，然后與文獻(xiàn)[11]中的CFD模擬計(jì)算結(jié)果進(jìn)行對(duì)比(圖7）。

圖7 葉片截面力對(duì)比Fig.7 Comparison of blade cross section force

由圖7可知：兩種不同風(fēng)力機(jī)模型計(jì)算得到的氣動(dòng)載荷趨勢一致，沿著展向葉片所受Fn逐漸增大，葉尖損失的存在使得在接近葉尖時(shí)Fn和Ft均有所減小，同時(shí)也從側(cè)面證明了本文葉片設(shè)計(jì)與計(jì)算結(jié)果的合理性；Fn和Ft在葉根處修正后的BEM與原始的BEM計(jì)算結(jié)果差別不大，而在葉尖處修正過后更加接近CFD模擬結(jié)果；相較于文獻(xiàn)[7]修正的BEM，文獻(xiàn)[8]修正的BEM對(duì)于Fn的預(yù)測更加準(zhǔn)確，對(duì)于Ft，二者預(yù)測結(jié)果差別不大。整體而言，修正后的BEM能夠較好地預(yù)測出氣動(dòng)載荷隨著葉片展向位置變化的趨勢，但是與CFD計(jì)算結(jié)果相比仍存在一定偏差。

3 結(jié)論

本文采用修正的BEM和三維CFD數(shù)值模擬兩種計(jì)算方法分析預(yù)測了1.5 MW風(fēng)力機(jī)的氣動(dòng)性能?？紤]了葉尖損失、失速延遲效應(yīng)以及均勻來流和非穩(wěn)態(tài)來流兩種流動(dòng)風(fēng)場。研究正弦振蕩和剪切流場下由風(fēng)力機(jī)三維旋轉(zhuǎn)效應(yīng)產(chǎn)生的流動(dòng)分離和動(dòng)態(tài)失速效應(yīng)，分析了脈動(dòng)風(fēng)條件下葉片表面壓力頻譜特性，最后利用QBlade軟件建立優(yōu)化設(shè)計(jì)的風(fēng)力機(jī)模型，將計(jì)算結(jié)果與CFD模擬結(jié)果進(jìn)行對(duì)比，得出以下結(jié)論。

①標(biāo)準(zhǔn)狀況下，葉片表面基本為附著流動(dòng)，在葉根和葉尖會(huì)產(chǎn)生分離渦。隨著風(fēng)速的增加，在離心力和葉片展向壓力梯度的作用下，葉片吸力面開始發(fā)生流動(dòng)分離，深失速區(qū)域逐漸增大，風(fēng)速增大到21 m/s時(shí)，吸力面氣流完全分離。

②通過葉片展向上不同位置截面壓力功率譜可以看出，測點(diǎn)壓力在1倍風(fēng)輪轉(zhuǎn)頻時(shí)出現(xiàn)峰值且前緣峰值最高，由于葉尖渦的影響，葉尖吸力面后緣峰值最高；葉尖部位翼型截面吸力面對(duì)失速旋渦的敏感程度較高，壓力面敏感程度較低，葉根部位翼型受轉(zhuǎn)頻影響較小。