李德順, 蘇進(jìn)誠, 李銀然, 郭 濤
(1. 蘭州理工大學(xué) 能源與動(dòng)力工程學(xué)院, 甘肅 蘭州 730050; 2. 蘭州理工大學(xué) 甘肅省風(fēng)力機(jī)工程技術(shù)研究中心, 甘肅 蘭州 730050; 3. 蘭州理工大學(xué) 甘肅省流體機(jī)械及系統(tǒng)重點(diǎn)實(shí)驗(yàn)室, 甘肅 蘭州 730050)
尾流效應(yīng)是影響風(fēng)電場發(fā)電效率的重要因素,而來流湍流度是影響水平軸風(fēng)力機(jī)尾流特性的重要參數(shù)之一.湍流度是風(fēng)速的標(biāo)準(zhǔn)偏差與平均風(fēng)速的比率.受大氣環(huán)境、地表粗糙度等影響,不同風(fēng)場來流湍流度大小及分布不盡相同.因此,研究湍流度對水平軸風(fēng)力機(jī)尾流特性的影響具有重要意義.
近年來,相關(guān)學(xué)者針對風(fēng)力機(jī)尾流做了大量研究.Sorensen等[1]考慮葉片展向和旋轉(zhuǎn)方位角的變化,提出了致動(dòng)線模型.Keck[2]和Ameur等[3]耦合RANS模型和致動(dòng)線模型對風(fēng)力機(jī)尾流進(jìn)行了相關(guān)研究.Ivanell等[4]采用大渦模擬結(jié)合致動(dòng)線的方法對風(fēng)力機(jī)的尾流葉尖渦進(jìn)行了穩(wěn)定性分析.MO等[5]研究了來流速度對風(fēng)力機(jī)尾流穩(wěn)定性的影響,發(fā)現(xiàn)尾渦破碎的位置與來流速度呈函數(shù)關(guān)系.HUN等[6]通過實(shí)驗(yàn)測量的方法,研究了湍流對水平軸風(fēng)力機(jī)尾流及功率特性的影響,發(fā)現(xiàn)湍流來流時(shí)下游風(fēng)力機(jī)的功率輸出小于均勻來流時(shí)的功率輸出,且與風(fēng)力機(jī)縱向間距呈非線性關(guān)系.Talavera等[7]通過風(fēng)洞實(shí)驗(yàn)分別研究了不同湍流來流對單臺及兩臺串列風(fēng)力機(jī)功率和尾流恢復(fù)的影響.王勝軍[8]基于ANSYS FLUENT平臺開發(fā)了風(fēng)力機(jī)風(fēng)輪的致動(dòng)線代碼,研究了風(fēng)力機(jī)不同排列情況下的尾流相互干涉特性.剛蕾等[9]基于Park模型尾流區(qū)線性膨脹假設(shè)和渦黏模型對尾流區(qū)徑向風(fēng)速呈高斯分布假設(shè),提出一種Park-Gauss新模型,對單臺風(fēng)力機(jī)尾流進(jìn)行了數(shù)值模擬.錢耀如等[10]利用致動(dòng)線結(jié)合大渦模擬(LES)的方法研究了低湍流度均勻來流中兩臺串列風(fēng)力機(jī)的氣動(dòng)性能、尾流干擾、功率系數(shù)和推力系數(shù)等,并驗(yàn)證了致動(dòng)線結(jié)合大渦模擬研究風(fēng)力機(jī)尾流的可靠性.章書成等[11]利用高頻PIV系統(tǒng)測量風(fēng)力機(jī)尾流數(shù)據(jù),發(fā)現(xiàn)葉尖渦的耗散速度比中心渦要快,且葉尖渦存在“交互跳躍”現(xiàn)象.侯亞麗等[12]采用大渦模擬方法研究了有、無風(fēng)切入流對風(fēng)力機(jī)尾流湍流特征的影響.張旭耀等[13]通過數(shù)值模擬研究了一臺水平軸風(fēng)力機(jī)尾流的自相似性及流場特性.
本文為進(jìn)一步研究來流湍流度對風(fēng)力機(jī)尾流及轉(zhuǎn)矩特性的影響,采用大渦模擬耦合致動(dòng)線的方法模擬均勻來流及不同湍流度來流條件下的風(fēng)力機(jī)尾流流場,并從尾流速度分布、尾流膨脹及風(fēng)輪轉(zhuǎn)矩三個(gè)方面研究湍流度對風(fēng)力機(jī)尾流及轉(zhuǎn)矩特性的影響規(guī)律.
致動(dòng)線模型由致動(dòng)盤模型改進(jìn)而來,將旋轉(zhuǎn)的葉片簡化為旋轉(zhuǎn)的致動(dòng)線,在葉片上布置若干計(jì)算點(diǎn),并利用BEM理論計(jì)算各計(jì)算點(diǎn)處的氣動(dòng)力,計(jì)算點(diǎn)處的氣動(dòng)力可以表示為
(1)
式中:W為葉片上某翼型的相對入流速度;c為計(jì)算點(diǎn)處翼型弦長;空氣密度為1.225 kg/m3;CL、CD分別為翼型的升力系數(shù)和阻力系數(shù);eL、eD分別為升力和阻力方向的單位向量.
為防止數(shù)值振蕩,需將體積力光滑分布到周圍的網(wǎng)格上,因此采用Gaussian光順函數(shù)將葉素上的作用力光順過渡到周圍網(wǎng)格上.Gaussian光順函數(shù)為
(2)
則計(jì)算點(diǎn)作用于流體的體積力可以表示為
fi=fe?ηε(d)
(3)
式中:fi將作為N-S方程中的源項(xiàng);d為體積力;fi為作用點(diǎn)與葉片力氣動(dòng)力fe作用點(diǎn)之間的距離;ε是光順函數(shù)長度尺度因子,表示高斯分布中調(diào)整葉片作用力集中程度的參數(shù),ε越大,體積力分布越平滑,體積力影響范圍越大.
相較于雷諾平均(RANS)方法,LES方法更容易捕捉各物理量的脈動(dòng)特征.Rethore[14]對比研究了RANS模型和LES模型,發(fā)現(xiàn)LES模型可以更好地捕捉多風(fēng)力機(jī)尾流的摻混及速度恢復(fù),計(jì)算得到的平均速度和湍流場更接近實(shí)驗(yàn)結(jié)果.因此本文基于LES方法研究風(fēng)力機(jī)尾流及載荷特性.
LES方法通過某種濾波函數(shù)將大尺度的渦和小尺度的渦分離,濾波一般采用積分運(yùn)算來實(shí)現(xiàn).在自然風(fēng)況下,流經(jīng)風(fēng)力機(jī)的風(fēng)速低于音速,因此,風(fēng)力機(jī)周圍的氣流可以看作不可壓縮流體.則過濾后的連續(xù)方程可以表示為
(4)
不可壓縮N-S方程為
(5)
(6)
建立亞格子應(yīng)力模型是大渦模擬的核心,本文采用Smagorinsky模型對亞格子應(yīng)力進(jìn)行模擬計(jì)算,該模型為
(7)
(8)
其中,濾波操作采用盒式濾波[15]表達(dá)式如下:
(9)
式中:x′i為小尺度空間坐標(biāo);xi為網(wǎng)格節(jié)點(diǎn)坐標(biāo);Δi為i方向上網(wǎng)格步長.
本文采用蘭州理工大學(xué)景泰實(shí)驗(yàn)基地的一臺33 kW水平軸風(fēng)力機(jī)為研究對象,其主要參數(shù)見表1.
表1 風(fēng)力機(jī)主要參數(shù)
本文采用長方體計(jì)算域,長224 m,寬、高均為36 m,風(fēng)輪位于距計(jì)算域入口24 m處的截面中心,如圖1所示.考慮到計(jì)算資源等各方面因素,將網(wǎng)格尺寸設(shè)為0.5 m,并將如圖所示的陰影區(qū)域加密一次,最小網(wǎng)格尺寸為0.25 m.
圖1 計(jì)算域示意圖Fig.1 Diagram of computational domain
采用OpenFOAM中turbulentIletFvPatchField庫設(shè)置入口速度,平均速度設(shè)為11 m/s,入口速度類型為turbulentIlet,采用Neumann壓力條件;出口邊界為Dirichlet壓力條件和Neumann速度條件;時(shí)間步長0.01 s.通過在來流方向添加隨機(jī)擾動(dòng)產(chǎn)生脈動(dòng)風(fēng)速,對應(yīng)的表達(dá)式為
(10)
式中:χP和χref分別表示擾動(dòng)速度添加網(wǎng)格位置和網(wǎng)格參考點(diǎn);n表示時(shí)間步;α為擾動(dòng)值,取值0.1;s表示擾動(dòng)強(qiáng)度,與不同網(wǎng)格點(diǎn)處的脈動(dòng)尺度有關(guān);CRMS表示REM系數(shù)
(11)
本文通過改變x,y,z方向上的脈動(dòng)尺度來模擬不同湍流度來流,脈動(dòng)尺度為(0.02,0.01,0.01)、(0.2,0.1,0.1)和(0.5,0.2,0.2)時(shí)對應(yīng)得到湍流度為0.85%、5.2%和10.3%.
為研究湍流度對風(fēng)力機(jī)尾流速度分布特性的影響,本文對均勻來流及來流湍流度為0.85%、5.2%及10.3%時(shí)的風(fēng)力機(jī)尾流軸向速度分布進(jìn)行分析.圖2為均勻來流及來流湍流度為0.85%、5.2%及10.3%時(shí)的風(fēng)力機(jī)尾流軸向速度云圖,圖中X為風(fēng)輪下風(fēng)向某斷面與風(fēng)輪平面的軸向距離,D為風(fēng)輪直徑.由圖可知,在均勻來流及湍流來流條件下,尾流區(qū)均產(chǎn)生了明顯的速度虧損(尾流軸向時(shí)均速度與來流軸向時(shí)均速度的差值,時(shí)均速度表示120 s內(nèi)軸向速度時(shí)均值).均勻來流時(shí),尾流邊界在整個(gè)流場內(nèi)比較規(guī)則,這是因?yàn)榫鶆騺砹鲿r(shí)尾流與周圍流動(dòng)摻混較慢;受風(fēng)輪旋轉(zhuǎn)影響,在風(fēng)輪后0~2D范圍內(nèi),尾流速度存在明顯的波動(dòng),隨后逐漸趨于平穩(wěn),這是因?yàn)轱L(fēng)輪后產(chǎn)生了周期性的脫落渦,隨后不斷耗散,使得周期性波動(dòng)逐漸消失.來流湍流度為0.85%時(shí),在0~1.5D范圍內(nèi)存在明顯波動(dòng),隨后逐漸趨于平穩(wěn),從7D開始,尾流邊界變得不規(guī)則,這是因?yàn)橥牧魑擦飨蚝蟀l(fā)展的過程中,尾渦逐漸耗散,尾流不斷與周圍流動(dòng)發(fā)生摻混.來流湍流度為5.2%時(shí),尾流邊界在距風(fēng)輪平面4D附近就已變得不規(guī)則.來流湍流度為10.3%時(shí),尾流邊界在距風(fēng)輪平面2D附近就已變得不規(guī)則,尾流速度在距風(fēng)輪平面13D處已基本恢復(fù)到來流風(fēng)速.可以看出,隨著來流湍流度增大,尾流與周圍流動(dòng)的能量交換加強(qiáng),摻混速度加快,尾流恢復(fù)加快.
圖2 軸向速度云圖Fig.2 Axial velocity contours
圖3為均勻來流與來流湍流度為0.85%、5.2%及10.3%時(shí)的風(fēng)力機(jī)尾流速度廓線及尾流邊界的對比圖,圖中Y為距風(fēng)輪軸線的徑向距離,R為風(fēng)輪半徑,v為尾流軸向時(shí)均速度,v0為來流速度,黑色、紅色虛線分別代表均勻來流和湍流來流時(shí)的尾流邊界.可以看出,在不同來流條件下,風(fēng)輪后均出現(xiàn)明顯的速度虧損,近尾流區(qū)速度廓線呈W型分布,隨著流動(dòng)向后發(fā)展,逐漸趨于V型分布[16].尾流速度最大虧損出現(xiàn)在風(fēng)輪后1D附近.均勻來流時(shí),尾流速度最大虧損率(尾流速度最大虧損值與來流速度的比值)為26%;來流湍流度為0.85%、5.2%及10.3%時(shí),尾流速度最大虧損率分別為23.8%、24.1%及23.9%.這說明湍流會(huì)使尾流速度最大虧損率減小,但湍流度大小對尾流速度最大虧損率影響很小.均勻來流時(shí),在風(fēng)輪后7D處尾流速度恢復(fù)為來流風(fēng)速的75%;來流湍流度分別為0.85%、5.2%及10.3%時(shí),在風(fēng)輪后7D處,尾流速度分別恢復(fù)至來流風(fēng)速的78%、82%及83%.這說明來流湍流度會(huì)加速尾流速度的恢復(fù),原因在于,湍流來流時(shí),風(fēng)力機(jī)尾流與周圍流場的能量交換比均勻來流時(shí)更強(qiáng),湍流摻混速度更快.定義尾流半徑為尾流中心線到尾流邊界(軸向速度為來流速度99%處)的距離,可以看到,尾流半徑在風(fēng)輪后0~1D范圍內(nèi)的增長速度逐漸減小,在1D之后尾流半徑基本呈線性增長.在風(fēng)輪后7D處,均勻來流時(shí),尾流半徑相對風(fēng)輪半徑增大19.5%;來流湍流度為0.85%、5.2%及10.3%時(shí),尾流半徑相對風(fēng)輪半徑增大26.8%、53.6%及78%.這說明隨著來流湍流度增大,尾流區(qū)的膨脹程度增大,主要原因在于,隨著來流湍流度的增大,尾流與周圍流動(dòng)的能量交換加強(qiáng),摻混速度加快,導(dǎo)致尾流徑向速度增大所引起.
圖3 均勻來流與不同湍流度來流尾流速度廓線及尾流邊界對比Fig.3 Comparisons of velocity profile and wake boundary between uniform inflow and turbulent inflow
來流湍流度對風(fēng)輪轉(zhuǎn)矩(風(fēng)輪軸承受的轉(zhuǎn)矩)影響較大,因此本節(jié)對不同來流湍流度條件下的風(fēng)輪轉(zhuǎn)矩特性開展研究.圖4為均勻來流,0.85%、5.2%及10.3%湍流度來流時(shí),風(fēng)輪轉(zhuǎn)矩的時(shí)程變化,紅色水平線表示時(shí)均值.由圖可知,均勻來流時(shí),風(fēng)輪轉(zhuǎn)矩隨時(shí)間呈現(xiàn)周期性波動(dòng),時(shí)均值為5 142 N·m,波動(dòng)幅值為95 N·m,波動(dòng)呈現(xiàn)明顯的周期性變化.來流湍流度為0.85%、5.2%及10.3%時(shí),風(fēng)輪轉(zhuǎn)矩時(shí)均值依次減小為4 807.6、4 762、4 701 N·m,這說明當(dāng)風(fēng)速為11 m/s時(shí),隨著湍流度增大,風(fēng)輪轉(zhuǎn)矩減小,風(fēng)輪軸功率減小,該結(jié)論與文獻(xiàn)[17]一致.來流時(shí)均風(fēng)速不變,風(fēng)輪轉(zhuǎn)速不變,風(fēng)輪軸功率降低,說明風(fēng)輪的風(fēng)能利用系數(shù)降低,即湍流度增大導(dǎo)致風(fēng)輪的風(fēng)能利用系數(shù)降低.來流湍流度為0.85%、5.2%及10.3%時(shí),風(fēng)輪轉(zhuǎn)矩波動(dòng)幅值依次為163、1 027、1 606 N·m,這說明隨著湍流度增大,風(fēng)輪轉(zhuǎn)矩波動(dòng)幅度增大,原因在于,隨著來流湍流度的增大,風(fēng)輪平面內(nèi)流動(dòng)各向異性增強(qiáng),導(dǎo)致風(fēng)力機(jī)葉片氣動(dòng)力波動(dòng)加劇.
圖4 風(fēng)輪轉(zhuǎn)矩時(shí)程圖Fig.4 Time-history diagram of wind turbine torque
功率譜是信號的自相關(guān)函數(shù)的傅里葉變換,描述信號的能量特征隨頻率的變化關(guān)系,可以用來分析不同頻率分量所攜帶的能量大小.圖5為均勻來流與不同湍流度來流條件下風(fēng)輪轉(zhuǎn)矩的功率譜對比.可以看出,相較于均勻來流,湍流來流時(shí)風(fēng)輪轉(zhuǎn)矩的低頻和高頻能量增大,且湍流度越大能量越大,這是由于湍流脈動(dòng)增大所導(dǎo)致.對比低頻與高頻能量分布,發(fā)現(xiàn)低頻段能量增大得更快,原因在于湍流度增大時(shí)湍流尺度增大.均勻來流時(shí)風(fēng)輪轉(zhuǎn)矩功率譜呈現(xiàn)出很強(qiáng)的周期性,原因是均勻來流條件下風(fēng)輪轉(zhuǎn)矩的波動(dòng)主要是由于尾流的誘導(dǎo)引起的.湍流來流時(shí),功率譜出現(xiàn)了對應(yīng)偶數(shù)倍轉(zhuǎn)頻的波峰:湍流度為0.85%時(shí),兩倍轉(zhuǎn)頻處出現(xiàn)波峰;湍流度為5.2%時(shí),兩倍、四倍轉(zhuǎn)頻處均出現(xiàn)波峰;湍流度為10.3%時(shí),兩倍、四倍轉(zhuǎn)頻處波峰對應(yīng)的能量更高.可見,隨著湍流度增大,偶數(shù)倍轉(zhuǎn)頻處出現(xiàn)波峰的現(xiàn)象更加明顯,原因是湍流度增大時(shí),風(fēng)輪范圍內(nèi)風(fēng)速分布的不對稱性增強(qiáng),由于風(fēng)輪為兩葉片,因此,風(fēng)輪轉(zhuǎn)矩呈現(xiàn)出偶數(shù)倍轉(zhuǎn)頻的波動(dòng)增強(qiáng).
圖5 風(fēng)輪轉(zhuǎn)矩功率譜對比圖Fig.5 Comparison of torque power spectrum of wind turbine
1) 均勻來流時(shí),尾流恢復(fù)較慢;湍流來流時(shí),風(fēng)力機(jī)尾流與周圍流場能量交換加強(qiáng),湍流摻混速度變快,從而使得尾流恢復(fù)加快,且湍流度越大,恢復(fù)越快.同時(shí),湍流會(huì)導(dǎo)致風(fēng)力機(jī)尾流速度最大虧損率減小,但湍流度的大小對尾流速度最大虧損率影響很小.
2) 均勻來流時(shí),尾流區(qū)膨脹程度較小,在風(fēng)輪后7倍直徑處,尾流半徑相比風(fēng)輪半徑增大19.5%.來流湍流度為0.85%、5.2%及10.3%時(shí),尾流半徑相比風(fēng)輪半徑分別增大26.8%、53.6%及78%,表明湍流會(huì)使尾流膨脹程度增大,且湍流度越大,膨脹程度越大.
3) 風(fēng)速為11 m/s時(shí),隨湍流度增大,風(fēng)輪轉(zhuǎn)矩時(shí)均值減小,波動(dòng)增大;來流湍流度增大,風(fēng)輪轉(zhuǎn)矩的高頻與低頻部分能量增大,且低頻部分能量增大地更快.