何煜平,楊建民
(上海交通大學(xué) 海洋工程國家重點(diǎn)實(shí)驗(yàn)室,上海 200240)
所謂海流,是風(fēng)應(yīng)力、密度梯度、引潮力等作用而形成的大規(guī)模相對穩(wěn)定的流動(dòng)。據(jù)統(tǒng)計(jì),全球海流能蘊(yùn)藏總量接近于3.5TW[1-2]。若能合理地將這些能源加以利用,則能為緩解能源危機(jī)提供一條新的道路。
海流能發(fā)電的方式類似于風(fēng)能發(fā)電,從概念上來說都是利用流體產(chǎn)生的升力或者阻力使獲能裝置運(yùn)動(dòng),并帶動(dòng)內(nèi)部的發(fā)電機(jī)發(fā)電。與風(fēng)能相比,海流能發(fā)電具有以下的優(yōu)勢:1)海水密度是空氣密度的830 倍,使得海流能具有比風(fēng)能更大的能量密度;2)海流,特別是潮流,其流向和流量相對穩(wěn)定并且容易預(yù)測[3-4]。
目前,海流能技術(shù)最為先進(jìn)、產(chǎn)業(yè)化進(jìn)程最為順利的國家為英國。2008年4月,英國MCT 公司完成了1.2 MW 的水平軸海流機(jī)發(fā)電裝置“SeaGen”的安裝和運(yùn)營,標(biāo)志著世界上第一個(gè)商業(yè)化規(guī)模的海流發(fā)電系統(tǒng)投入使用[2]。國內(nèi)的部分學(xué)校和研究所也都在積極開展海流能發(fā)電的研究,但尚處于起步階段。
目前海流能領(lǐng)域的研究熱點(diǎn)主要集中在以下方面:1)海流機(jī)葉片的水動(dòng)力設(shè)計(jì)方法;2)葉片可調(diào)距機(jī)構(gòu)的設(shè)計(jì);3)安裝與支撐方式的設(shè)計(jì)研究;4)葉片氣蝕的研究和抑制;5)機(jī)組布局的研究;6)海流能發(fā)電裝置新概念的設(shè)計(jì)[2-4]。
本文利用了葉素-動(dòng)量理論[5-6],對給定額定工況的海流能發(fā)電機(jī)的葉片進(jìn)行了設(shè)計(jì),并分別利用理論方法和CFD 方法對葉片的水動(dòng)性能進(jìn)行了計(jì)算和分析。
葉素-動(dòng)量理論是同時(shí)利用動(dòng)量定理和葉素理論進(jìn)行葉片水動(dòng)力外型設(shè)計(jì)的方法。若葉輪遠(yuǎn)前方流速為U1,遠(yuǎn)后方流速為U2,葉輪轉(zhuǎn)速為Ω,旋轉(zhuǎn)尾流角速度為ω,可定義軸向誘導(dǎo)因子a 和周向誘導(dǎo)因子b:
由動(dòng)量定理,葉輪受到的軸向力T 和轉(zhuǎn)矩M 等于單位時(shí)間內(nèi)通過葉輪的流體動(dòng)量和角動(dòng)量的變化量。對于葉輪上一個(gè)半徑為r、寬度為dr 的圓環(huán),其受到的軸向力dT 和轉(zhuǎn)矩dM 可表示為:
葉素理論則是將葉片沿徑向分成若干葉素,對每個(gè)葉素分別進(jìn)行分析,然后將作用在每個(gè)葉素上的力疊加,得到整個(gè)葉輪上的流體作用力。位于半徑r 處葉素的速度和受力情況如圖1 所示。
葉素上的流體相對速度U 為:
葉輪旋轉(zhuǎn)面與流體相對速度之間的夾角φ 為:
式中:α 為攻角;β 為扭角;λ 為局部速度比,λ =Ωr/U1,其中葉尖處局部速度比稱尖速比,記為λ0。
設(shè)在攻角α 下,翼型的升力系數(shù)和阻力系數(shù)分別為CL和CD。則對于弦長為c,寬度為dr 的葉素,其受到的升力dL 和阻力dD 可分別表示為:
若葉片數(shù)為N,則半徑為r、寬度為dr 的葉輪圓環(huán)上受到的軸向力dT 和轉(zhuǎn)矩dM 可以表示為:
令式(2)和式(7)相等,式(3)和式(8)相等,可得
Tony Burton 通過比較考慮阻力和忽略阻力兩種情形下的最優(yōu)葉片設(shè)計(jì),發(fā)現(xiàn)阻力對葉片的最優(yōu)設(shè)計(jì)影響很?。?]。因此,忽略葉片阻力,則式(9)變?yōu)?
根據(jù)式(5)和式(10),得到a、b 和λ 三者之間的關(guān)系:
圖1 葉素的速度和作用力分析Fig.1 Velocities and forces of a blade element
表征海流機(jī)葉片性能的參數(shù)主要包括推力系數(shù)CT、轉(zhuǎn)矩系數(shù)CM和獲能系數(shù)CP:
海流機(jī)葉片水動(dòng)力外形設(shè)計(jì),是在確定的設(shè)計(jì)流速U1、尖速比λ0和輸出功率P0的條件下,使葉片的獲能系數(shù)CP盡可能大。設(shè)計(jì)時(shí),需要首先確定葉片數(shù)目N、葉輪半徑R、葉轂半徑RH,并選定合適的翼型。目前比較常見的葉片數(shù)主要有2 葉和3 葉。葉片數(shù)越多,獲能系數(shù)CP越大,但軸向力系數(shù)CT也隨之大幅提高,對支撐系統(tǒng)的設(shè)計(jì)有較高的要求。葉輪半徑R 可以通過下式估算:
式中的獲能系數(shù)CP是通過經(jīng)驗(yàn)確定的預(yù)估值。葉轂半徑RH可根據(jù)設(shè)計(jì)經(jīng)驗(yàn),取葉輪半徑合適的百分?jǐn)?shù)。葉片翼型的選取可以參照各翼型的水動(dòng)力性能曲線,選取最大升阻比L/Dmax較大的翼型,同時(shí)還需保證具有較好的失速特性。通常葉尖處翼型較薄,葉根處翼型較厚,并逐漸過渡為輪轂處的圓形截面。
確定了上述設(shè)計(jì)參數(shù)后,可利用BEM 理論確定弦長c 和扭角β 沿葉輪徑向的分布。步驟如下:
1)選取一控制截面,其徑向位置為r。計(jì)算局部速度比λ=Ωr / U1。
2)確定a 和b,使獲能系數(shù)CP最大。根據(jù)式(8),葉輪圓環(huán)處的功率可表示為:
因此,(1 -a)b 的取值越大,dP 就越大。但是,a、b 和λ 三者之間受到了式(11)的約束。為了確定a 和b 的最優(yōu)值,需要求解以下單目標(biāo)帶約束的非線性優(yōu)化問題:
3)根據(jù)控制截面處翼型升力系數(shù)曲線,找到對應(yīng)最大升阻比的最優(yōu)攻角α 以及相應(yīng)的升力系數(shù)CL。
4)根據(jù)式(5)確定φ,則扭角β=φ-α。根據(jù)式(9)確定弦長c。
5)重復(fù)步驟1)~4),直到所有控制截面計(jì)算完畢。
在葉片設(shè)計(jì)時(shí),需要葉片截面翼型的升力系數(shù)曲線和阻力系數(shù)曲線作為輸入?yún)?shù)。對于NACA 翼型族相關(guān)性能曲線,可通過XFoil 等翼型分析軟件獲得,也可通過CFD 方法求得[7]。本文所使用的翼型性能曲線是通過Fluent 軟件求得,其中NACA 2414、NACA 2416 和NACA 2418 翼型的升力系數(shù)曲線和阻力系數(shù)曲線如圖2 所示。翼型失速角約為16°到18°,最大升阻比對應(yīng)的最佳攻角約為7°。
圖2 NACA 24 系列翼型升力系數(shù)和阻力系數(shù)曲線Fig.2 Lift coefficient and drag coefficient curves of NACA 24 series
設(shè)有如表1 所示的設(shè)計(jì)要求,利用BEM 理論進(jìn)行水平軸海流機(jī)葉片的設(shè)計(jì)。設(shè)計(jì)時(shí),為了簡化問題,認(rèn)為海流機(jī)所處水深較深,海流不受自由表面波浪的影響,且認(rèn)為海流流速在整個(gè)葉輪處均勻分布。由于海流相對比較穩(wěn)定,在通常情形下流速能保持在較小范圍內(nèi)浮動(dòng),因此設(shè)計(jì)流速可以從該流速范圍內(nèi)選取。
表1 海流機(jī)葉片設(shè)計(jì)參數(shù)Tab.1 Parameters of blade designing of a current turbine
取獲能系數(shù)預(yù)估值為0.4,根據(jù)式(13)計(jì)算葉片半徑,得R≈5.5 m。葉片翼型選取NACA 24 系列,葉根處為了保證足夠的結(jié)構(gòu)強(qiáng)度,選取葉片厚度t 為18%,然后向葉尖處逐漸過渡為16%和14%。利用前述設(shè)計(jì)步驟計(jì)算截面弦長c 和扭角β,c 和β 沿著葉輪徑向的分布曲線如圖3 所示。利用Pro/E 進(jìn)行葉片實(shí)體建模,如圖4 所示。
圖3 150 kW 海流機(jī)葉片幾何參數(shù)徑向分布Fig.3 Radial distribution of geometric parameters of the 150 kW current turbine's blade
圖4 150 kW 海流機(jī)葉片實(shí)體模型Fig.4 Model of the 150 kW current turbine's blade
為了檢驗(yàn)葉片設(shè)計(jì)是否合理,以及預(yù)測不同工況下的性能,需對海流機(jī)的水動(dòng)力性能進(jìn)行計(jì)算。BEM理論結(jié)合普朗特(Prandtl)修正以及葛勞渥(Glauert)修正是較常用的海流機(jī)水動(dòng)力性能理論計(jì)算方法[5-6]。
葉素-動(dòng)量理論忽略了順著翼展方向的速度分量,然而實(shí)際上水流在葉片上存在著徑向流動(dòng)。尤其在葉尖和輪轂處,這種三維流動(dòng)效應(yīng)更是明顯,因此葉尖損失和輪轂損失對葉片性能的影響不容忽略。Prandtl針對該現(xiàn)象提出了Prandtl 漸進(jìn)法。根據(jù)該理論,葉尖損失和輪轂損失可用損失因子Ft(r)和Fh(r)來表示:
總的損失因子F (r)可表示為:
經(jīng)過修正后的軸向誘導(dǎo)因子a 和軸向誘導(dǎo)因子b 的表達(dá)式變?yōu)?
當(dāng)軸向誘導(dǎo)因子a 較大時(shí),根據(jù)動(dòng)量定理,葉輪后方的尾流速度將很低,甚至發(fā)生反向。實(shí)際上這種情況不可能發(fā)生,取而代之發(fā)生的是尾流變成了湍流,而動(dòng)量定理不再適用。同時(shí)通過混合過程,湍流尾流從周圍流體中重新獲得能量。針對這種現(xiàn)象,Glauert 提出了一種利用經(jīng)驗(yàn)公式對a 進(jìn)行修正的方法:
利用BEM 理論進(jìn)行海流機(jī)水動(dòng)力性能計(jì)算的步驟如下:
1)確定計(jì)算工況:流速U1、尖速比λ0、槳距角β0。
2)選取計(jì)算截面。該截面翼型弦長為c(r),扭角為β(r)。
3)對截面軸向誘導(dǎo)因子a 和周向誘導(dǎo)因子b 賦初值??闪頰0= 0,b0= 0。
4)利用式(5)計(jì)算確定φ,則攻角α = φ - β(r)- β0。
5)利用攻角α 查得升力系數(shù)CL和阻力系數(shù)CD。利用式(16)和(17)計(jì)算損失因子F(r)。利用式(18)對a0和b0進(jìn)行修正,修正后的軸向誘導(dǎo)因子和周向誘導(dǎo)因子記為a 和b。
6)判斷a 是否大于0.38。若a 大于0.38,利用式(19)對a 進(jìn)行修正。若不滿足,則直接進(jìn)入步驟g。
7)計(jì)算Δa= |a-a0|和Δb= |b-b0|。判斷Δa 和Δb 是否小于收斂精度ε:若滿足,進(jìn)入步驟h;否則令a0=a,b0=b,并重復(fù)步驟c 至步驟g,直至收斂。
8)利用式(4)和式(5)計(jì)算U 和φ,然后利用式(7)和式(8)計(jì)算截面dT/dr 和dM/dr。
9)判斷是否所有計(jì)算截面計(jì)算完畢:若滿足,進(jìn)入步驟10;否則,選取新的計(jì)算截面,返回步驟2。
10)對各個(gè)截面求得的dF 和dM 進(jìn)行積分,求得軸向推力T 和轉(zhuǎn)矩M,捕獲的功率P =MΩ。為了便于比較不同尺寸海流機(jī)的性能,對載荷與功率進(jìn)行無因次化,求得推力系數(shù)CT、轉(zhuǎn)矩系數(shù)CM和獲能系數(shù)CP。
利用BEM 理論,對所設(shè)計(jì)的150 kW 水平軸海流機(jī)在定槳距運(yùn)行狀態(tài)和變槳距運(yùn)行狀態(tài)下的水動(dòng)力性能分別進(jìn)行了計(jì)算。
3.4.1 固定槳距
設(shè)海流機(jī)葉片槳距角β0固定為0°,來流U1保持為2 m/s 不變,計(jì)算不同尖速比λ0下海流機(jī)的水動(dòng)力性能。圖5 和圖6 所示的是不同尖速比λ0下,α、dCT/dr 和dCM/dr 沿葉片徑向的分布情況。其中dCT/dr 和dCM/dr 為無因次化的單位長度葉片推力和轉(zhuǎn)矩,定義為:
圖5 不同尖速比下攻角徑向分布Fig.5 Radial distribution of attack angle under different tip speed ratios
由圖5 可知,隨著尖速比λ0的增大,葉片的攻角隨之減小。當(dāng)λ0等于2 或3 時(shí),整個(gè)葉片的攻角均大于失速角,葉片工作于失速狀態(tài)中。當(dāng)λ0等于6 時(shí),即處于設(shè)計(jì)尖速比時(shí),整個(gè)葉片的攻角基本處于6°到8°之間,接近于翼型最大升阻比對應(yīng)的最佳攻角。此外,由于葉尖損失的影響,導(dǎo)致葉梢處的攻角明顯減小。
由圖6 可知,隨著尖速比λ0的增大,葉片產(chǎn)生的軸向推力也增大,最大轉(zhuǎn)矩則出現(xiàn)在尖速比λ0等于4附近。當(dāng)λ0等于2 或3 時(shí),推力分布曲線與轉(zhuǎn)矩分布曲線的變化趨勢不同于λ0較大時(shí)的情況,這是由葉片失速所導(dǎo)致。當(dāng)λ0增大時(shí),葉尖處首先離開失速狀態(tài),因此葉片載荷分布也由葉尖處開始逐漸恢復(fù)正常。對于正常工作的葉片,在同一尖速比λ0下,最大推力出現(xiàn)在r =5.00 m 附近,最大轉(zhuǎn)矩出現(xiàn)在r =4.75 m附近。
圖6 不同尖速比下推力和轉(zhuǎn)矩徑向分布Fig.6 Radial distribution of thrust force and torque under different tip speed ratios
圖7 CT、CM 和CP 隨尖速比λ0 變化曲線Fig.7 Curve of CT,CM and CP to λ0
圖7 所示的是水平軸海流機(jī)推力系數(shù)CT、轉(zhuǎn)矩系數(shù)CM和獲能系數(shù)CP隨著尖速比λ0變化的曲線。如前文所述,隨著尖速比λ0增大,推力系數(shù)CT增大。轉(zhuǎn)矩系數(shù)CM在λ0=4 時(shí)取得最大值CMmax=0.092。而獲能系數(shù)CP則在設(shè)計(jì)尖速比λ0=6 時(shí)取得最大值CPmax=0.420,且在λ0=6 附近CP數(shù)值變化幅度相對較小。說明利用BEM 理論進(jìn)行水平軸海流機(jī)葉片的設(shè)計(jì)基本能實(shí)現(xiàn)預(yù)期要求。
3.4.2 可變槳距
圖8 不同槳距角下CP 隨尖速比λ0 變化曲線Fig.8 Curves of CP to λ0 with different pitch angles
葉片設(shè)計(jì)完畢之后,葉片截面的各項(xiàng)參數(shù)也隨之確定,因此在相同的工況下,定槳距海流機(jī)的水動(dòng)力性能也是確定的。對于變槳距海流機(jī),則可以通過調(diào)節(jié)葉片的槳距角β0改變來流攻角以調(diào)節(jié)葉片的水動(dòng)力性能,從而實(shí)現(xiàn)過載功率控制、提高啟動(dòng)轉(zhuǎn)矩等要求[8]。
圖8 所示的是槳距角β0為0°、±5°和±10°時(shí),獲能系數(shù)CP隨尖速比λ0變化的曲線??梢园l(fā)現(xiàn),當(dāng)葉片運(yùn)行在低尖速比條件下時(shí),正的槳距角能提高裝置的獲能系數(shù)CP,而負(fù)的槳距角將減小裝置的獲能系數(shù)CP。而當(dāng)葉片運(yùn)行在高尖速比條件下時(shí),無論正的槳距角還是負(fù)的槳距角,都將減小裝置的獲能系數(shù)CP。槳距角越大,對獲能系數(shù)CP的影響越是明顯。另外,正的槳距角下獲能系數(shù)曲線較為平坦,負(fù)的槳距角下獲能系數(shù)曲線則比較陡峭,說明運(yùn)行在負(fù)槳距角下的葉片對尖速比λ0的變化較為敏感。
除了使用相關(guān)理論計(jì)算海流機(jī)水動(dòng)力性能之外,利用CFD 方法對海流機(jī)進(jìn)行數(shù)值計(jì)算也是驗(yàn)證設(shè)計(jì)結(jié)果的重要方法。本文使用Fluent 軟件,利用多重參考系(MRF)模型[9]模擬了海流機(jī)在均勻來流U1、固定槳距β0、不同尖速比λ0下的工作情況。
在進(jìn)行數(shù)值計(jì)算前,首先使用Gambit 軟件建立數(shù)值計(jì)算模型,劃分計(jì)算網(wǎng)格。
計(jì)算模型如圖9 所示。整個(gè)流場劃分為兩個(gè)相嵌套的圓柱形流體區(qū)域A 和B。其中區(qū)域A 是固定區(qū)域,長150 m,半徑50 m,左側(cè)底面的邊界條件設(shè)置為速度入口(velocity inlet)類型,右側(cè)底面和側(cè)面邊界條件設(shè)置為出流(outflow)類型。區(qū)域B 是旋轉(zhuǎn)區(qū)域,長10 m,半徑10 m,兩個(gè)區(qū)域的接觸面設(shè)置為交界面(interface)類型。海流機(jī)葉片位于較小的圓柱形流體區(qū)域B 中,設(shè)置為壁面(wall)類型,并將隨著區(qū)域B 一起轉(zhuǎn)動(dòng)。
圖9 水平軸海流機(jī)水動(dòng)力性能數(shù)值計(jì)算模型Fig.9 Numerical model of horizontal axis current turbine
本文計(jì)算模擬了150 kW 水平軸海流機(jī)在固定來流U1=2 m/s,不同尖速比λ0下的工作情況。
圖10 所示的是額定工況(U1=2 m/s,λ0=6)下海流機(jī)葉片的迎流面和背流面的壓力分布等值線。可以發(fā)現(xiàn),葉片迎流面為高壓區(qū),背流面為低壓區(qū),且這種壓力差在葉尖處尤為明顯。
圖11 所示的是額定工況(U1=2 m/s,λ0=6)下葉輪(y=0 m)近后方截面(y= + 0.2 m)的流體合成速度等值線云圖。從徑向來看,隨著半徑的增大,流體的合成速度首先逐漸減小,說明葉根處吸收的流體動(dòng)能較少,葉尖處吸收的流體動(dòng)能較大。當(dāng)半徑增大到靠近葉稍處時(shí),由于葉尖損失現(xiàn)象,流體動(dòng)能的吸收率相對減小,故流體合成速度有所增大。
圖12 比較了利用BEM 方法與CFD 方法計(jì)算獲得的獲能系數(shù)曲線。可以發(fā)現(xiàn),利用CFD 方法計(jì)算獲得的獲能系數(shù)曲線峰值同樣位于設(shè)計(jì)尖速比λ0=6 時(shí),再次驗(yàn)證了水平軸海流機(jī)葉片設(shè)計(jì)方法的有效性。通過CFD 方法計(jì)算求得的獲能系數(shù)最大值CPmax=0.372,為BEM 方法計(jì)算所得值的88.6%。
圖10 額定工況下海流機(jī)葉片壓力分布等值線Fig.10 Pressure contour of current turbine under rated condition
圖11 額定工況葉輪近后方速度等值線云圖Fig.11 Velocity contour behind the turbine under rated condition
圖12 BEM 理論與CFD 方法CP 曲線比較Fig.12 Comparison of CP curves by BEM and CFD methods
葉片的獲能系數(shù)CP決定了水平軸海流機(jī)的發(fā)電效率,葉片設(shè)計(jì)的目標(biāo)就是使其獲能系數(shù)CP盡可能大。BEM 理論是用于海流機(jī)葉片水動(dòng)力外型設(shè)計(jì)的常用方法。通過使用BEM 理論結(jié)合普朗特修正和葛勞渥修正的理論計(jì)算方法以及CFD 方法,分別對150 kW 水平軸海流機(jī)的水動(dòng)力性能進(jìn)行了計(jì)算,兩種方法的結(jié)果都表明該海流機(jī)的最大獲能系數(shù)CPmax位于設(shè)計(jì)尖速比處,說明本文的葉片設(shè)計(jì)方法是有效的。
通過對以上兩種方法的計(jì)算結(jié)果進(jìn)行分析,發(fā)現(xiàn)葉片從海流中截獲的能量和受到的載荷主要集中在葉尖區(qū)域,但是由于葉尖損失的原因,葉梢處吸收的能量和受到的載荷有明顯的下降。另外,槳距角對海流機(jī)水動(dòng)力性能的影響明顯,通過調(diào)節(jié)槳距角可以實(shí)現(xiàn)海流機(jī)的功率和載荷控制。
[1]Blunden L S,Bahaj A S.Tidal energy resource assessment for tidal stream generators[J].Journal of Power and Energy,2007,221(2):137-146.
[2]劉美琴,仲 穎,鄭 源,等.海流能利用技術(shù)研究進(jìn)展與展望[J].可再生能源,2009,27(5):78-81.
[3]Ponta F L,Jacovkis P M.Marine-current power generation by diffuser-augmented floating hydro-turbines[J].Renewable Energy,2008,33 (4):665-673.
[4]李 偉,林勇剛,劉宏偉,等.水平軸螺旋槳式海流能發(fā)電技術(shù)研究[C]∥中國可再生能源學(xué)會(huì)海洋能專業(yè)委員會(huì)第一屆學(xué)術(shù)討論會(huì)文集.2008:81-90.
[5]Burton T.風(fēng)能技術(shù)[M].北京:科學(xué)出版社.2007:37-55.
[6]馬 舜,李 偉,劉宏偉,等.水平軸潮流能發(fā)電系統(tǒng)能量捕獲機(jī)構(gòu)研究[J].機(jī)械工程學(xué)報(bào),2010,46(18):150-156.
[7]王 立,董志勇,韓 偉,等.海流能轉(zhuǎn)換器葉片翼型的失速和水動(dòng)力特性的數(shù)值研究[J].水動(dòng)力學(xué)研究與進(jìn)展:A輯,2011,26(1):58-64.
[8]馬 舜,李 偉,劉宏偉,等.海流能發(fā)電系統(tǒng)的最大功率跟蹤控制研究[J].太陽能學(xué)報(bào),2011,32(4):577-582.
[9]Wang J,Muller N.Performance prediction of array arrangement on ducted composite material marine current turbines (CMMCTs)[J].Ocean Engineering,2012,41:21-26.