風力發電機葉片氣動性能數值模擬

王博，祁文軍，孫文磊，姜超

(1.新疆大學機械工程學院，新疆烏魯木齊830000;2.三峽新能源有限公司，新疆烏魯木齊830000)

風力發電機葉片是整個風力機中最為關鍵的部分，其氣動性能和結構性能的優劣直接決定風力機的工作效率和運行壽命。因此，風力機葉片的氣動設計和氣動性能計算是非常關鍵的工作。為了不斷提高設計和計算的準確性，國內外許多專家和學者在不斷研究和探索的基礎上提出了很多有效的理論和方法，主要包括葉素動量理論、風洞實驗、外場實驗及CFD數值模擬。其中CFD數值模擬方法是近年逐漸發展起來的理論和方法，在計算流體力學和計算機硬件的不斷發展推動下，逐步發展并不斷走向成熟，為認識風力機的流動性能提供了新的途徑。與前3種理論方法相比，有許多方面的優勢和更好的發展前景。CFD數值模擬技術可以方便快捷地計算出流場流動的各個定量結果，便于對各種流動參數進行分析，對設計方案進行評估，對后續的實驗驗證提供指導，避免了設計過程的盲目重復性，提高了風力機的整體設計水平，大幅度縮短了研制周期，降低了研制成本，代替了大量的實驗，成為風力機設計與分析的不可替代的工具［1-6］。

CFD是在計算機上實現一個特定的計算，就好像是在計算機上做一次全尺寸的風洞實驗。目前的應用主要集中在風電葉片翼型氣動特性數值模擬、風力機葉片氣動性能模擬及氣動彈性計算等方面［7］。然而，要將CFD直接應用于工程計算，并獲得可信賴的數值模擬計算結果，還存在計算域尺度與邊界條件的合理選取、能夠滿足計算精度要求的網格質量與網格數目的選擇以及合適的湍流模型選擇等問題。

作者應用UG軟件完成750 kW風力機葉輪和用于數值模擬的風力機流場的三維幾何建模，利用專用網格劃分軟件GAMBIT劃分網格，應用計算流體力學CFD中已被廣泛應用的Fluent軟件，對風力機建模葉片在額定風速和12個非額定風速工況下的氣動性能進行數值模擬，分析多個風速工況下的風力機建模葉片的氣動性能特點和差異，模擬仿真葉片氣動流場，計算葉輪的受力、扭轉力矩、輸出軸功率和風能利用效率等性能參數 ，觀察葉輪表面的壓強分布、流速分布、湍流強度、流速矢量等流態圖，驗證風力機氣動性能數值模擬的可靠性，并以現場750 kW葉片實際測得的功率數據曲線作為依據，驗證葉片建模的合理性及計算機仿真葉片的氣動性能是否滿足現場的要求。

1 數值模擬計算基礎——計算流體動力學方程

文中研究對象屬于湍流模型，通過分析模型屬于重整化k-ε模型 (RNG k-ε模型)。

RNG k-ε模型的輸運方程為:

其中:

式中:Gk表示平均速度梯度引起的湍動能的產生; Gb表示浮力影響引起的湍動能的產生;YM表示可壓湍流脈動膨脹對總的耗散率的影響。

FLUENT中，默認常數值為，C1ε=1.44，C2ε= 1.92，Cμ=0.09，湍動能k和耗散率ε的湍流普朗特數分別為σk=1.0，σε=1.3，αk和αε分別為湍動能k和耗散率ε的有效湍流普朗特數的倒數。

湍流黏性系數計算公式為:

在FLUENT中，用RNG k-ε模型時，默認設置為高雷諾數流動問題。如果是低雷諾數問題，必須進行相應的設置。

RNG模型在ε方程中加了一個條件，又考慮到了湍流漩渦，還提供了一個考慮低雷諾數流動黏性的解析公式，這些特點使得RNG k-ε模型比標準k-ε模型有更高的可信度和精度。

2 建模、網格劃分和邊界條件設定

2.1 幾何建模

風力發電機基本參數為:額定功率為750 kW，額定風速為15.0 m/s，啟動風速為4.0 m/s，切出風速為25.0 m/s，極大風速為70.0 m/s，功率調節方式為失速，風輪直徑為48.4 m，輪轂高度為50 m，葉片材料為玻璃纖維增強樹脂，葉片端線速度為57.4 m/s，風力機尖速比為3.83，葉輪額定轉速為22.5 r/min，葉輪直徑為48.4 m，掃風面積為1 840.0 m2，三葉片，選用NACA 4412翼型進行葉片建模。

風力機幾何建模包括:葉片建模，3葉片復制和旋轉體建模，輪轂建模，葉片周圍流場建模，風力機流場建模和整機流場建模，見圖1、圖2。

圖1 葉輪幾何建模

圖2 整機流場幾何建模

葉輪旋轉流場是包羅葉輪的扁圓盤，旋轉中心指向進風口方向4.0 m，指向出風口方向6.0 m，半徑30.0 m。整機流場是包羅整個葉輪旋轉流場的流場，整個風場由一個長方體和一個半圓體構成，長方體的長、寬、高分別為350、250、50 m，半圓體的半徑為125 m。旋轉中心到進風口 150 m，到出風口200 m。

網格劃分采用4面體網格。其中:旋轉流場節點209 568個，網格1 150 464個;靜止流場節點54 464個，網格282 124個;總共節點264 032個，網格1 432 588個，見圖3、圖4。

圖3 葉輪旋轉流場網格劃分

圖4 整機流場網格劃分

2.2 邊界條件設定

此次三維模型的氣動數值模擬計算采用分離式求解器 (Segregated)，三維穩態流動 (3D，Steady)，速度是絕對速度 (Absolute)，隱式求解 (Implicit)，RNG k-ε模型，流體材料設為空氣。動量采用二階迎風離散格式 (Second Order Up-wind)。壓力和速度耦合采用SIMPLE算法。

速度進口為15.0 m/s，同時計算變工況下風力機氣動數值模擬。風速變化為:4.0，6.0，8.0，10.0，12.0，14.0，15.0，16.0，18.0，20.0，22，24.0，25 m/s。自由出口。進出口之間壁面模擬為風洞壁面，靜止無滑移。葉輪旋轉流體采用MRF模型，額定轉速22.55 r/min，方向右手定則確定。

對于葉片壁面和近壁面區域采用動參考系(MRF)模型 (Multiple Reference Frame Model)。其中葉片壁面為旋轉壁面，相對旋轉流場旋轉速度為0，方向右手定則確定。對于近壁面區域，采用壁面函數法。在黏性流體流場中壁面處默認為非滑移邊界條件。

3 數值模擬計算結果與實測結果的對比

在額定工況和變工況下對風力機進行氣動數值模擬，計算收斂結束后，輸出各個工況下風力機的繞軸旋轉扭矩及軸力，并對記錄的旋轉數據和葉片已知數據進行整理計算，得到對應的轉矩功率，將其與現場測試的750 kW風力機實驗功率曲線對比，二者的功率曲線對比如圖5所示。

圖5 750 kW風力機計算機模擬與實測輸出功率曲線對比

從圖5中可以看出:

(1)利用FLUENT軟件對750 kW建模葉片在額定風速和啟動風速到切出風速的各個工況下的氣動性能進行模擬計算，因為模擬時不考慮風力機變槳距以及風力機內部的功率損耗等影響，可以說是在理想狀態下對風力機進行的數值模擬，風力機輸出功率與風速基本上成3次方曲線增長，這和理論上的風力機輸出功率計算公式相符。通過對風力機建模葉片進行氣動數值模擬計算，得出風力機氣動旋轉扭矩功率，在15 m/s額定風速情況下，模擬得到的轉矩功率為768.4 kW，與額定功率750 kW誤差只有2.5%，基本達到預期的數值，說明基于FLUENT的風力發電機氣動性能數值模擬計算和仿真是可靠的，也是可行的。

(2)在低于額定風速15 m/s時，風力發電機的實際輸出功率和建模葉片氣動性能數值模擬理想轉矩功率相差比較小，主要原因是在低于額定風速時風力機變槳距不對風機進行調節，而風速較小時對應的功率損耗也較小，這就使風力機實際工作狀態接近于數值模擬時的邊界條件，即理想狀態。

(3)在風速大于15 m/s時，實際輸出功率和模擬計算的理想輸出功率相差較遠，原因是風速大于額定風速時，風力機變槳距開始工作，會通過槳距調節，限制風輪的轉速，以保護風機，延長使用壽命，提高風力機工作穩定性，其結果必然降低風力機的功率和風能利用率;同時也因風速加大，變速箱及傳動軸間的摩擦等功率損耗會加大，這些都導致風速大于額定風速時實際輸出功率幾乎不增大，甚至減小。而數值模擬中并不會考慮這些因素，還按照理論上的功率輸出關系模擬計算其模擬輸出功率，按照3次方曲線方式增長。

(4)結合風力機的實際工作狀態，不可能讓風力機功率隨著風速不斷增大，因此數值模擬計算的理想輸出功率在小于額定輸出功率時會更接近于實際，而風力機絕大部分時間工作在小于或接近額定風速的工作狀態下，這樣，對于風力機氣動性能的研究在小于額定風速時會更有實際意義。而從圖3中可以看到，小于額定風速時模擬結果與實際結果非常接近，從而證明了建模葉片氣動性能較好地滿足了實際的需要，也證明了利用FLUENT軟件研究風力機的氣動性能振動特性是完全可行的。

4 風力機葉輪及周圍流場壓力、速度的分布特性

圖6—17分別為葉輪、流場出入口、葉輪周圍旋轉流場這3個區域的壓強與速度的分布圖與跡線圖，能直觀地、感性地認識和觀察葉輪及周圍流場的氣動特性和流動現象，以及壓力、流速等重要氣動參數的分布規律，并能很好地觀察葉輪以及風力機周圍的流場分布。

從圖6、圖7可以看到:葉輪正面為迎風面，在遇到葉輪阻礙之前，沒有發生流場的較大變化，因而氣動壓強分布較為均勻;而葉輪背面卻因為氣流通過旋轉葉輪時，紊流現象急劇加強，導致葉片表面流場發生了較大變化，故氣動壓強分布表現出明顯的不均勻性。

圖6 葉輪正面氣動壓強分布圖

圖7 葉輪背面氣動壓強分布圖

從圖8、圖9可以看到:葉輪正面和背面的速度分布幾乎沒有變化，但是由于風輪的轉動帶動了周圍氣流的旋轉，從而在接近葉輪周圍的圓柱形區域，空氣流場幾乎和葉輪同步旋轉。因而在選擇風力機氣動數值模擬邊界條件時，設定這一區域是和風輪同步旋轉的動流場，計算結果表明是正確合理的。

圖8 葉輪正面氣動速度分布圖

圖9 葉輪背面氣動速度分布圖

從圖10、圖11可以看到:在風力機葉輪周圍空氣受阻，氣動壓強大，遠離風力機區域空氣受阻力小，氣動壓強小，符合實際流場的空氣流動特性。

從圖12、圖13可以看到:進出口的空氣速度分布有較大的不同，這主要是因為空氣氣流在通過葉片風輪前紊流較弱，而通過風力機葉輪后，空氣氣流發生了明顯的紊流現象，流場中速度分布也隨之發生明顯變化。

圖10 流場入口氣動壓強分布圖

圖11 流場出口氣動壓強分布圖

圖12 流場入口氣動速度分布圖

圖13 流場出口的氣動速度分布圖

從圖14、圖15可以看到:風力機葉輪及周圍旋轉流場中壓強和速度的分布規律，葉輪處于旋轉中，故靠近葉輪旋轉區域速度大，遠離葉輪旋轉區域的速度減小;同樣，由于空氣氣流速度的影響，葉輪及周圍區域的氣流壓強也隨之發生變化，呈現和速度相呼應的分布特點。

圖14 葉輪周圍旋轉流場氣動壓強圖

圖15 葉輪周圍旋轉流場速度分布圖

從圖16、圖17中可以看到:通過FLUENT軟件 建立風力機葉輪旋轉流場，并進行數值模擬計算，得到的葉輪表面的壓強和速度跡線圖非常直觀、形象地描述了葉輪表面的壓強和速度流場的動態變化和分布，總的來看，按照旋轉方向呈流線型分布。

圖16 葉輪旋轉流場氣動壓強跡線圖

圖17 葉輪旋轉流場速度跡線圖

5 結論

(1)以NACA63XXX翼型葉片為基礎，通過UG軟件完成750 kW風力機葉輪和用于數值模擬的風力機流場的三維幾何建模，利用專用網格劃分軟件GAMBIT劃分網格，最后在FLUENT軟件中針對額定風速和12個非額定風速的工況進行氣動數值模擬計算，分析不同風速工況下的風力機氣動性能的特點和差異。通過風力機氣動數值模擬計算發現:風力機轉矩輸出功率與風速成3次方曲線增長，符合風力機輸出功率公式;得出的風力機氣動旋轉扭矩功率和750 kW風力機額定功率基本一致，得到的功率曲線在低于額定功率時也與現場實測750 kW葉輪得到的功率曲線相符，說明建模葉片的氣動性能基本符合750 kW機組的技術要求，基于FLUENT的風力發電機氣動性能數值模擬計算和仿真是可靠的、可行的。

(2)觀察葉輪表面以及整機周圍流場的壓強分布、流速分布、湍流強度、速度矢量等動態云圖，可以直觀和感性地認識和觀察葉輪和整機周圍流場的氣動特性和流動現象，以及其壓力、流速等氣動參數的分布特點，并且能很好地觀察葉輪以及整機周圍的流場分布。

(3)將模擬數據與實驗數據結合起來進行分析，有利于研究風力機葉輪及整機氣動性能和運行情況，數值模擬仿真研究在一定程度上可以替代實際風力機的運行實驗，在成本、效率方面有不可替代的優勢，也豐富了對風力機進行氣動分析的手段和方法，為風力機葉片及整機的設計、改型和研發工作提供技術參數和指導意見。

【1】SIMMS D，SCHREEK S，HAND M，et al.NREL Unsteady Aerodynamics Experiment in the NASA-Ames Wind Tunnel:A Comparison of Predictions to Measurements［R］.NREL/TP-500-29494，NREL，2001.

【2】BENJANIRAT S，SANKAR L N，XU G.Evaluation of turbulence models for the prediction of wind turbine aerodynamics［C］//41st AIAA Aerospace Sciences Meeting and Exhibit，Reno，American，2003.

【3】DUQUE E P N，BURKLUND M D，JOHNSON W.Navier-Stokes and comprehensive analysis performance predictions of the NREL Phase VI experiment［C］//41st ALAA Aerospace Sciences Meeting and Exhibit，Reno，America，2003.

【4】SEZER-UZOL Nilay，LONG Lyle N.3-D Time-Accurate CFD Simulations of Wind Turbine Rotor Flow Fields［C］//44th AIAA Aerospace Sciences Meeting and Exhibit，Reno，America，2006.

【5】ZAHLE Frederik，JPHANSEN Jeppe，SORENSEN Niels N.Wind TurbineAerodynamicsUsingan Incompressible Overset Grid Method［C］//Proceedings of the European Wind Energy Conference＆Exhibition，Milan，Italy，2007.

【6】van ROOIJ R P J O M，ARENS E A.Analysis of the Experimental and Computational Flow characteristics with Respect to the Augmented Lift Phenomenon Caused by Blade Rotation［C］//Journal of Physics:Conference Series，Lyngby，Denmark，2007.

【7】劉丹，李軍向，薛忠民，等.風力發電機葉片數值模擬綜述［J］.電氣技術，2010(7):7-13.