風力機翼型動態失速的數值模擬

宗濤

（湖北能源集團鄂州發電有限公司，湖北 鄂州 436032）

風力機翼型動態失速的數值模擬

宗濤

（湖北能源集團鄂州發電有限公司，湖北 鄂州 436032）

利用計算流體力學軟件FLUENT，對二維翼型NACA23012在非定常流動中的動態失速進行數值模擬研究，分析了翼型在振蕩周期內升力系數隨攻角的變化，并進一步研究了俯仰運動的幾個關鍵參數對翼型動態失速的影響，得出了翼型升力系數循環曲線在不同攻角、不同振幅、不同衰減頻率下的變化規律，并與實驗結果相比較。

翼型；動態失速；遲滯效應；數值模擬

0 引言

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

1 計算網格及邊界條件

對數值模擬流場進行網格劃分是求解控制方程的基礎，網格質量的好壞將直接影響數值計算的結果，而網格的生成關鍵在于網格點的合理分布。基于不同的網格生成方法，可使用結構化和非結構化網格來離散計算域。本文采用非結構化網格，并采用局部加密的方法，圖1所示為網格的局部放大圖。網格點在翼型前緣分布較密，整個流場計算域的網格數為69 380個。

圖1 網格分布圖Fig.1 Grid profile

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

2 結果分析

2.1 模擬值與實驗值比較

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

圖2 靜態時升力系數隨來流攻角的變化Fig.2 Lift coefficient with flow static angle of attack

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

2.2 遲滯效應的影響因素

從俯仰運動的數學表達式可知：簡諧俯仰運動的幾個關鍵參數為振幅Δα，振蕩角頻ω，平均攻角α0，以及俯仰運動的扭轉點位置。如果氣動計算條件不變，則俯仰運動的遲滯效應由這幾個參數決定[8]。

2.2.1 平均攻角的影響

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

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

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

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

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

2.2.2 振幅的影響

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

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

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

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

從圖7-圖9可見，當振幅為2°時，翼型完全失速；當振幅為6°時，動態失速圓環曲線再次出現，此時出現輕失速現象，失速由后緣分離引起；當振幅為10°時，遲滯效應更加明顯，發生深失速現象，上仰階段發生前緣分離，升力系數明顯增加。上仰過程模擬結果與實驗結果變化趨勢較吻合，下俯過程差異較大。

2.2.3 衰減頻率的影響

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

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

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

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

從圖11-圖13可見，當衰減頻率為0.1時，典型的動態失速圓環曲線即已出現。隨著衰減頻率的增加，循環曲線變窄，遲滯效應越來越不明顯，流動分離現象出現延遲，失速攻角增加，上仰過程模擬結果與實驗結果吻合較好。

3 結論

通過對以上動態失速現象的數值模擬，可以看出：

（1）翼型發生動態失速時，其升力系數峰值大于靜態時升力系數峰值，失速造成風力機實際輸出功率大于理論計算值。

（2）翼型俯仰振動的幾個關鍵參數對動態失速有很大的影響，平均攻角較小時，流體幾乎附著在翼型附近流動，遲滯環不明顯。

（3）振幅越大，動態失速越明顯，振幅較小時失速多為輕失速，由后緣分離引起；振幅較大時失速多為深失速，首先形成很大的前緣分離渦，該分離渦在翼型表面運動誘發出二次流動，引起升力系數的顯著變化。

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

模擬值較準確地計算了升力系數遲滯曲線的變化趨勢，基本趨勢與實驗結果較符合，但在具體值上與實測值還存在有一定的差異。

（References）

[1]尚麗娜，夏品奇.可變翼型動態失速氣動模型研究[J].中國科學:技術科學，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]劉雄，梁濕，陳嚴等.風力機翼型動態失速氣動特性仿真[J].工程力學,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]陳旭，郝輝，田杰等.水平軸風力機翼型動態失速特性的數值研究[J].太陽能學報,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]錢煒祺，符松，蔡金獅等.翼型動態失速的數值研究[J].空氣動力學學報，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]張正秋，鄒正平，劉火星等.振蕩翼型非定常流動數值模擬研究[J].燃氣渦輪試驗與研究，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