屋頂H型垂直軸風力機氣動性能CFD計算方法研究

■ 梅毅 曲建俊

（1.中國電力工程顧問集團華北電力設計院有限公司；2.哈爾濱工業(yè)大學機電工程學院）

0 引言

風能是一種清潔的可再生能源，風力發(fā)電提供的大量清潔電力可替代常規(guī)火力發(fā)電，具有巨大的生態(tài)環(huán)境效益和社會效益。我國風電產(chǎn)業(yè)自2006年以來保持強勁的發(fā)展勢頭，裝機容量屢創(chuàng)新高。然而，在風電裝機高速增長的同時，棄風限電問題也日益嚴重[1]?！讹L電發(fā)展“十三五”規(guī)劃》中明確指出，在負荷中心發(fā)展分布式風力發(fā)電，是解決風電消納問題的重要途徑之一。容量相對較小、分散布置的分布式風電，已成為電力系統(tǒng)的重要組成部分[2,3]。近年來，城市樓頂風能的利用引起了國內(nèi)外學者的廣泛關注[4-6]。與應用在野外風場的水平軸余量不同，垂直軸風力機具有無需對風裝置、結構相對簡單、安裝維修方便等優(yōu)勢，尤其是H型垂直軸風力機，結構簡單緊湊、占地面積小、噪音小，更適合作為城市屋頂?shù)娘L能利用設備。因此，各國研究人員對其進行了針對性研究[7-9]。

風輪氣動性能對風電機組的運行性能有重要影響。計算流體動力學(CFD)方法是研究垂直軸風力機氣動性能的重要工具。Qin[10]和Untaroiu[11]比較了不同時間步長應用于H型垂直軸風力機CFD計算的結果，分別認為選用0.0002 s和0.001 s的時間步長合適，但研究結論缺乏通用性。Wang[12]和Edwards[13]分別用幾種常用湍流模型計算了俯仰振蕩翼型的升阻力系數(shù)，指出湍流模型SSTk-ω模型更適用于H型垂直軸風力機數(shù)值計算，但單個翼型俯仰振蕩并不能反映旋轉(zhuǎn)中的垂直軸風力機的工作狀態(tài)。國家標準規(guī)定功率系數(shù)是反映風力機氣動性能的重要參數(shù)之一[14]，為保證H型垂直軸風力機氣動性能數(shù)值計算結果可信，有必要考察CFD方法用于計算風輪功率系數(shù)的精度。本文以一種屋頂H型垂直軸風力機為研究對象，分析不同湍流模型、網(wǎng)格單元形狀和時間步長對功率系數(shù)計算適用性的影響，以期為工程設計提供參考。

1 實驗風力機

文獻[15,16]所述的為一臺用于城市屋頂風力發(fā)電的H型垂直軸風力機(如圖1所示)，在風洞中進行了全尺寸實驗，與其他垂直軸風力機實驗相比，該實驗修正了由風輪支撐臂阻力及測試系統(tǒng)誤差等因素造成的功率損失，獲得了理想條件下的風力機功率系數(shù)值，實驗功率系數(shù)值被各國學者用作H型垂直軸風力機數(shù)值模擬技術驗證研究[17,18]。風力機風輪直徑為2.5 m，旋轉(zhuǎn)主軸直徑為0.1 m，葉片為NACA0015翼型，當葉片長3 m，弦長0.4 m。當風速為10 m/s時，風洞實驗測得尖速比λ=1.6時有最大功率系數(shù)0.34。該尖速比下，葉片的雷諾數(shù)在1×105～7.5×105之間變化，屬于典型的垂直軸風力機工作的雷諾數(shù)范圍[19]。本文采用不同湍流模型、網(wǎng)格單元形狀和時間步長計算該實驗風力機在各尖速比下的功率系數(shù)。鑒于現(xiàn)代垂直軸風力機常采用功率控制方法使之在最大功率系數(shù)下工作[20]，本文還將通過比較風輪在最大功率系數(shù)時葉片的瞬時力矩來分析不同方法計算時差異產(chǎn)生的原因。

圖1 實驗風力機

2 評價指標

風力機流場是非定常流場，CFD計算獲得相關氣動性能數(shù)據(jù)前，須使流場充分發(fā)展。CFD方法是否適用要考慮計算效率和計算精度兩個方面，即考察風輪旋轉(zhuǎn)一周平均力矩達到收斂所需運算時間和風輪功率系數(shù)的計算精度。為使流場充分發(fā)展，規(guī)定風輪旋轉(zhuǎn)一周平均力矩收斂的標準為相鄰兩周風輪平均力矩變化百分比r＜1%，如式(1)所示：

式中，Tave(i)為風輪旋轉(zhuǎn)第i周(i≥1)時的平均力矩。

設待考察的尖速比為n個，第j(j=1，2，...，n)個尖速比下功率系數(shù)CP的計算精度用εj來表示，如式(2)所示：

式中，CPej為第j個尖速比下風輪功率系數(shù)實驗值；CPsj為第j個尖速比下風輪功率系數(shù)CFD計算值。CP可根據(jù)文獻[21]中所述方法計算得到。

n個尖速比下CP的計算精度為εave，可根據(jù)式(3)計算：

3 數(shù)值計算方法

式(4)為由連續(xù)性方程和動量方程組成的控制方程組。

式中，ρ為空氣密度，kg/m3；p為壓強，Pa；μe為有效粘性系數(shù)；ui、uj分別為各坐標方向上的速度分量，m/s；t為時間，s。

由風輪二維模型組成的計算域ABCD如圖2所示。設風輪直徑為Φ，AD=BC=10Φ，AB=CD=20Φ，風輪回轉(zhuǎn)中心距AD為5Φ，距BC為15Φ，風輪旋轉(zhuǎn)域直徑為2Φ。

圖2 計算域示意圖

用Fluent求解控制方程組， AD為速度進口邊界，來流速度為10 m/s，壓力為大氣壓力，BC為壓力出口邊界。葉片、轉(zhuǎn)軸及計算域邊界AB和CD均為無滑移壁面，旋轉(zhuǎn)域和外部靜止域結合面為interface邊界條件。用壓力基求解方法，速度與壓力耦合采用PISO算法，壓力插值項為PRESTO格式，動量、湍動能和耗散率均采用QUICK格式，各項殘差控制在10-5。

4 結果與分析

4.1 不同湍流模型比較

本文主要比較RNGk-ε模型、Realizablek-ε模型和SSTk-ω模型3種常用的兩方程模型。計算域網(wǎng)格劃分采用四邊形單元，設置時間步長為風輪旋轉(zhuǎn)0.5°所需時間。圖3為尖速比λ=1.6時3種湍流模型計算的風輪平均力矩變化百分比。由圖3可知，3種模型計算的風輪流場均在旋轉(zhuǎn)至第6周時達到充分發(fā)展(r≤1%)。

表1為3種湍流模型在不同尖速比下的平均計算時間。由表1可知，RNGk-ε模型和Realizablek-ε模型計算的流場達到充分發(fā)展的平均耗時分別為8.46 h和8.28 h，而SSTk-ω模型為8.94 h。SSTk-ω模型計算耗時分別比RNGk-ε和Realizablek-ε模型多5.67%和7.97%。

圖3 λ=1.6時3種湍流模型計算的風輪平均力矩變化百分比

表1 3種湍流模型在不同尖速比下的平均計算時間

圖4為3種湍流模型計算得到的風輪功率系數(shù)曲線，各尖速比下3種湍流模型的計算精度如表2所示。由表2可知，SSTk-ω模型的模擬精度最高。

圖4 3種湍流模型計算的風輪功率系數(shù)曲線

表2 3種湍流模型的計算精度

圖5為尖速比λ=1.6時，3種湍流模型計算的葉片瞬時力矩曲線圖。從圖5中可以看到，葉片瞬時力矩Tbinst在風輪上游區(qū)域達到峰值后逐漸降低，進入下游區(qū)域后維持在較低水平，原因是風輪吸收的大部分空氣能量主要來自于風輪上游，下游葉片獲取空氣來流的能量較少。圖5中3條瞬時力矩曲線相差較大處在60°＜θ＜150°區(qū)間內(nèi)，SSTk-ω模型在吸收空氣能量較大的上游區(qū)域模擬值大于RNGk-ε模型和Realizablek-ε模型，導致其功率系數(shù)計算結果高于后兩者。

圖5 3種湍流模型計算的葉片瞬時力矩曲線

4.2 不同網(wǎng)格單元形狀比較

如圖6所示，采用兩種網(wǎng)格單元劃分計算域，葉片周圍用O型網(wǎng)格加密，葉片壁面網(wǎng)格數(shù)均為984，網(wǎng)格增長率均為1.07。用SSTk-ω模型封閉控制方程，設置非定常計算時間步長為風輪旋轉(zhuǎn)0.5°所需時間。通過網(wǎng)格無關性檢驗確定計算域網(wǎng)格單元數(shù)量如表3所示。

圖6 四邊形和三角形網(wǎng)格劃分計算域

表3 計算域網(wǎng)格數(shù)量對比

圖7為尖速比λ=1.6時兩種形狀網(wǎng)格計算的風輪平均力矩變化百分比。由圖7可知，采用四邊形網(wǎng)格計算時，風輪旋轉(zhuǎn)至第6周流場達到充分發(fā)展；而用三角形網(wǎng)格計算時，風輪旋轉(zhuǎn)至第11周才滿足收斂標準。

圖7 λ=1.6時兩種形狀網(wǎng)格計算的風輪平均力矩變化百分比

表4為兩種形狀網(wǎng)格在不同尖速比下的平均計算時間。表4中數(shù)據(jù)表明，三角形網(wǎng)格模擬風輪旋轉(zhuǎn)一周平均需0.87 h，耗時僅為四邊形網(wǎng)格的58%；但四邊形網(wǎng)格計算的收斂速度比三角形網(wǎng)格快5周，總耗時約為三角形網(wǎng)格的93%，計算效率相對更高。

表4 兩種形狀網(wǎng)格在不同尖速比下的平均計算時間

兩種形狀網(wǎng)格模擬所得風輪功率曲線如圖8所示。

表5為兩種形狀網(wǎng)格的計算精度。從表5中可以看到，四邊形網(wǎng)格的計算精度為7.40%，而三角形網(wǎng)格為7.89%，兩者較為接近。

圖8 兩種形狀網(wǎng)格計算的風輪功率系數(shù)曲線

表5 兩種形狀網(wǎng)格的計算精度

圖9為尖速比λ=1.6時，兩種形狀網(wǎng)格最大功率系數(shù)點時計算所得的葉片瞬時力矩曲線，其中相差較大之處主要位于80°＜θ＜120°區(qū)間內(nèi)，而在其他位置兩條曲線相對吻合。在網(wǎng)格尺寸和數(shù)量相同的條件下，四邊形網(wǎng)格的節(jié)點數(shù)多于三角形網(wǎng)格，因此其模擬風輪旋轉(zhuǎn)一周所需時間比三角形網(wǎng)格長。有限體積法將物理量存儲于網(wǎng)格單元中心點上，沿網(wǎng)格邊界線計算流場參數(shù)。垂直軸風力機風輪旋轉(zhuǎn)時，流場會沿葉片壁面彎曲。四邊形網(wǎng)格單元由于有4條邊界線，在某些位置流體流動的矢量方向會沿著其中兩條網(wǎng)格邊界線；而三角形網(wǎng)格只有3條邊界線，流體流動的矢量方向最多與其中一條網(wǎng)格邊界線一致，因而采用四邊形網(wǎng)格計算時產(chǎn)生的離散誤差小于三角形網(wǎng)格，求解穩(wěn)定性優(yōu)于三角形網(wǎng)格，計算時流場達到穩(wěn)定耗時更少。

圖9 兩種形狀網(wǎng)格計算的葉片瞬時力矩曲線

4.3 不同時間步長比較

選取時間步長為風輪旋轉(zhuǎn) 0.25°、0.5°、1°和2°所需時間用于CFD計算，記為0.25°ω-1、0.5°ω-1、1°ω-1和 2°ω-1。選用 SSTk-ω模型，計算域用四邊形網(wǎng)格劃分。圖10為尖速比λ=1.6時4種時間步長計算的風輪平均力矩變化百分比，表6為4種時間步長在不同尖速比下的平均計算時間。

從圖10和表6可以看出，較大的時間步長有利于提高Tave的收斂速度，減小時間步長導致計算時間增加，時間步長由2°ω-1縮短為1°ω-1、0.5°ω-1和 0.25°ω-1后，Tave平均收斂時間分別增至原收斂時間的1.55倍、2.98倍和4.66倍。

圖10 λ=1.6時4種時間步長計算的風輪平均力矩變化百分比

表6 4種時間步長在不同尖速比下的平均計算時間

由圖11和表7可知，減小時間步長有助于得到較準確的計算結果，但時間步長減小到一定程度后計算精度不會再大幅提高，且會使計算耗時增加。

圖11 4種時間步長計算的風輪功率系數(shù)曲線

表7 4種時間步長的計算精度

從圖12中可以看出，尖速比λ=1.6時，4種時間步長計算的葉片瞬時力矩差異主要發(fā)生在吸收風能較多的風輪上游區(qū)域，時間步長為2°ω-1和1°ω-1時計算的葉片瞬時力矩大于0.25°ω-1和0.5°ω-1時的結果，導致時間步長為2°ω-1和1°ω-1時計算的功率系數(shù)更大。

圖12 4種時間步長計算的葉片瞬時力矩曲線

5 結論

改變湍流模型、網(wǎng)格單元類型和時間步長對屋頂H型垂直軸風力機功率系數(shù)的計算有較大影響，與實驗數(shù)據(jù)對比驗證的結果表明：

1)與 RNGk-ε和 Realizablek-ε模型相比，SSTk-ω模型計算耗時略長，但功率系數(shù)計算更接近實驗值；

2)四邊形網(wǎng)格的功率系數(shù)的總體計算精度與三角形網(wǎng)格相近，但計算效率更高；

3)采用較大的時間步長會降低風輪功率系數(shù)的計算精度，較小的時間步長可提高模擬精度，但會增加計算時間，將風輪旋轉(zhuǎn)0.5o所需時間作為非定常計算的時間步長比較合適。

[1]熊敏鵬,張嚴,袁家海,等.我國風電的經(jīng)濟性評價及政策建議[J].中國能源,2016,38(10):20－26.

[2]張旭,呂新良,宋曉林,等.小型分布式風力發(fā)電系統(tǒng)設計及控制技術[J].陜西電力,2012,(1):75－78.

[3]張敏吉,梁嘉,孫洋洲,等.分布式風電——電池儲能系統(tǒng)可用性分析[J].電力建設,2016,37(11):29－34.

[4]袁行飛,余俊偉.屋頂安裝型風力機塔架風振反應分析[J].浙江大學學報(工學版),2013,17(11):1911－1916.

[5]潘雷.建筑環(huán)境中的風能利用[D].濟南:山東建筑大學,2006.

[6]楊蓉,彭興黔.高層建筑屋頂風能利用的數(shù)值模擬[J].華僑大學學報(自然科學版),2012,33(1):69－73.

[7]Battisti L,Zanne L,Anna S D,et al.Aerodynamic measurements on a vertical axis wind turbine in a large scale wind tunnel[J].Journal of Energy Resources Technology,2011,133:0312011－0312019.

[8]李巖,田川公太郎.葉片附著物對直線翼垂直軸風力機性能影響的風洞試驗[J].可再生能源,2008,26(5):21－23.

[9]Keiana M.A modi fi ed streamtube model for vertical axis wind turbines[J].Wind Engineering,2012,36(2):145－180.

[10]Qin N,Howell R,Durrani N,et al.Unsteady fl ow simulation and dynamic stall behaviour of vertical axis wind turbine blades[J].Wind Engineering,2011,35(4):511－527.

[11]Untaroju A,Wood H G,Allaire P E,et al.Investigation of self-starting capability of vertical axis wind turbines using a computational fl uid dynamics approach[J].Journal of Solar Energy Engineering,2011,133:0410101－0410108.

[12]Wang S Y,Ingham D B,Ma L,et al.Numerical investigations on dynamic stall of low Reynolds number fl ow around oscillating airfoils[J].Computers and Fluids,2010,39(9):1529－1541.

[13]Edwards J M,Danao A L,Howell R J.Novel experimental power curve determination and computational methods for the performance analysis of vertical axis wind turbines[J].Journal of Solar Energy Engineering,2012,134(3):0310081－03100811.

[14]GB/T 19068.3-2003,離網(wǎng)型風力機發(fā)電機組 第3部分:風洞實驗方法[S].

[15]Kooiman S J,Tullis S.Response of a vertical axis wind turbine to time varying wind conditions found within the urban environment[J].Wind Engineering,2010,34(40):389－401.

[16]Fiedler A J，Tullis S.Blade offset and pitch effects on a high solidity vertical axis wind turbine[J].Wind Engineering,2009,33(3):237－246.

[17]Almohammadi K M,Ingham D B,Ma L,et al.Computational fl uid dynamics (CFD) mesh independency techniques for a straight blade vertical axis wind turbine[J].Energy,2013,58:483－493.

[18]Islam M,Amin M R,Ting D S K,et al.Performance Analysis of a Small Capacity Straight Balded VAWT with Prospective Airfoils[A].46th AIAA Aerospace Science Meeting and Exhibit[C].Reno,2008.

[19]Islam M,Ting D S K,Fartaj A.Desirable airfoil features for small capacity straight bladed VAWT[J].Wind Engineering,2007,31(3):165－196.

[20]茅靖峰,吳愛華,吳國慶,等.垂直軸永磁直驅(qū)風力發(fā)電系統(tǒng)全風速功率控制[J].可再生能源,2014,32 (9) :1339－1345.

[21]梅毅,曲建俊,許明偉.垂直軸風力機葉片動態(tài)失速數(shù)值模擬[J].農(nóng)業(yè)機械學報,2014,45(3):184－190．