何煜平，楊建民

(上海交通大學 海洋工程國家重點實驗室，上海 200240)

所謂海流，是風應力、密度梯度、引潮力等作用而形成的大規(guī)模相對穩(wěn)定的流動。據統(tǒng)計，全球海流能蘊藏總量接近于3.5TW［1-2］。若能合理地將這些能源加以利用，則能為緩解能源危機提供一條新的道路。

海流能發(fā)電的方式類似于風能發(fā)電，從概念上來說都是利用流體產生的升力或者阻力使獲能裝置運動，并帶動內部的發(fā)電機發(fā)電。與風能相比，海流能發(fā)電具有以下的優(yōu)勢:1)海水密度是空氣密度的830 倍，使得海流能具有比風能更大的能量密度;2)海流，特別是潮流，其流向和流量相對穩(wěn)定并且容易預測［3-4］。

目前，海流能技術最為先進、產業(yè)化進程最為順利的國家為英國。2008年4月，英國MCT 公司完成了1.2 MW 的水平軸海流機發(fā)電裝置“SeaGen”的安裝和運營，標志著世界上第一個商業(yè)化規(guī)模的海流發(fā)電系統(tǒng)投入使用［2］。國內的部分學校和研究所也都在積極開展海流能發(fā)電的研究，但尚處于起步階段。

目前海流能領域的研究熱點主要集中在以下方面:1)海流機葉片的水動力設計方法;2)葉片可調距機構的設計;3)安裝與支撐方式的設計研究;4)葉片氣蝕的研究和抑制;5)機組布局的研究;6)海流能發(fā)電裝置新概念的設計［2-4］。

本文利用了葉素－動量理論［5-6］，對給定額定工況的海流能發(fā)電機的葉片進行了設計，并分別利用理論方法和CFD 方法對葉片的水動性能進行了計算和分析。

1 葉素－動量(BEM)理論

葉素－動量理論是同時利用動量定理和葉素理論進行葉片水動力外型設計的方法。若葉輪遠前方流速為U1，遠后方流速為U2，葉輪轉速為Ω，旋轉尾流角速度為ω，可定義軸向誘導因子a 和周向誘導因子b:

由動量定理，葉輪受到的軸向力T 和轉矩M 等于單位時間內通過葉輪的流體動量和角動量的變化量。對于葉輪上一個半徑為r、寬度為dr 的圓環(huán)，其受到的軸向力dT 和轉矩dM 可表示為:

葉素理論則是將葉片沿徑向分成若干葉素，對每個葉素分別進行分析，然后將作用在每個葉素上的力疊加，得到整個葉輪上的流體作用力。位于半徑r 處葉素的速度和受力情況如圖1 所示。

葉素上的流體相對速度U 為:

葉輪旋轉面與流體相對速度之間的夾角φ 為:

式中:α 為攻角;β 為扭角;λ 為局部速度比，λ =Ωr/U1，其中葉尖處局部速度比稱尖速比，記為λ0。

設在攻角α 下，翼型的升力系數和阻力系數分別為CL和CD。則對于弦長為c，寬度為dr 的葉素，其受到的升力dL 和阻力dD 可分別表示為:

若葉片數為N，則半徑為r、寬度為dr 的葉輪圓環(huán)上受到的軸向力dT 和轉矩dM 可以表示為:

令式(2)和式(7)相等，式(3)和式(8)相等，可得

Tony Burton 通過比較考慮阻力和忽略阻力兩種情形下的最優(yōu)葉片設計，發(fā)現阻力對葉片的最優(yōu)設計影響很?。?］。因此，忽略葉片阻力，則式(9)變?yōu)?

根據式(5)和式(10)，得到a、b 和λ 三者之間的關系:

圖1 葉素的速度和作用力分析Fig.1 Velocities and forces of a blade element

2 基于BEM 理論的海流機葉片設計

2.1 基本設計步驟

表征海流機葉片性能的參數主要包括推力系數CT、轉矩系數CM和獲能系數CP:

海流機葉片水動力外形設計，是在確定的設計流速U1、尖速比λ0和輸出功率P0的條件下，使葉片的獲能系數CP盡可能大。設計時，需要首先確定葉片數目N、葉輪半徑R、葉轂半徑RH，并選定合適的翼型。目前比較常見的葉片數主要有2 葉和3 葉。葉片數越多，獲能系數CP越大，但軸向力系數CT也隨之大幅提高，對支撐系統(tǒng)的設計有較高的要求。葉輪半徑R 可以通過下式估算:

式中的獲能系數CP是通過經驗確定的預估值。葉轂半徑RH可根據設計經驗，取葉輪半徑合適的百分數。葉片翼型的選取可以參照各翼型的水動力性能曲線，選取最大升阻比L/Dmax較大的翼型，同時還需保證具有較好的失速特性。通常葉尖處翼型較薄，葉根處翼型較厚，并逐漸過渡為輪轂處的圓形截面。

確定了上述設計參數后，可利用BEM 理論確定弦長c 和扭角β 沿葉輪徑向的分布。步驟如下:

1)選取一控制截面，其徑向位置為r。計算局部速度比λ=Ωr / U1。

2)確定a 和b，使獲能系數CP最大。根據式(8)，葉輪圓環(huán)處的功率可表示為:

因此，(1 －a)b 的取值越大，dP 就越大。但是，a、b 和λ 三者之間受到了式(11)的約束。為了確定a 和b 的最優(yōu)值，需要求解以下單目標帶約束的非線性優(yōu)化問題:

3)根據控制截面處翼型升力系數曲線，找到對應最大升阻比的最優(yōu)攻角α 以及相應的升力系數CL。

4)根據式(5)確定φ，則扭角β=φ－α。根據式(9)確定弦長c。

5)重復步驟1)～4)，直到所有控制截面計算完畢。

2.2 翼型水動力性能

在葉片設計時，需要葉片截面翼型的升力系數曲線和阻力系數曲線作為輸入參數。對于NACA 翼型族相關性能曲線，可通過XFoil 等翼型分析軟件獲得，也可通過CFD 方法求得［7］。本文所使用的翼型性能曲線是通過Fluent 軟件求得，其中NACA 2414、NACA 2416 和NACA 2418 翼型的升力系數曲線和阻力系數曲線如圖2 所示。翼型失速角約為16°到18°，最大升阻比對應的最佳攻角約為7°。

圖2 NACA 24 系列翼型升力系數和阻力系數曲線Fig.2 Lift coefficient and drag coefficient curves of NACA 24 series

2.3 150 kW 水平軸海流機葉片設計實例

設有如表1 所示的設計要求，利用BEM 理論進行水平軸海流機葉片的設計。設計時，為了簡化問題，認為海流機所處水深較深，海流不受自由表面波浪的影響，且認為海流流速在整個葉輪處均勻分布。由于海流相對比較穩(wěn)定，在通常情形下流速能保持在較小范圍內浮動，因此設計流速可以從該流速范圍內選取。

表1 海流機葉片設計參數Tab.1 Parameters of blade designing of a current turbine

取獲能系數預估值為0.4，根據式(13)計算葉片半徑，得R≈5.5 m。葉片翼型選取NACA 24 系列，葉根處為了保證足夠的結構強度，選取葉片厚度t 為18%，然后向葉尖處逐漸過渡為16%和14%。利用前述設計步驟計算截面弦長c 和扭角β，c 和β 沿著葉輪徑向的分布曲線如圖3 所示。利用Pro/E 進行葉片實體建模，如圖4 所示。

圖3 150 kW 海流機葉片幾何參數徑向分布Fig.3 Radial distribution of geometric parameters of the 150 kW current turbine's blade

圖4 150 kW 海流機葉片實體模型Fig.4 Model of the 150 kW current turbine's blade

3 海流機水動力性能計算

為了檢驗葉片設計是否合理，以及預測不同工況下的性能，需對海流機的水動力性能進行計算。BEM理論結合普朗特(Prandtl)修正以及葛勞渥(Glauert)修正是較常用的海流機水動力性能理論計算方法［5-6］。

3.1 Prandtl 葉尖/輪轂損失修正

葉素-動量理論忽略了順著翼展方向的速度分量，然而實際上水流在葉片上存在著徑向流動。尤其在葉尖和輪轂處，這種三維流動效應更是明顯，因此葉尖損失和輪轂損失對葉片性能的影響不容忽略。Prandtl針對該現象提出了Prandtl 漸進法。根據該理論，葉尖損失和輪轂損失可用損失因子Ft(r)和Fh(r)來表示:

總的損失因子F (r)可表示為:

經過修正后的軸向誘導因子a 和軸向誘導因子b 的表達式變?yōu)?

3.2 Glauert 軸向誘導因子修正

當軸向誘導因子a 較大時，根據動量定理，葉輪后方的尾流速度將很低，甚至發(fā)生反向。實際上這種情況不可能發(fā)生，取而代之發(fā)生的是尾流變成了湍流，而動量定理不再適用。同時通過混合過程，湍流尾流從周圍流體中重新獲得能量。針對這種現象，Glauert 提出了一種利用經驗公式對a 進行修正的方法:

3.3 海流機水動力性能理論計算步驟

利用BEM 理論進行海流機水動力性能計算的步驟如下:

1)確定計算工況:流速U1、尖速比λ0、槳距角β0。

2)選取計算截面。該截面翼型弦長為c(r)，扭角為β(r)。

3)對截面軸向誘導因子a 和周向誘導因子b 賦初值?？闪頰0= 0，b0= 0。

4)利用式(5)計算確定φ，則攻角α = φ － β(r)－ β0。

5)利用攻角α 查得升力系數CL和阻力系數CD。利用式(16)和(17)計算損失因子F(r)。利用式(18)對a0和b0進行修正，修正后的軸向誘導因子和周向誘導因子記為a 和b。

6)判斷a 是否大于0.38。若a 大于0.38，利用式(19)對a 進行修正。若不滿足，則直接進入步驟g。

7)計算Δa= |a－a0|和Δb= |b－b0|。判斷Δa 和Δb 是否小于收斂精度ε:若滿足，進入步驟h;否則令a0=a，b0=b，并重復步驟c 至步驟g，直至收斂。

8)利用式(4)和式(5)計算U 和φ，然后利用式(7)和式(8)計算截面dT/dr 和dM/dr。

9)判斷是否所有計算截面計算完畢:若滿足，進入步驟10;否則，選取新的計算截面，返回步驟2。

10)對各個截面求得的dF 和dM 進行積分，求得軸向推力T 和轉矩M，捕獲的功率P =MΩ。為了便于比較不同尺寸海流機的性能，對載荷與功率進行無因次化，求得推力系數CT、轉矩系數CM和獲能系數CP。

3.4 海流機水動力性能計算結果與分析

利用BEM 理論，對所設計的150 kW 水平軸海流機在定槳距運行狀態(tài)和變槳距運行狀態(tài)下的水動力性能分別進行了計算。

3.4.1 固定槳距

設海流機葉片槳距角β0固定為0°，來流U1保持為2 m/s 不變，計算不同尖速比λ0下海流機的水動力性能。圖5 和圖6 所示的是不同尖速比λ0下，α、dCT/dr 和dCM/dr 沿葉片徑向的分布情況。其中dCT/dr 和dCM/dr 為無因次化的單位長度葉片推力和轉矩，定義為:

圖5 不同尖速比下攻角徑向分布Fig.5 Radial distribution of attack angle under different tip speed ratios

由圖5 可知，隨著尖速比λ0的增大，葉片的攻角隨之減小。當λ0等于2 或3 時，整個葉片的攻角均大于失速角，葉片工作于失速狀態(tài)中。當λ0等于6 時，即處于設計尖速比時，整個葉片的攻角基本處于6°到8°之間，接近于翼型最大升阻比對應的最佳攻角。此外，由于葉尖損失的影響，導致葉梢處的攻角明顯減小。

由圖6 可知，隨著尖速比λ0的增大，葉片產生的軸向推力也增大，最大轉矩則出現在尖速比λ0等于4附近。當λ0等于2 或3 時，推力分布曲線與轉矩分布曲線的變化趨勢不同于λ0較大時的情況，這是由葉片失速所導致。當λ0增大時，葉尖處首先離開失速狀態(tài)，因此葉片載荷分布也由葉尖處開始逐漸恢復正常。對于正常工作的葉片，在同一尖速比λ0下，最大推力出現在r =5.00 m 附近，最大轉矩出現在r =4.75 m附近。

圖6 不同尖速比下推力和轉矩徑向分布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 所示的是水平軸海流機推力系數CT、轉矩系數CM和獲能系數CP隨著尖速比λ0變化的曲線。如前文所述，隨著尖速比λ0增大，推力系數CT增大。轉矩系數CM在λ0=4 時取得最大值CMmax=0.092。而獲能系數CP則在設計尖速比λ0=6 時取得最大值CPmax=0.420，且在λ0=6 附近CP數值變化幅度相對較小。說明利用BEM 理論進行水平軸海流機葉片的設計基本能實現預期要求。

3.4.2 可變槳距

圖8 不同槳距角下CP 隨尖速比λ0 變化曲線Fig.8 Curves of CP to λ0 with different pitch angles

葉片設計完畢之后，葉片截面的各項參數也隨之確定，因此在相同的工況下，定槳距海流機的水動力性能也是確定的。對于變槳距海流機，則可以通過調節(jié)葉片的槳距角β0改變來流攻角以調節(jié)葉片的水動力性能，從而實現過載功率控制、提高啟動轉矩等要求［8］。

圖8 所示的是槳距角β0為0°、±5°和±10°時，獲能系數CP隨尖速比λ0變化的曲線。可以發(fā)現，當葉片運行在低尖速比條件下時，正的槳距角能提高裝置的獲能系數CP，而負的槳距角將減小裝置的獲能系數CP。而當葉片運行在高尖速比條件下時，無論正的槳距角還是負的槳距角，都將減小裝置的獲能系數CP。槳距角越大，對獲能系數CP的影響越是明顯。另外，正的槳距角下獲能系數曲線較為平坦，負的槳距角下獲能系數曲線則比較陡峭，說明運行在負槳距角下的葉片對尖速比λ0的變化較為敏感。

4 基于Fluent 的水平軸海流機性能計算

除了使用相關理論計算海流機水動力性能之外，利用CFD 方法對海流機進行數值計算也是驗證設計結果的重要方法。本文使用Fluent 軟件，利用多重參考系(MRF)模型［9］模擬了海流機在均勻來流U1、固定槳距β0、不同尖速比λ0下的工作情況。

4.1 計算模型

在進行數值計算前，首先使用Gambit 軟件建立數值計算模型，劃分計算網格。

計算模型如圖9 所示。整個流場劃分為兩個相嵌套的圓柱形流體區(qū)域A 和B。其中區(qū)域A 是固定區(qū)域，長150 m，半徑50 m，左側底面的邊界條件設置為速度入口(velocity inlet)類型，右側底面和側面邊界條件設置為出流(outflow)類型。區(qū)域B 是旋轉區(qū)域，長10 m，半徑10 m，兩個區(qū)域的接觸面設置為交界面(interface)類型。海流機葉片位于較小的圓柱形流體區(qū)域B 中，設置為壁面(wall)類型，并將隨著區(qū)域B 一起轉動。

圖9 水平軸海流機水動力性能數值計算模型Fig.9 Numerical model of horizontal axis current turbine

4.2 數值計算結果與分析

本文計算模擬了150 kW 水平軸海流機在固定來流U1=2 m/s，不同尖速比λ0下的工作情況。

圖10 所示的是額定工況(U1=2 m/s，λ0=6)下海流機葉片的迎流面和背流面的壓力分布等值線?？梢园l(fā)現，葉片迎流面為高壓區(qū)，背流面為低壓區(qū)，且這種壓力差在葉尖處尤為明顯。

圖11 所示的是額定工況(U1=2 m/s，λ0=6)下葉輪(y=0 m)近后方截面(y= + 0.2 m)的流體合成速度等值線云圖。從徑向來看，隨著半徑的增大，流體的合成速度首先逐漸減小，說明葉根處吸收的流體動能較少，葉尖處吸收的流體動能較大。當半徑增大到靠近葉稍處時，由于葉尖損失現象，流體動能的吸收率相對減小，故流體合成速度有所增大。

圖12 比較了利用BEM 方法與CFD 方法計算獲得的獲能系數曲線?？梢园l(fā)現，利用CFD 方法計算獲得的獲能系數曲線峰值同樣位于設計尖速比λ0=6 時，再次驗證了水平軸海流機葉片設計方法的有效性。通過CFD 方法計算求得的獲能系數最大值CPmax=0.372，為BEM 方法計算所得值的88.6%。

圖10 額定工況下海流機葉片壓力分布等值線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

5 結 語

葉片的獲能系數CP決定了水平軸海流機的發(fā)電效率，葉片設計的目標就是使其獲能系數CP盡可能大。BEM 理論是用于海流機葉片水動力外型設計的常用方法。通過使用BEM 理論結合普朗特修正和葛勞渥修正的理論計算方法以及CFD 方法，分別對150 kW 水平軸海流機的水動力性能進行了計算，兩種方法的結果都表明該海流機的最大獲能系數CPmax位于設計尖速比處，說明本文的葉片設計方法是有效的。

通過對以上兩種方法的計算結果進行分析，發(fā)現葉片從海流中截獲的能量和受到的載荷主要集中在葉尖區(qū)域，但是由于葉尖損失的原因，葉梢處吸收的能量和受到的載荷有明顯的下降。另外，槳距角對海流機水動力性能的影響明顯，通過調節(jié)槳距角可以實現海流機的功率和載荷控制。

［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］劉美琴，仲 穎，鄭 源，等.海流能利用技術研究進展與展望［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ā)電技術研究［C］∥中國可再生能源學會海洋能專業(yè)委員會第一屆學術討論會文集.2008:81-90.

［5］Burton T.風能技術［M］.北京:科學出版社.2007:37-55.

［6］馬 舜，李 偉，劉宏偉，等.水平軸潮流能發(fā)電系統(tǒng)能量捕獲機構研究［J］.機械工程學報，2010，46(18):150-156.

［7］王 立，董志勇，韓 偉，等.海流能轉換器葉片翼型的失速和水動力特性的數值研究［J］.水動力學研究與進展:A輯，2011，26(1):58-64.

［8］馬 舜，李 偉，劉宏偉，等.海流能發(fā)電系統(tǒng)的最大功率跟蹤控制研究［J］.太陽能學報，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.