風(fēng)力機(jī)翼型動(dòng)態(tài)失速的數(shù)值模擬

宗濤

（湖北能源集團(tuán)鄂州發(fā)電有限公司，湖北 鄂州 436032）

風(fēng)力機(jī)翼型動(dòng)態(tài)失速的數(shù)值模擬

宗濤

（湖北能源集團(tuán)鄂州發(fā)電有限公司，湖北 鄂州 436032）

利用計(jì)算流體力學(xué)軟件FLUENT，對(duì)二維翼型NACA23012在非定常流動(dòng)中的動(dòng)態(tài)失速進(jìn)行數(shù)值模擬研究，分析了翼型在振蕩周期內(nèi)升力系數(shù)隨攻角的變化，并進(jìn)一步研究了俯仰運(yùn)動(dòng)的幾個(gè)關(guān)鍵參數(shù)對(duì)翼型動(dòng)態(tài)失速的影響，得出了翼型升力系數(shù)循環(huán)曲線在不同攻角、不同振幅、不同衰減頻率下的變化規(guī)律，并與實(shí)驗(yàn)結(jié)果相比較。

翼型；動(dòng)態(tài)失速；遲滯效應(yīng)；數(shù)值模擬

0 引言

動(dòng)態(tài)失速是指一個(gè)振蕩（或做其他非定常運(yùn)動(dòng)）的壓力面在超過其臨界迎角時(shí)繞流流場(chǎng)發(fā)生非定常分離和失速的現(xiàn)象。對(duì)于風(fēng)力機(jī)翼型來說，在大攻角、低雷諾數(shù)來流的工況下，當(dāng)翼型產(chǎn)生振顫時(shí)極易出現(xiàn)動(dòng)態(tài)失速現(xiàn)象。當(dāng)風(fēng)輪處于失速狀態(tài)運(yùn)行時(shí)，最大輸出功率和最大葉片載荷將會(huì)同時(shí)出現(xiàn)，而目前不計(jì)失速影響的理論計(jì)算載荷只有實(shí)測(cè)值的50%-70%。失速風(fēng)輪在高風(fēng)速下所產(chǎn)生的這種很大輸出功率可能引起發(fā)電機(jī)的損壞，因此，在設(shè)計(jì)時(shí)必須認(rèn)真考慮動(dòng)態(tài)失速的作用。近年來，許多學(xué)者針對(duì)這一問題進(jìn)行了大量的研究[1-3]，如劉雄[4]等人通過實(shí)驗(yàn)初步分析了NACA4415翼型的動(dòng)態(tài)氣動(dòng)特性，陳旭[5]等人對(duì)NREL s809翼型的二維動(dòng)態(tài)流場(chǎng)進(jìn)行了數(shù)值模擬等。本文應(yīng)用CFD軟件FLUENT對(duì)二維翼型NACA23012的動(dòng)態(tài)失速現(xiàn)象進(jìn)行數(shù)值模擬。

1 計(jì)算網(wǎng)格及邊界條件

對(duì)數(shù)值模擬流場(chǎng)進(jìn)行網(wǎng)格劃分是求解控制方程的基礎(chǔ)，網(wǎng)格質(zhì)量的好壞將直接影響數(shù)值計(jì)算的結(jié)果，而網(wǎng)格的生成關(guān)鍵在于網(wǎng)格點(diǎn)的合理分布?；诓煌木W(wǎng)格生成方法，可使用結(jié)構(gòu)化和非結(jié)構(gòu)化網(wǎng)格來離散計(jì)算域。本文采用非結(jié)構(gòu)化網(wǎng)格，并采用局部加密的方法，圖1所示為網(wǎng)格的局部放大圖。網(wǎng)格點(diǎn)在翼型前緣分布較密，整個(gè)流場(chǎng)計(jì)算域的網(wǎng)格數(shù)為69 380個(gè)。

圖1 網(wǎng)格分布圖Fig.1 Grid profile

動(dòng)態(tài)失速中一個(gè)關(guān)鍵的問題是翼型在流場(chǎng)中發(fā)生了振蕩運(yùn)動(dòng)，因此，翼型振蕩運(yùn)動(dòng)的形式對(duì)流場(chǎng)的變化起著關(guān)鍵性作用。研究表明，風(fēng)力機(jī)葉片發(fā)生的振顫可認(rèn)為是按照正弦規(guī)律變化的運(yùn)動(dòng)。當(dāng)葉片發(fā)生振蕩時(shí)，其截面翼型與流場(chǎng)中來流的相對(duì)方位是隨時(shí)間不斷變化的，這個(gè)變化所引起的主要影響是來流對(duì)翼型的攻角也產(chǎn)生周期性的振蕩變化。由流場(chǎng)運(yùn)動(dòng)的相對(duì)性，可以設(shè)翼型振蕩運(yùn)動(dòng)過程中其攻角的變化形式為α=α0+Δα?sin(ω?t)，其中α0為平均攻角，Δα為振幅。取翼型的振蕩軸心為弦長的四分之一處。邊界條件設(shè)為無窮遠(yuǎn)壓力邊界，空氣視為理想氣體，雷諾數(shù)Re=1.5×106，計(jì)算采用k-ω SST模式。

2 結(jié)果分析

2.1 模擬值與實(shí)驗(yàn)值比較

文獻(xiàn)[6]給出了翼型NACA23012靜態(tài)時(shí)、繞振蕩軸心規(guī)律振蕩時(shí)升力系數(shù)（翼型所受到的升力與氣流動(dòng)壓和參考面積的乘積之比）隨攻角變化的實(shí)驗(yàn)值。圖2所示為靜態(tài)時(shí)升力系數(shù)隨攻角的變化。由圖知，升力系數(shù)隨攻角的增加，先增加后減小，失速攻角在15°左右，升力系數(shù)峰值約為1.3。在實(shí)驗(yàn)條件與數(shù)值模擬條件一致的情況下，模擬值與實(shí)驗(yàn)測(cè)得的結(jié)果吻合。

圖2 靜態(tài)時(shí)升力系數(shù)隨來流攻角的變化Fig.2 Lift coefficient with flow static angle of attack

圖3 振蕩翼型升力系數(shù)隨來流攻角的變化Fig.3 Variation of lift coefficient of oscillating airfoil with inflow angle of attack

圖4 不同攻角下翼型周圍速度矢量圖Fig.3 Velocity vector around an airfoil at different angles of attack

圖3所示為翼型俯仰振蕩時(shí)升力系數(shù)隨攻角的變化，圖4為不同攻角下翼型周圍速度矢量圖。采用的數(shù)值計(jì)算條件為：自由來流Ma=0.12，衰減頻率，其中c為弦長，U∞為遠(yuǎn)場(chǎng)來流速度；平均攻角α0=15°，振幅Δα=10°。由圖易知，振蕩上仰階段升力系數(shù)隨著攻角的增加而增加，當(dāng)上仰至最大值時(shí)，翼型向下振蕩，此時(shí)升力系數(shù)急速下降，而后有小幅上升。1點(diǎn)為失速延遲點(diǎn)，由于翼型的振蕩作用延遲了失速現(xiàn)象的發(fā)生，甚至不發(fā)生失速，使得升力系數(shù)峰值遠(yuǎn)大于靜態(tài)時(shí)峰值，幾乎為靜態(tài)時(shí)峰值的2倍；當(dāng)翼型上仰至2點(diǎn)時(shí)，發(fā)生前緣分離，使升力系數(shù)明顯增加。隨著攻角的進(jìn)一步增加，主分離渦的位置升高，向后緣移動(dòng)，并在前緣和后緣誘導(dǎo)出二次分離渦；進(jìn)入下俯過程以后，后緣的分離渦變大，并和主分離渦一起脫離翼型表面，此時(shí)翼型上方的渦為前緣誘導(dǎo)出的二次渦，當(dāng)下俯至4點(diǎn)時(shí)該二次渦到達(dá)翼型尾緣，并隨著翼型的下俯振蕩逐漸脫離翼型表面；當(dāng)下俯至5點(diǎn)時(shí)，翼型表面的流動(dòng)又恢復(fù)為附著流動(dòng)。模擬值與實(shí)驗(yàn)所得結(jié)果變化趨勢(shì)一致，但具體值存在有差異，遲滯環(huán)小于實(shí)驗(yàn)所得結(jié)果，數(shù)值精度的不足和湍流模式的缺陷是造成這一差異的主要原因，尤其是湍流模式方面[7]。同時(shí)上述數(shù)值模擬所選擇的雷諾數(shù)Re、衰減頻率k的取值與實(shí)際工況相符。

2.2 遲滯效應(yīng)的影響因素

從俯仰運(yùn)動(dòng)的數(shù)學(xué)表達(dá)式可知：簡(jiǎn)諧俯仰運(yùn)動(dòng)的幾個(gè)關(guān)鍵參數(shù)為振幅Δα，振蕩角頻ω，平均攻角α0，以及俯仰運(yùn)動(dòng)的扭轉(zhuǎn)點(diǎn)位置。如果氣動(dòng)計(jì)算條件不變，則俯仰運(yùn)動(dòng)的遲滯效應(yīng)由這幾個(gè)參數(shù)決定[8]。

2.2.1 平均攻角的影響

數(shù)值計(jì)算條件為：自由來流Ma=0.12，衰減頻率k=0.10，振幅Δα=10°，取平均攻角分別為6°、12°、20°，如圖5-圖7所示。

圖5 平均攻角為6°時(shí)振蕩翼型升力系數(shù)隨攻角的變化Fig.5 The lift coefficient of oscillating airfoil varies with angle of attack when the mean angle of attack is 6°

圖6 平均攻角為12°時(shí)振蕩翼型升力系數(shù)隨攻角的變化Fig.6 The lift coefficient of oscillating airfoil varies with angle of attack when the mean angle of attack is 12°

圖7 平均攻角為20°時(shí)振蕩翼型升力系數(shù)隨攻角的變化Fig.7 The lift coefficient of oscillating airfoil varies with angle of attack when the mean angle of attack is 20°

從圖5-圖7可見，當(dāng)平均攻角最小時(shí)，流體幾乎附著在翼型附近流動(dòng)，升力系數(shù)曲線為特殊的橢圓形。當(dāng)平均攻角為12°時(shí)，出現(xiàn)動(dòng)態(tài)失速，由實(shí)驗(yàn)值知直到翼型下俯至4°時(shí)升力系數(shù)又恢復(fù)為附著流動(dòng)時(shí)的值，這表明了下俯振蕩延遲了流動(dòng)恢復(fù)為附著流動(dòng)的過程。當(dāng)平均攻角為20°時(shí)，最大升力系數(shù)幾乎是靜態(tài)時(shí)升力系數(shù)峰值的2倍，在上仰階段，升力曲線斜率呈現(xiàn)比較大的非線性變化。失速攻角在25°左右，升力系數(shù)達(dá)到峰值后，先下降后又升高，再次出現(xiàn)峰值。二次峰值的出現(xiàn)和二次渦的形成有關(guān)。模擬所得結(jié)果與實(shí)驗(yàn)結(jié)果變化趨勢(shì)一致，具體值存在有一定差異，模擬所得遲滯效應(yīng)沒有實(shí)驗(yàn)結(jié)果明顯。

2.2.2 振幅的影響

數(shù)值計(jì)算條件為：自由來流Ma=0.12，衰減頻率k=0.10，平均攻角α0=18°，取振幅分別為2°、6°、10°，如圖8-圖10所示。

圖8 振幅為2°時(shí)振蕩翼型升力系數(shù)隨攻角的變化Fig.8 The lift coefficient of oscillating airfoil varies with angle of attack when the amplitude is 2°

圖9 振幅為6°時(shí)振蕩翼型升力系數(shù)隨攻角的變化Fig.9 The lift coefficient of oscillating airfoil varies with angle of attack when the amplitude is 6°

圖10 振幅為10°時(shí)振蕩翼型升力系數(shù)隨攻角的變化Fig.10 The lift coefficient of oscillating airfoil varies with angle of attack when the amplitude is 10°

從圖7-圖9可見，當(dāng)振幅為2°時(shí)，翼型完全失速；當(dāng)振幅為6°時(shí)，動(dòng)態(tài)失速圓環(huán)曲線再次出現(xiàn)，此時(shí)出現(xiàn)輕失速現(xiàn)象，失速由后緣分離引起；當(dāng)振幅為10°時(shí)，遲滯效應(yīng)更加明顯，發(fā)生深失速現(xiàn)象，上仰階段發(fā)生前緣分離，升力系數(shù)明顯增加。上仰過程模擬結(jié)果與實(shí)驗(yàn)結(jié)果變化趨勢(shì)較吻合，下俯過程差異較大。

2.2.3 衰減頻率的影響

數(shù)值計(jì)算條件為：自由來流Ma=0.12，平均攻角α0=10°，振幅Δα=10°，取衰減頻率分別為0.10、0.16、0.20，如圖11-圖13所示。

圖11 衰減頻率為0.10時(shí)翼型升力系數(shù)隨攻角的變化Fig.11 Variation of lift coefficient of airfoil with angle of attack when attenuation frequency is 0.10

圖12 衰減頻率為0.16時(shí)翼型升力系數(shù)隨攻角的變化Fig.12 Variation of lift coefficient of airfoil with angle of attack when attenuation frequency is 0.16

圖13 衰減頻率為0.20時(shí)翼型升力系數(shù)隨攻角的變化Fig.13 Variation of lift coefficient of airfoil with angle of attack when attenuation frequency is 0.20

從圖11-圖13可見，當(dāng)衰減頻率為0.1時(shí)，典型的動(dòng)態(tài)失速圓環(huán)曲線即已出現(xiàn)。隨著衰減頻率的增加，循環(huán)曲線變窄，遲滯效應(yīng)越來越不明顯，流動(dòng)分離現(xiàn)象出現(xiàn)延遲，失速攻角增加，上仰過程模擬結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好。

3 結(jié)論

通過對(duì)以上動(dòng)態(tài)失速現(xiàn)象的數(shù)值模擬，可以看出：

（1）翼型發(fā)生動(dòng)態(tài)失速時(shí)，其升力系數(shù)峰值大于靜態(tài)時(shí)升力系數(shù)峰值，失速造成風(fēng)力機(jī)實(shí)際輸出功率大于理論計(jì)算值。

（2）翼型俯仰振動(dòng)的幾個(gè)關(guān)鍵參數(shù)對(duì)動(dòng)態(tài)失速有很大的影響，平均攻角較小時(shí)，流體幾乎附著在翼型附近流動(dòng)，遲滯環(huán)不明顯。

（3）振幅越大，動(dòng)態(tài)失速越明顯，振幅較小時(shí)失速多為輕失速，由后緣分離引起；振幅較大時(shí)失速多為深失速，首先形成很大的前緣分離渦，該分離渦在翼型表面運(yùn)動(dòng)誘發(fā)出二次流動(dòng)，引起升力系數(shù)的顯著變化。

（4）衰減頻率的大小決定了翼型振蕩角頻的大小，衰減頻率較小時(shí)失速延遲現(xiàn)象明顯。

模擬值較準(zhǔn)確地計(jì)算了升力系數(shù)遲滯曲線的變化趨勢(shì)，基本趨勢(shì)與實(shí)驗(yàn)結(jié)果較符合，但在具體值上與實(shí)測(cè)值還存在有一定的差異。

（References）

[1]尚麗娜，夏品奇.可變翼型動(dòng)態(tài)失速氣動(dòng)模型研究[J].中國科學(xué):技術(shù)科學(xué)，2015,(09):975-983.Shang Lina，Xia Pingqi.Study of aerodynamic model for morphing airfoils in dynamic stall[J].Scientia Sinica(Technologica),2015,(09):975-983.

[2]B.Stoevesandt,Joachim Peinke,A.Shishkin,Claus Wag?ner.Numerical Simulation of Dynamic Stall using Spectral/hp Method[M].Wind Energy,2004.

[3]O.Frederich,U.Bunge,C.Mockett,F.Thiele.FlowPre?diction Around an Oscillating NACA0012 Airfoil at Re= 1,000,000[M],IUTAM Bookseries,IUTAM Symposium on‘Unsteady Separated Flows and their Control’,Corfu,Greece,18-22 June 2007.

[4]劉雄，梁濕，陳嚴(yán)等.風(fēng)力機(jī)翼型動(dòng)態(tài)失速氣動(dòng)特性仿真[J].工程力學(xué),2015,(03):203-211.Liu Xiong，Liang Xiang，Chen Yan,et al.Dynamic stallsimulation ofwind turbine airfoils[J].Engi?neering Mechanics,2015,(03):203-211.

[5]陳旭，郝輝，田杰等.水平軸風(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.

[6]J.G.Leishman.Dynamic stall experiment on the NACA23012 aerofoil[J].Experiments in Fluids,1990, 9(1-2):49-58.

[7]錢煒祺，符松，蔡金獅等.翼型動(dòng)態(tài)失速的數(shù)值研究[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2001,19(4):427-433.Qian Weiqi,Fu song,Cai Jinshi,et al.Numerical study of airfoil dynamic stall[J].Acta Aerodynamica Sinica,2001,19(4):427-433.

[8]張正秋，鄒正平，劉火星等.振蕩翼型非定常流動(dòng)數(shù)值模擬研究[J].燃?xì)鉁u輪試驗(yàn)與研究，2009,22(3):1-8.Zhang Zhengqiu,Zhou Zhengping,Liu Huoxing,et al.Numerical Study of Two Dimensional Flow on an Oscillating Airfoil[J].Gas Turbine Experiment and Research,2009,22(3):1-8.

Numerical Simulation of Airfoil Dynamic Stall of Wind Turbine

ZONG Tao
（Hubei Energy Group Ezhou Power Generation Company Limited，Ezhou Hubei436032，China）

In this paper,the dynamic stall of the 2-d airfoil NACA23012 is simulated by using the computational fluid mechanics software FLUENT,the change of lift coefficient with power An?gle in the oscillating cycle is analyzed,and the influence of several key parameters of pitch motion on the dynamic stall of airfoil is studied.The variation rules of the cyclic curve of the wing lift co?efficient in different Angle,amplitude and attenuation frequency are obtained and compared with the experimental results.

airfoil;dynamic stall;hysteresis effect;numerical simulation

TK83

B

1006-3986(2017)04-0050-05

2017-02-21

宗 濤(1986),男,湖北武漢人，碩士，工程師。

10.19308/j.hep.2017.04.011