四葉片小型水平軸風(fēng)力機尾流特征的數(shù)值模擬

辛大波, 李佳宇, 張洪福, 薛志成, 王 健

(1. 東北林業(yè)大學(xué) 土木工程學(xué)院, 哈爾濱 150040; 2. 廣東石油化工學(xué)院 建筑工程學(xué)院, 廣東 茂名 525000)

風(fēng)能是一種清潔、豐富的可再生資源,在全球變暖和資源枯竭的背景下,風(fēng)力發(fā)電變得越來越重要.在實際工程中,風(fēng)力機總是在某一風(fēng)能資源豐富的區(qū)域分布式排列.當(dāng)來流風(fēng)受到上游風(fēng)力機的干擾作用時,風(fēng)速降低,脈動成分增強,將嚴(yán)重影響下游風(fēng)力機的性能,主要體現(xiàn)在兩個方面:1)上游尾流區(qū)域的速度損失降低了下游風(fēng)力機的輸出功率;2)尾流區(qū)域的高湍流會導(dǎo)致下游風(fēng)力機的動態(tài)負(fù)荷增加,影響風(fēng)力機的使用壽命.因此,研究風(fēng)力機尾流特征對于有效收集風(fēng)能以及保證風(fēng)力機的結(jié)構(gòu)安全至關(guān)重要.

目前,針對風(fēng)力機尾流的研究主要有兩種方法:風(fēng)洞試驗和數(shù)值模擬.風(fēng)洞試驗可以準(zhǔn)確地獲得速度損失和尾流湍流特性[1-3].但由于實驗室尺寸條件的限制,只能開展風(fēng)力機近尾流區(qū)域的流動特性的相關(guān)研究,對遠(yuǎn)尾流區(qū)域的速度損失和湍流特性的研究鮮有文獻報道.與風(fēng)洞試驗相比,基于計算流體動力學(xué)(CFD)的數(shù)值模擬方法[4-6]不受尺寸的影響,且可同時獲取全區(qū)域的速度場分布,發(fā)展?jié)摿薮?Sedaghatizadeh等[7]研究了LES數(shù)值模擬方法在風(fēng)力機流場分析方面的可行性,結(jié)果表明,LES可以精確地模擬風(fēng)力機尾流特征;Hornsh?j-M?ller等[8]采用CFD方法精確地模擬了風(fēng)力機尾流速度虧損和湍流強度;Mo等[9]結(jié)合大渦模擬研究了風(fēng)力機尾流特性,指出5倍風(fēng)力機直徑是近、遠(yuǎn)尾流的分界線,且風(fēng)速越大,近尾流區(qū)的長度越長.盡管CFD方法廣泛用于預(yù)測風(fēng)力機的空氣動力學(xué)性能及風(fēng)力機尾流特征的研究中,但由于高精度的模擬需要解析風(fēng)力機葉片的邊界層特征,計算時間長并且成本較高[10].鑒于此,眾多學(xué)者提出了無需實體風(fēng)力機的致動盤CFD源項模型[11-12].Amini等[13]使用致動盤方法對風(fēng)力機的三維修正空氣動力系數(shù)進行了數(shù)值模擬;任會來等[14]針對現(xiàn)有致動盤模型高估尾流速度的問題,對動量源項進行修正,以正確反映風(fēng)輪對來流的作用,結(jié)果表明,對耗散率源項進行修正可以提高模擬的精度并且可以準(zhǔn)確模擬不同風(fēng)力機排列方式下的尾流特征;劉鑫等[15]將CFD與致動盤方法相結(jié)合,研究了風(fēng)力機尾流發(fā)展特點,驗證了該方法可以精準(zhǔn)模擬風(fēng)力機尾流.

目前,風(fēng)力機尾流的相關(guān)研究主要集中在大型風(fēng)力機上,而很少關(guān)注小型風(fēng)力機設(shè)備.由于小型風(fēng)力機可以從低風(fēng)速中獲取風(fēng)能,在城市市區(qū)、鄉(xiāng)村近郊區(qū)等低風(fēng)速區(qū)域逐漸被廣泛使用.此外,當(dāng)小型水平軸風(fēng)力機在工作狀態(tài)時,葉片旋轉(zhuǎn)會產(chǎn)生順流向渦.該類型渦結(jié)構(gòu)已被證實可有效抑制大跨建筑、大跨橋梁風(fēng)效應(yīng),提高工程結(jié)構(gòu)抗風(fēng)能力.本課題組采用風(fēng)洞試驗驗證了小型水平軸葉片對于抑制大跨橋梁主梁渦激振動的有效性[16].基于此,本文以四葉片小型水平軸風(fēng)力發(fā)電機為研究對象,建立真實風(fēng)力機的全三維模型,并采用大渦模擬研究該風(fēng)力機的近、遠(yuǎn)尾流定量化特征,提出能夠進行快速計算的風(fēng)力機尾渦致動盤模型,通過風(fēng)洞試驗驗證數(shù)值模擬與模型的準(zhǔn)確性.

1 致動盤模型和數(shù)值模擬方法

1.1 致動盤模型

風(fēng)力機的掃掠范圍是圓形的,可將風(fēng)力機簡化為一個圓柱形的致動盤模型.當(dāng)風(fēng)通過風(fēng)輪時風(fēng)速會減小,風(fēng)力機后面的壓力階梯發(fā)生變化,理論示意圖如圖1所示.假設(shè)來流速度為v0,根據(jù)風(fēng)輪機的阻力系數(shù),可以得到軸向力[17]為

圖1 一維動量理論示意圖

(1)

式中:ρ0為空氣密度,ρ0=1.225 kg/m3;CT為阻力系數(shù).應(yīng)用于AD模型的單位體積源項的表達(dá)式為

(2)

式中,Δx為致動盤厚度,Δx=0.2 m.

1.2 數(shù)值模擬方法

大渦模擬方法可直接求解大尺度湍流,對于小尺度湍流,僅利用模型考慮其對大尺度運動的作用,而不直接進行求解.因此,LES需對N-S方程進行空間過濾,過濾后的連續(xù)性方程和N-S方程表達(dá)式[7]為

(3)

(4)

(5)

其中:μ為動力黏度;δij為克羅內(nèi)克符號,當(dāng)i=j時δij=1,i≠j時δij=0.τij為亞格子應(yīng)力(subgrid-scale streese,簡稱SGS應(yīng)力),體現(xiàn)了小尺度渦的運動對所求解運動方程的影響,定義為

(6)

本文采用工程中常用的Smagorinsky-Lilly亞格子模型,計算亞格子應(yīng)力最常用的方法是利用Boussinesq渦粘假設(shè)[18].將τij表示為正應(yīng)力與偏向應(yīng)力,即

(7)

(8)

亞格子應(yīng)力表示為

(9)

式中,ξSGS為亞格子尺度粘度,在Smagorinsky模型中,其表達(dá)式為

(10)

將式(10)代入式(9),則有

(11)

式中:Ck、Ce為Smagorinsky模型無量綱系數(shù),Ck=0.067 6,Ce=0.93;Δ為濾波器過濾尺度.

Smagorinsky模式是根據(jù)唯象論推出的剪切湍流亞格子模型,該模型概念簡單、易于實施,在風(fēng)力機模擬中應(yīng)用較為廣泛[18].故本文采用工程中常見的Smagorinsky-Lilly亞格子模型.

2 模型建立及網(wǎng)格劃分

本文以四葉片小型水平軸風(fēng)力機為研究對象,重點研究該風(fēng)力機葉片的尾流特征,僅對風(fēng)力機葉片(忽略輪轂影響)進行三維建模及數(shù)值模擬分析.尖速比是影響風(fēng)力機氣動特性的關(guān)鍵參數(shù),其表達(dá)式為

(12)

式中:ω為角速度;R為葉片的半徑.本文中風(fēng)力機尖速比取為3.0,則風(fēng)力機葉片的角速度為50 rad/s.風(fēng)力機的尺寸和形狀如圖2所示(單位:m).

圖2 風(fēng)力機尺寸及模型圖

為了有足夠的內(nèi)部空間供尾流發(fā)展,計算區(qū)域為長方體區(qū)域,如圖3所示.該區(qū)域的來流方向沿Z軸計算,寬度沿X軸方向,Y軸沿垂直方向,XOY為旋轉(zhuǎn)平面,YOZ為水平平面.基于風(fēng)力機的中心,上游區(qū)域長度為5D(D為風(fēng)輪機的直徑),下游區(qū)域長度為20D,徑向直徑為5D.計算區(qū)域可分為兩部分,具有動態(tài)網(wǎng)格的風(fēng)力機旋轉(zhuǎn)區(qū)域和長方體的內(nèi)部空間.圓柱區(qū)域和長方體區(qū)域采用“Interface”連接.由于適當(dāng)?shù)木W(wǎng)格質(zhì)量既可以減少計算時間,又可以提高精度,可將立方體分為兩部分,如圖4所示.對于風(fēng)力機旋轉(zhuǎn)區(qū)域,采用混合四面體網(wǎng)格,數(shù)量為180萬個網(wǎng)絡(luò),而對于剩余部分,采用六面體網(wǎng)格劃分,約為481萬個網(wǎng)格.為了保證數(shù)值模擬的準(zhǔn)確性,入口定義為Velocity-inlet,出口定義為均勻壓力出口Outlet,其他表面定義為Symmetry.網(wǎng)格使用Ansys FLUENT17.0進行計算,壓力和速度耦合使用SIMPLEC方法,動量方程采用半隱式Adams-Bashforth和Crank-Nicolson方法進行求解.為保證計算過程的穩(wěn)定性,時間步長取為0.000 1 s,風(fēng)力機和旋轉(zhuǎn)區(qū)域的運動保持靜態(tài).

圖3 計算域示意圖

圖4 網(wǎng)格劃分示意圖

3 模擬驗證

在來流風(fēng)速為10 m/s、風(fēng)機轉(zhuǎn)速為50 rad/s時,選取三種不同網(wǎng)格數(shù)量的模型進行網(wǎng)格無關(guān)性驗證.網(wǎng)格數(shù)分別為275萬、481萬和643萬,以風(fēng)力機扭矩系數(shù)Cm驗證網(wǎng)格無關(guān)性,結(jié)果如表1所示.綜合考慮計算效率和成本,本文最終選用481萬的網(wǎng)格模型進行分析.

表1 網(wǎng)格無關(guān)性驗證

驗證試驗在東北林業(yè)大學(xué)風(fēng)洞實驗室中進行,風(fēng)洞試驗段截面為0.8 m×1.0 m,長度為5 m,風(fēng)速范圍為2～50 m/s,自由湍流強度小于0.5%.為了滿足阻塞率的要求,風(fēng)力機葉片縮尺比為1∶6.7,風(fēng)力機葉片直徑為18 cm.風(fēng)力機尾流測試裝置及結(jié)果分別如圖5～6所示,試驗結(jié)果與模擬結(jié)果基本吻合.圖6中,Y/D為風(fēng)力機位置的相對高度,v/v0為軸向速度的無量綱量.

圖5 風(fēng)力機尾流測試布置圖

圖6 風(fēng)速變化圖

4 結(jié)果與分析

4.1 尾流速度分布和渦結(jié)構(gòu)

圖7為渦結(jié)構(gòu)示意圖.由圖7可知,葉片旋轉(zhuǎn)區(qū)域的流動結(jié)構(gòu)比較復(fù)雜,包括葉尖渦、附著渦和中心渦,其中,葉尖渦流發(fā)生在葉片的頂端,并且隨著風(fēng)力機的旋轉(zhuǎn)逐漸減少,旋轉(zhuǎn)方向與葉片的旋轉(zhuǎn)方向相反,同時不斷向下游移動.附著渦是在葉片表面形成并附著在葉片上的渦結(jié)構(gòu),在邊緣逐漸脫落.由葉片根部和輪轂產(chǎn)生的中心渦以螺旋狀軌跡從輪轂中心向下游擴散,在8D位置后完全消失.圖8為瞬態(tài)軸向渦量云圖.從圖8中可以看出,風(fēng)力機后方沿Z方向有明顯的渦量分布.在Z/D=1時可清晰識別出尾渦結(jié)構(gòu),隨著距離的增加,渦結(jié)構(gòu)逐漸分解(對稱性消失).

圖7 渦結(jié)構(gòu)示意圖

圖8 瞬態(tài)軸向渦量云圖

圖9為風(fēng)力機尾流的軸向速度虧損分布.隨著葉尖渦旋的發(fā)展,尾流速度逐漸降低,之后隨著渦旋消失,尾流速度逐漸增加.

圖9 尾流軸向速度分布

圖10為葉片后方典型位置的風(fēng)速.由圖10可以看出,在葉片后方,分布基本上是軸對稱的,由于葉片上游的風(fēng)速下降,下游的尾流能量減少.隨著尾流渦旋的消失,速度逐漸增大.在Z/D=8之后,隨著下游距離的增加,速度虧損減少,橫截面積擴大,尾流寬度相應(yīng)增加.

圖10 尾流速度分布

4.2 軸向速度分布

圖11為致動盤模型軸向速度分布.在-2.0

圖11 致動盤模型軸向速度分布

與Jiménez等[19]研究的三葉片風(fēng)力機相比,在近尾流區(qū)域,四葉片風(fēng)力機首先達(dá)到最低值,且該值低于三葉片模擬值,但在遠(yuǎn)尾流(Z/D>5)風(fēng)速虧損值比較接近,這表明四葉片風(fēng)力機在風(fēng)環(huán)境改善和結(jié)構(gòu)風(fēng)振控制中更具有優(yōu)勢.從圖11中可以看出,風(fēng)力機從來流中提取能量后恢復(fù)過程所造成的速度損失在-0.58時,尾流效應(yīng)仍然存在,但不明顯,軸向速度基本恢復(fù)到了來流風(fēng)速的95%以上,隨著尾流的速度逐漸恢復(fù)至來流風(fēng)速,尾流完全擴散.

在致動盤的近尾流區(qū)域,特別是在尾流最低風(fēng)速下的中心區(qū)域,速度虧損有很大不同.在Z/D=1處,致動盤將葉片組合成一個簡化的圓盤,忽略了葉片的旋轉(zhuǎn)效應(yīng),致動盤的擾動效應(yīng)使風(fēng)速達(dá)到最低值,虧損最大.在Z/D=3處,速度恢復(fù)至40%,但風(fēng)力機尾流速度卻沒有發(fā)生變化,這說明風(fēng)力機對尾流的干擾作用比致動盤更明顯.而在遠(yuǎn)處的尾流區(qū)域(Z/D>5),由于流場的變化主要取決于自身的發(fā)展,兩者結(jié)果基本一致.綜上所述,致動盤對近尾流區(qū)域的速度有較大影響,在遠(yuǎn)尾流區(qū)域幾乎吻合,因此可以模擬四葉片風(fēng)力機的尾流分布.

4.3 湍流強度分布

圖12為在-2.0

圖12 湍流強度分布

(13)

式中:u′為脈動風(fēng)速;U為平均風(fēng)速.

根據(jù)Chamorro等[20]的研究,低湍流強度的來流風(fēng)通過旋轉(zhuǎn)的風(fēng)力機時,湍流強度將在尾流中明顯增加.從圖12中可以看出,在Z/D=-1處,湍流強度保持在0.01左右.然而,在Z/D=1處,湍流強度最大值達(dá)到了0.3且有兩個湍流強度峰值.在Z/D=3處,也可觀察到兩個湍流強度峰值,最大值約為入流湍流強度的30倍,達(dá)到了0.36.當(dāng)Z/D>5時,葉尖附近區(qū)域的湍流強度峰值開始消失.此外,在Z/D<3的近尾流區(qū)域,湍流強度值隨著與風(fēng)力機后部距離的增加而逐漸增大,可以看出,當(dāng)風(fēng)場中相鄰風(fēng)力機之間的距離小于3D時,計算風(fēng)力機功率時必須考慮湍流對風(fēng)力機輸出的影響;在遠(yuǎn)尾流區(qū)域,由于尾流的不穩(wěn)定性,尾流的流動結(jié)構(gòu)在Z/D=3以后的區(qū)域不斷波動,當(dāng)Z/D>5時,湍流強度在湍流擾動造成的速度損失恢復(fù)過程中開始迅速降低,最終在Z/D=10處,最大湍流強度降至0.18.

5 結(jié) 論

本文以四葉片小型水平軸風(fēng)力機為研究對象,采用致動盤源項模型并結(jié)合大渦模擬獲取風(fēng)力機尾流渦特征.主要結(jié)論如下:

1) 大渦模擬可以較為準(zhǔn)確地獲得四葉片風(fēng)力機的尾流分布,無葉片實物的致動盤模型可以有效預(yù)測尾流流速的虧損過程,大大減少計算資源.

2) 四葉片小型水平軸風(fēng)力機在尾流區(qū)域有明顯速度虧損,與三葉片風(fēng)力機不同的是,在近尾流區(qū)四葉片風(fēng)力機首先達(dá)到最低值,且該值低于三葉片風(fēng)力機.對于湍流度,在近尾流區(qū)出現(xiàn)了兩個湍流強度峰值,遠(yuǎn)尾流區(qū)湍流度逐漸減小.

3) 四葉片小型風(fēng)力機在工作時,葉片旋轉(zhuǎn)區(qū)域的流動結(jié)構(gòu)比較復(fù)雜,在近尾流區(qū)的渦結(jié)構(gòu)主要為葉尖渦與附著渦,而遠(yuǎn)尾流區(qū)只有中心渦.