基于CFD的旋轉(zhuǎn)風(fēng)輪氣動(dòng)性能分析

李少華, 岳巍澎, 王東華, 岳 征

(1.東北電力大學(xué)能源與動(dòng)力工程學(xué)院,吉林 132012;2.張家口出入境檢驗(yàn)檢疫局,張家口 075000)

風(fēng)力發(fā)電機(jī)是將風(fēng)能轉(zhuǎn)換為電能的機(jī)械裝置,葉片以及風(fēng)輪的氣動(dòng)性能直接決定了風(fēng)力發(fā)電機(jī)組的效率.由葛勞渥旋渦理論可知:空氣流過(guò)旋轉(zhuǎn)葉片時(shí)會(huì)形成中心渦、附著渦和螺旋后緣尾渦[1].由于渦系對(duì)風(fēng)輪氣動(dòng)性能有較大干擾,因此風(fēng)輪中葉片所處的流場(chǎng)特性分布與單葉片有明顯區(qū)別,旋轉(zhuǎn)葉片渦的脫落對(duì)旋轉(zhuǎn)下游葉片造成很大影響,導(dǎo)致風(fēng)輪下風(fēng)向流場(chǎng)周圍環(huán)境發(fā)生變化.

目前,國(guó)內(nèi)外有關(guān)大型風(fēng)力發(fā)電機(jī)風(fēng)輪流場(chǎng)的研究有試驗(yàn)和數(shù)值模擬兩種方法:試驗(yàn)方面,有關(guān)研究者[2-4]多采用PIV粒子圖像等技術(shù)測(cè)量風(fēng)力發(fā)電機(jī)葉尖及尾流區(qū)的流場(chǎng)工況,但試驗(yàn)成本較高,對(duì)環(huán)境要求嚴(yán)格,且試驗(yàn)研究本身針對(duì)性強(qiáng)、局限性大、不適合觀察和分析旋轉(zhuǎn)風(fēng)輪的整體流場(chǎng)特性;數(shù)值模擬方面,郭婷婷等[5]和吳殿文[6]分析了千瓦級(jí)風(fēng)力發(fā)電機(jī)二維翼型失速特性,并獲得了旋渦產(chǎn)生、發(fā)展和脫落規(guī)律;包飛[7]采用SST k-ω模擬單葉片周圍流場(chǎng),對(duì)葉片截面所在平面的流場(chǎng)域進(jìn)行了分析,并系統(tǒng)研究了壓力面與吸力面的壓差以及升力系數(shù);S?RENSEN 等[8]和SEZER-UZOL等[9]對(duì)千瓦級(jí)風(fēng)力發(fā)電機(jī)風(fēng)輪進(jìn)行了模擬,并重點(diǎn)描述了風(fēng)輪旋轉(zhuǎn)平面的流動(dòng)特性.

到目前為止,國(guó)內(nèi)外系統(tǒng)分析兆瓦級(jí)風(fēng)力發(fā)電機(jī)旋轉(zhuǎn)風(fēng)輪流場(chǎng)特性以及介紹旋轉(zhuǎn)葉片間相互影響方面的文獻(xiàn)較少.筆者對(duì)1.2 MW風(fēng)力發(fā)電機(jī)旋轉(zhuǎn)風(fēng)輪和單葉片流場(chǎng)進(jìn)行了數(shù)值模擬,經(jīng)迭代計(jì)算收斂后,截取葉片周向曲面,等長(zhǎng)度拉伸成平面,并按周期性邊界進(jìn)行了排列,重點(diǎn)研究了旋轉(zhuǎn)葉片運(yùn)動(dòng)軌跡所在平面的流場(chǎng)分布,并通過(guò)觀察湍動(dòng)能場(chǎng)渦脫落云圖分析了葉片間的相互影響.

1 數(shù)學(xué)模型

假設(shè)葉片為剛體且模擬過(guò)程不考慮葉片表面的變形.基于穩(wěn)態(tài)不可壓縮流動(dòng)三維定常雷諾時(shí)均NS方程(RANS)進(jìn)行了數(shù)值模擬,采用Segregated隱式求解器三維穩(wěn)態(tài)算法,紊流模型采用SST k-ω模型,壓力-速度耦合使用Simple算法,對(duì)流項(xiàng)差分格式則采用二階迎風(fēng)格式[10].

式中:φ為通用變量;u為速度;ρ為流質(zhì)密度;Γ為廣義擴(kuò)散系數(shù);S為廣義源項(xiàng).

2 數(shù)值計(jì)算方法

2.1 幾何模型建立與網(wǎng)格劃分

采用Fluent軟件的前處理軟件Gambit進(jìn)行了幾何建模,并對(duì)基于Matlab自主編程設(shè)計(jì)的1.2 MW風(fēng)力發(fā)電機(jī)葉片進(jìn)行了模擬[11],其設(shè)計(jì)風(fēng)速為11.26 m/s,風(fēng)輪額定轉(zhuǎn)速為18.44 r/min.4組葉片均采用相同的流場(chǎng)區(qū)域.圖1為風(fēng)輪流場(chǎng)分區(qū)的主要尺寸及其網(wǎng)格劃分.在圖1中,葉片長(zhǎng)為35 m,取輪轂直徑為3.25 m.

筆者對(duì)本文中風(fēng)力發(fā)電機(jī)小域中葉片與輪轂進(jìn)行了網(wǎng)格局部加密處理,并采用非結(jié)構(gòu)網(wǎng)格劃分,而對(duì)流場(chǎng)的其他部分則采用了相對(duì)較稀疏的結(jié)構(gòu)網(wǎng)格.為了保證網(wǎng)格質(zhì)量,筆者還進(jìn)行了分區(qū)劃分.圖2為不同葉片數(shù)量的風(fēng)輪網(wǎng)格示意圖.由于在對(duì)比不同數(shù)量葉片的葉輪時(shí),若保持每個(gè)葉片近壁面貼體網(wǎng)格數(shù)量一致,則會(huì)導(dǎo)致整體網(wǎng)格的數(shù)量相差很多.在模擬過(guò)程中,將速度和湍流度變化梯度較大的部分采用自適應(yīng)網(wǎng)格處理,每次提高6～8萬(wàn)網(wǎng)格后繼續(xù)迭代計(jì)算,直至輸出功率不再明顯提高為止.

圖1 風(fēng)輪流場(chǎng)分區(qū)的主要尺寸及其網(wǎng)格劃分(單位:m)Fig.1 Partition,main dimension and g rid of flow field(unit:m)

圖2 風(fēng)輪的網(wǎng)格劃分Fig.2 Grid division of wind wheel

2.2 邊界條件的設(shè)置

進(jìn)口:選擇速度進(jìn)口邊界條件,假定進(jìn)口處具有相同的風(fēng)速,不考慮風(fēng)速變化的影響,進(jìn)氣方向垂直于進(jìn)口面,給出進(jìn)口風(fēng)速.

出口:假定出口流動(dòng)充分發(fā)展,此處選出口邊界為自由出流.

葉片及輪轂壁面無(wú)滑移且絕熱,設(shè)定為旋轉(zhuǎn)壁面,相對(duì)附近區(qū)域的旋轉(zhuǎn)速度為零.

3 性能計(jì)算與誤差分析

根據(jù)下式計(jì)算風(fēng)力發(fā)電機(jī)轉(zhuǎn)子的輸出功率.

式中:P為輸出功率,W;T為轉(zhuǎn)矩,N?m;n為葉輪的轉(zhuǎn)速,r/min;b為葉片數(shù).

在得到各葉片的輸出功率后,根據(jù)式(3)可獲得相對(duì)誤差,再根據(jù)相對(duì)誤差的大小分析風(fēng)輪各葉片及單葉片對(duì)風(fēng)力發(fā)電機(jī)出力的影響.

式中:γ為相對(duì)誤差,%;P0=1.2 MW,為設(shè)計(jì)功率.

根據(jù)式(2)和式(3),結(jié)合數(shù)值模擬結(jié)果可獲得轉(zhuǎn)矩.表1為不同葉片數(shù)的模擬結(jié)果.從表1可看出:隨著葉片數(shù)的增多,風(fēng)力發(fā)電機(jī)轉(zhuǎn)子的輸出功率逐漸減小.這是由于旋轉(zhuǎn)上游的葉片發(fā)生渦脫落,使得旋轉(zhuǎn)下游葉片周圍流場(chǎng)的環(huán)境發(fā)生改變,導(dǎo)致其輸出功率減小.可見(jiàn),葉片數(shù)越多,葉片間距離越近,葉片間的相互影響就越大.其中,三葉片風(fēng)輪的輸出功率最接近設(shè)計(jì)輸出功率,這是由于葉片翼型的所有參數(shù)均是按目前普遍使用的兆瓦級(jí)三葉片風(fēng)輪機(jī)設(shè)計(jì)的.四葉片風(fēng)輪的總輸出功率最大,但是由于其成本較昂貴,因此在實(shí)際中不考慮使用.

表1 不同葉片數(shù)的模擬結(jié)果Tab.1 Simulated results for wind wheels with different numbers of blade

4 后處理與流場(chǎng)分析

4.1 后處理技術(shù)

圖3為風(fēng)輪截面的后處理示意圖.在圖3中,筆者在單葉片以及風(fēng)輪流場(chǎng)中分別截取r/R為46%截面的柱面,并利用 Tecplot后處理軟件將圓周曲面拉伸成長(zhǎng)度與曲面周長(zhǎng)相等的平面,然后將單葉片圓周平面按周期性邊界條件復(fù)制,得到不考慮葉片間損失的三葉片翼型后流場(chǎng),再將相同半徑的單葉片和風(fēng)輪流場(chǎng)導(dǎo)入Tecplot后處理軟件進(jìn)行比較和分析.這種后處理方法的優(yōu)點(diǎn)是:沿葉片的運(yùn)動(dòng)軌跡取截面研究,可清晰地看到葉片之間的相互影響.

圖3 風(fēng)輪截面的后處理示意圖Fig.3 Wind turbine section for post processing

4.2 結(jié)果與分析

4.2.1 不同數(shù)目葉片的風(fēng)輪速度比較

圖4為不同數(shù)目葉片風(fēng)輪的旋轉(zhuǎn)平面速度比較云圖.從圖4可看出:單葉片和雙葉片風(fēng)輪的葉片周圍速度分布基本相似,其葉片吸力面均存在同樣的低速區(qū)和增長(zhǎng)區(qū),且壓力面附近隨轉(zhuǎn)動(dòng)軌跡的速度變化不明顯.與其他3種風(fēng)輪葉片相比,三葉片風(fēng)輪葉片間的速度分布最理想,其兩葉片間的速度大小分布均勻,隨著速度由快變慢,速度帶逐漸變寬.四葉片風(fēng)輪的速度分布不均勻,靠近葉片部分的速度變化梯度較大,而兩葉片的中間部分速度變化梯度較小.

圖4 不同葉片數(shù)風(fēng)輪的旋轉(zhuǎn)平面速度比較云圖Fig.4 Velocity comparison of rotating flat among wind wheels with different numbers of blade

4.2.2 不同葉片風(fēng)輪的湍動(dòng)能比較

圖5為r/R等于46%曲面處不同葉片數(shù)目風(fēng)輪旋轉(zhuǎn)截面的湍動(dòng)能對(duì)比,其截面生成方式及截面范圍見(jiàn)圖3,即取風(fēng)輪上風(fēng)向5 m到下風(fēng)向15 m之間部分進(jìn)行分析.將圖5中三維環(huán)面區(qū)域流場(chǎng)簡(jiǎn)化為二維平面問(wèn)題進(jìn)行分析比較,可清楚觀察到不同情況下湍流強(qiáng)度分布規(guī)律及差別.

圖5 不同葉片數(shù)風(fēng)輪的旋轉(zhuǎn)截面湍動(dòng)能對(duì)比Fig.5 Comparison of the turbulent kinetic energy on circular cut of wind wheels with different numbers of blade

單葉片流場(chǎng)域獨(dú)立,無(wú)其他葉片尾流影響.尾流湍動(dòng)能在2 m2/s2以上,其長(zhǎng)度為10 m,與旋轉(zhuǎn)平面約成25°,之后出現(xiàn)湍動(dòng)能較穩(wěn)定的帶狀區(qū)域.在(15,5)處,開(kāi)始出現(xiàn)變化梯度較大的低湍動(dòng)能尾跡渦,并向下風(fēng)向擴(kuò)散延伸.3種多葉片風(fēng)輪的翼型流場(chǎng)相似,但它們的湍流強(qiáng)度與影響區(qū)域的面積不同.雙葉片風(fēng)力發(fā)電機(jī)葉片周圍高湍區(qū)域的湍動(dòng)能大于單葉片,其影響面積略小于單葉片,葉片避開(kāi)了另一葉片的尾流區(qū),即各葉片間的影響較小,尾流遠(yuǎn)場(chǎng)出現(xiàn)低湍動(dòng)能的尾跡渦.三葉片風(fēng)輪葉片間的相互影響很大,葉片尾緣的湍動(dòng)能帶出現(xiàn)分離,即附著渦平穩(wěn)脫落,壓力面脫落湍動(dòng)能沿x軸切向延伸至旋轉(zhuǎn)下游的葉片吸力面,湍動(dòng)能在2～3 m2/s2,吸力面脫落湍動(dòng)能與旋轉(zhuǎn)上游葉片切向延伸的湍動(dòng)能匯合后偏y軸延伸,此段湍動(dòng)能變化明顯,由5 m2/s2降為3 m2/s2,影響范圍在y=0到y(tǒng)=10 m間.四葉片風(fēng)輪葉片均處于旋轉(zhuǎn)上游的葉片尾流陰影中,葉輪間影響作用最大,但由于間隔較小,風(fēng)輪上游及下風(fēng)向y＞14 m流場(chǎng)區(qū)域流動(dòng)平穩(wěn),沒(méi)有出現(xiàn)明顯的旋渦.

圖6為風(fēng)輪旋轉(zhuǎn)平面湍動(dòng)能的比較.由圖6可知:隨著葉片數(shù)的增加,葉片間湍動(dòng)能逐漸增大.單葉片和雙葉片旋轉(zhuǎn)風(fēng)輪葉片附近的流場(chǎng)基本一致,即雙葉片風(fēng)輪葉片間影響可以忽略不計(jì).三葉片風(fēng)輪葉片處于旋轉(zhuǎn)上游葉片脫落附著渦衰減后的平穩(wěn)部分,受影響較小.四葉片風(fēng)輪的流場(chǎng)較混亂,葉片周圍湍動(dòng)能較大.

圖6 不同葉片數(shù)風(fēng)輪旋轉(zhuǎn)平面湍動(dòng)能的比較Fig.6 Comparison of the turbulent kinetic energy on rotating flat of wind wheels with different numbers of blade

5 結(jié) 論

(1)由于旋轉(zhuǎn)上游的葉片發(fā)生渦脫落,導(dǎo)致旋轉(zhuǎn)下游葉片周圍流場(chǎng)的環(huán)境發(fā)生變化.隨著葉片數(shù)的增加,旋轉(zhuǎn)風(fēng)輪中單個(gè)葉片的輸出功率逐漸減小,風(fēng)力發(fā)電機(jī)效率與運(yùn)行質(zhì)量都受到影響.

(2)多葉片旋轉(zhuǎn)風(fēng)輪背風(fēng)側(cè)流場(chǎng)的平均湍動(dòng)能比單葉片和雙葉片風(fēng)輪大.脫落后的附著渦通過(guò)旋轉(zhuǎn)下游葉片時(shí)沒(méi)有形成類似單葉片后的較大旋渦,因而對(duì)風(fēng)輪下風(fēng)向的穩(wěn)定性有利.

[1]HANSEN M O L.Areodynamics of wind turbines[M].London,UK:Earthscan Science Publishers,2008:18-26.

[2]胡丹梅,田杰,杜朝輝.水平軸風(fēng)力機(jī)尾跡的測(cè)量與分析[J].動(dòng)力工程,2006,26(5):751-760.HU Danmei,TIAN Jie,DU Zhaohui.Measurement and analysis of wake behind horizontally orientated air turbines[J].Journal of Power Engineering,2006,26(5):751-760.

[3]BART HELEMIE R J,RATHMANN O,FRANDSEN S T,et al.Modelling and measurements of wakes in large wind farms[J].Journal of Physics:Conference Series,2007,75(1):12-49.

[4]DOBREV I,MAALOUF B,TROLDBORG N,et al.Investigation of the wind turbine vortex structure[C]//14th Int Symp on Applications of Laser Techniques to Fluid Mechanics.Lisbon,Portugal:[s.n.],2008.

[5]郭婷婷,吳殿文,王成蔭,等.風(fēng)力發(fā)電機(jī)葉片預(yù)彎設(shè)計(jì)及其數(shù)值研究[J].動(dòng)力工程學(xué)報(bào),2010,30(6):450-455.GUO Tingting,WU Dianwen,WANG Chengyin,et al.Prebend design and numerical simulation of wind turbo-generator[J].Journal of Chinese Society of Power Engineering,2010,30(6):450-455.

[6]吳殿文.風(fēng)力機(jī)優(yōu)化設(shè)計(jì)及氣動(dòng)性能的數(shù)值模擬[D].吉林:東北電力大學(xué)能源與動(dòng)力工程學(xué)院,2010.

[7]包飛.風(fēng)力機(jī)葉片幾何設(shè)計(jì)與空氣動(dòng)力學(xué)仿真[D].大連:大連理工大學(xué)機(jī)械工程學(xué)院,2008.

[8]S?RENSEN N N,MICHELSEN J A.Drag prediction for blades at high angle of attack using CFD[J].Journal of Solar Energy Engineering,2004,126(4):1011-1016.

[9]SEZER-UZOL N,LONG LN.3-D time-accurate CFD simulations of wind turbine rotor flow fields[C]//44th AIAA Aerospace Sciences Meeting and Exhibit.Reno,Nevada:American Institute of Aeronautics and Astronautics,2006:6-19.

[10]何顯富,盧霞,楊躍進(jìn),等.風(fēng)力機(jī)設(shè)計(jì)、制造與運(yùn)行[M].北京:化學(xué)工業(yè)出版社,2009:189-196.

[11]GUO Tingting,WU Dianwen,XU Jihui.The method of large-scale wind turbine blades design based on Matlab programming[C]//2009 Supergen Conference.Nanjing,China:[s.n.],2009.