水平軸潮流能葉輪設(shè)計(jì)與水動(dòng)力特性分析

張之陽，王曉航，劉葳興，紀(jì)仁瑋，郭廣廓

(1.江蘇海洋大學(xué) 機(jī)械與海洋工程學(xué)院，江蘇 連云港 222005；2.哈爾濱 大電機(jī)研究所，黑龍江 哈爾濱 150001；3.哈爾濱工程大學(xué) 海洋可再生能源研究所，黑龍江 哈爾濱 150001；4.工業(yè)和信息化部電子第五研究所，廣東 廣州 510000)

0 引 言

隨著能源問題日益嚴(yán)峻，作為一種清潔海洋可再生能源，潮流能具有儲(chǔ)量豐富、分布集中且可預(yù)測性強(qiáng)等優(yōu)點(diǎn)。如何高效開發(fā)和利用潮流能已受到國內(nèi)外學(xué)者的關(guān)注[1]。水平軸葉輪是一種常見的潮流能發(fā)電裝置，具有效率較高，功率波動(dòng)較小，自啟動(dòng)性能良好等特點(diǎn)[2]。目前，已經(jīng)在海上部署了一些商用前的水平軸潮流能葉輪原型，典型代表有：英國MCT公司的SeaFlow系列，英國TGL公司的Alstom，挪威Hammerfest公司的HS系列，新加坡Altantis公司的AR，AK系列[3]。哈爾濱工程大學(xué)設(shè)計(jì)的水平軸葉輪10 kW的“海明 I號”，2×100 kW 的“海能 II號”[4]。

葉片是水平軸潮流能裝置一級能量轉(zhuǎn)換的核心部件，葉片型線設(shè)計(jì)對能量轉(zhuǎn)換效率和運(yùn)行穩(wěn)定性至關(guān)重要，直接影響發(fā)電效率。葉片的性能主要取決選取的翼型形狀以及沿展長方向的形狀變化。前者根據(jù)翼型氣動(dòng)力性能選擇，后者主要由沿展長方向各葉素截面的弦長和槳距角確定[5-6]。水平軸葉輪的葉片設(shè)計(jì)方法主要有：基于圓盤理論的簡化風(fēng)車模型，基于渦流理論的Schmitz模型、Glauert模型和Wilson模型。其中Glauert模型考慮了軸向、切向誘導(dǎo)因子，設(shè)計(jì)結(jié)果具有較高的精度[7]。

本文基于哈爾濱工程大學(xué)承擔(dān)的國家公益性項(xiàng)目要求，完成水平軸潮流能葉輪的設(shè)計(jì)，研究并發(fā)展一套水平軸潮流能葉輪水動(dòng)力性能有效的預(yù)報(bào)方法，為今后水平軸潮流能葉輪的性能設(shè)計(jì)建立可靠的理論分析方法，積累準(zhǔn)確的實(shí)驗(yàn)數(shù)據(jù)。針對項(xiàng)目要求的技術(shù)參數(shù)與環(huán)境參數(shù)，應(yīng)用Glauert渦流設(shè)計(jì)理論，對水平軸定槳距葉輪的結(jié)構(gòu)形式進(jìn)行設(shè)計(jì)，并采用BEM葉素動(dòng)量理論和CFD數(shù)值模擬方法對設(shè)計(jì)的水輪機(jī)進(jìn)行載荷與性能的預(yù)報(bào)?？偨Y(jié)出系統(tǒng)的設(shè)計(jì)方法與預(yù)報(bào)方法，歸納出水平軸葉輪的載荷特點(diǎn)。

1 葉片設(shè)計(jì)的 Glauert渦流設(shè)計(jì)理論

對于有限展長的葉片，其葉輪尾流中存在葉尖渦和葉根渦。美國Amherst大學(xué)改進(jìn)的Glauert渦流理論，考慮了葉輪引起的渦流影響，是目前應(yīng)用較為廣泛的理論之一[8]。

對于每個(gè)葉素來說，考慮渦流的影響，假設(shè)軸向和切向誘導(dǎo)速度子分別為a和b，則葉素的相對來流速度為：

葉素的入流角和槳距角可表示為：

根據(jù)動(dòng)量定理，作用在 dr段圓環(huán)處的推力為：

圖1 葉素受力示意圖Fig.1 The force diagram of blade element

根據(jù)葉素理論，可得：

其中：

聯(lián)立上式，并根據(jù)三角關(guān)系轉(zhuǎn)化可得：

假設(shè)每一個(gè)葉素均在理想狀態(tài)下運(yùn)行，定義Cd=0，即 ε =0°，可得：

其中： co tφ由下式確定， λ為半徑r處的尖速比，則有：

葉素圓環(huán)的功率和能量利用率表達(dá)為：

將式（14）代入上式，則有：

對Cp求導(dǎo)取最大值，可知a是 λ的函數(shù)，可表示為：

若葉輪半徑、葉輪轉(zhuǎn)速和來流速度給定，則 λ即為已知，軸向誘導(dǎo)因子可由上述函數(shù)關(guān)系求出，再代入式（14）即可求出切向誘導(dǎo)因子。根據(jù)軸向和切向誘導(dǎo)因子可得入流角，因此入流角也只與速比λ 有關(guān)，表示為：

將式（8）進(jìn)行無量綱處理，由于變量a與 φ均可以表示成 λ的相關(guān)函數(shù)，定義葉片形狀參數(shù)lc為：

上述推導(dǎo)表明形狀參數(shù)lc與 入流角 φ只與工作速比有關(guān)，與翼型的氣動(dòng)性能并無關(guān)系。在變槳距角，變弦長的葉片設(shè)計(jì)中，需要確定的參數(shù)即沿展長的弦長與槳距角分布情況。圖2和圖3是翼型形狀參數(shù)lc和入流角 φ與葉尖速比 λ的關(guān)系。

2 葉片型線確定

葉片弦長l、葉片截面翼型的攻角 α分別如下式：

圖2 葉片形狀參數(shù)與速比 λ的關(guān)系曲線Fig.2 Relationship curve between blade shape parameters lc and speed ratio λ

圖3 葉片入流角與速比 λ的關(guān)系曲線Fig.3 Relationship curve between blade inflow angle and speed ratio λ

式中：r為葉片不同位置的半徑，形狀參數(shù)lc可通過上圖查找；Cl為最優(yōu)升阻比下的升力系數(shù)；N為葉片數(shù)目； α0為升力系數(shù)為零時(shí)對應(yīng)的翼型攻角，一般情況下為負(fù)值；Rz為展弦比；Lm為 葉片平均弦長；Kl為翼型的升力曲線平均斜率；CLmax為升力曲線在失速前的最大值； αLmax為此時(shí)對應(yīng)的攻角。

葉片沿展長的型線通過葉素的弦長和槳距角確定，葉片設(shè)計(jì)的流程如下：

1）根據(jù)額定輸出功率確定葉輪的掃掠面積S；

2）確定葉輪直徑D；

3）根據(jù)水輪機(jī)工作速比為 λ ，確定葉輪轉(zhuǎn)速（主要用于設(shè)計(jì)軸承機(jī)械結(jié)構(gòu)以及匹配發(fā)電機(jī)）；

4）計(jì)算不同葉素的尖速比 λ*；

5）計(jì)算不同葉素的入流角 φ；

6）確定不同葉素的葉片形狀參數(shù)lc；

7）根據(jù)形狀參數(shù)lc確定葉片各葉素位置弦長l；

8）計(jì)算葉片的平均弦長Lm，升力曲線平均斜率Kl，葉片的展弦比Rz，葉片截面位置的翼型攻角α ；

9）確定葉片各葉素的槳距角；

10）驗(yàn)證設(shè)計(jì)葉片是否滿足設(shè)計(jì)要求。

采用NREL風(fēng)機(jī)標(biāo)準(zhǔn)翼型S809，該翼型最優(yōu)攻角為6.08°，根據(jù)上述流程計(jì)算葉素截面的弦長與槳距角，圖4為設(shè)計(jì)完成的三維葉片圖。

圖4 葉片三維示意圖Fig.4 Three dimensional diagram of turbine blade

3 水平軸葉輪水動(dòng)力性能分析方法

采用BEM葉素動(dòng)量理論與CFD數(shù)值模擬的2種方法，對本文設(shè)計(jì)的模型進(jìn)行載荷與性能的預(yù)報(bào)。

3.1 BEM葉素動(dòng)量理論方法

通過Matlab將BEM葉素動(dòng)量理論編程，即可用于葉輪載荷與性能的計(jì)算。BEM方法的求解流程如下[9]：

1）首先給出誘導(dǎo)速度因子a，b的初值（可設(shè)為0）；

2）根據(jù)公式計(jì)算每個(gè)葉素翼型的來流角度；

3）計(jì)算每個(gè)葉素翼型的攻角 α；

4）根據(jù)上一步得到的攻角 α ，找出其對應(yīng)的升力系數(shù)Cl與阻力系數(shù)Cd；

5）計(jì)算當(dāng)前葉素的推力系數(shù)CT與轉(zhuǎn)矩系數(shù)CM；

6）根據(jù)公式重新計(jì)算a和b的值；

7）返回步驟2，重新迭代，直至滿足容差要求。

3.2 CFD數(shù)值模擬方法

水平軸葉輪CFD數(shù)值模擬主要使用滑移網(wǎng)格方法，設(shè)置葉輪距離入口和兩側(cè)壁面均為3～4D（D為葉輪直徑），距離出口8～10D。旋轉(zhuǎn)域采用圓柱體，靜止域采用長方體或圓柱體均可。靜止域采用結(jié)構(gòu)化網(wǎng)格，旋轉(zhuǎn)域采用非結(jié)構(gòu)網(wǎng)格，網(wǎng)格效果如圖5所示。為充分模擬邊界層效應(yīng)，使得湍流模擬較為準(zhǔn)確，需保證葉片表面y+＜20[10]。

邊界條件的設(shè)置為：大氣壓為參考壓力，給定重力加速度的方向。入口邊界為速度入口，給定均勻來流速度、湍流參數(shù)。流體計(jì)算域的左右兩側(cè)和底面為自由滑動(dòng)壁面。流體計(jì)算域的出口和頂部為開放的壓力邊界，相對壓力設(shè)為0。葉片和輪轂表面設(shè)置為不可滑移壁面。給定旋轉(zhuǎn)域旋轉(zhuǎn)角速度，靜止域和旋轉(zhuǎn)域之間通過滑移交界面連接。計(jì)算中湍流模型采用SST模型，求解器為瞬態(tài)求解器，時(shí)間步長為葉輪旋轉(zhuǎn)3°所用的時(shí)間。

4 計(jì)算結(jié)果與分析

基于上述2種方法的計(jì)算，圖6對比了葉輪能量利用率Cp與葉尖速比λ 的關(guān)系。

可以看出，2種方法得到的能量利用率曲線隨速比的變化趨勢相同，均先增大后減小。 λ =4為最優(yōu)速比，峰值Cp約為40%。當(dāng) λ=10時(shí)，能量利用趨于零，此時(shí)葉輪處于空載狀態(tài)下的最高轉(zhuǎn)速，即飛逸轉(zhuǎn)速。

通過對比圖中的結(jié)果可以發(fā)現(xiàn)，BEM方法因忽略流體沿展向的流動(dòng)，以及粘性摩擦等，計(jì)算結(jié)果偏高。BEM方法和CFD方法的誤差在可接受的范圍內(nèi)，對于葉輪水動(dòng)力性能的預(yù)報(bào)均有較高的精度。因此考慮到時(shí)間成本，基于BEM方法繼續(xù)對此葉輪的載荷特性進(jìn)行研究，得到的葉輪的轉(zhuǎn)矩系數(shù)、軸向載荷系數(shù)隨速比的變化規(guī)律。

圖5 計(jì)算域網(wǎng)格示意圖Fig.5 Schematic diagram of computational domain mesh

圖6 水平軸潮流能水輪機(jī) C p-λ曲線Fig.6 The C p -λ curve of horizontal-axis tidal current turbine

從圖7可以看出，在低速比時(shí)，葉輪轉(zhuǎn)速系數(shù)較低，即葉輪啟動(dòng)時(shí)的主動(dòng)力矩較小。因此在設(shè)計(jì)軸系，水倉密封等時(shí)，需考慮軸系間的摩擦不易過大，否則會(huì)出現(xiàn)較難啟動(dòng)的問題。隨著轉(zhuǎn)速的增大，葉輪的主動(dòng)轉(zhuǎn)矩迅速增大，當(dāng) λ=3.5時(shí)達(dá)到最大。水輪機(jī)在最優(yōu)速比 λ =4時(shí)的轉(zhuǎn)矩并不是最大的，葉輪轉(zhuǎn)矩在沒有達(dá)到最優(yōu)速比時(shí)就已經(jīng)達(dá)到最大，提前了0.5個(gè)速比。當(dāng)速比繼續(xù)增大時(shí)，葉片的攻角隨即降低，并逐漸偏離最優(yōu)攻角，流體動(dòng)力性能下降，轉(zhuǎn)矩降低。葉片在速比10時(shí)，轉(zhuǎn)矩系數(shù)趨于0，能量利用趨于0，此時(shí)如果未加任何負(fù)載葉輪也不會(huì)再繼續(xù)做加速旋轉(zhuǎn)，即達(dá)到飛逸轉(zhuǎn)速。從圖8可以看出，軸向載荷系數(shù)隨著速比的增大而增大，當(dāng)超過最優(yōu)速比時(shí)，速比繼續(xù)增大，葉片處于失速狀態(tài)，軸向載荷系數(shù)增大的速度逐漸減緩，但仍然較大。

圖7 葉輪轉(zhuǎn)矩系數(shù)隨速比變化曲線Fig.7 Variation curve of turbine torque coefficient with speed ratio

圖8 葉輪軸向載荷系數(shù)隨速比變化曲線Fig.8 Variation curve of turbine axial load coefficient with speed ratio

為進(jìn)一步分析葉片表面的載荷分布情況，選取葉片沿展長方向3個(gè)位置進(jìn)行對比分析（30%，60%和90%葉展位置的葉素）。

從圖9可以看出，當(dāng)?shù)退俦葧r(shí)，葉輪根部葉片的攻角較大，處于失速狀態(tài)，效率較低。隨著速比的增大，葉片的轉(zhuǎn)速增大，葉片的攻角減小。為了使葉輪在工作速比時(shí)效率達(dá)到最優(yōu)，葉輪設(shè)計(jì)時(shí)應(yīng)盡量保證不同位置處的葉素，在工作速比時(shí)的攻角均達(dá)到最優(yōu)。在最優(yōu)速比 λ =4時(shí)，不同葉展位置處的攻角均達(dá)到該葉素翼型的最優(yōu)攻角6.08°左右，該計(jì)算結(jié)果證明本文葉輪設(shè)計(jì)方法的可靠性。

圖9 不同展長位置處攻角 α隨速比 λ變化曲線Fig.9 Variation curve of angle of attack α with speed ratio λ at different span positions

圖10 不同展長位置處軸向推力Ft隨速比 λ變化曲線Fig.10 Variation curve of axial thrust Ft with speed ratio λ at different span positions

圖11 不同展長位置處轉(zhuǎn)矩M隨速比 λ變化曲線Fig.11 Variation curve of torque M with speed ratio λ at different span positions

從圖10和圖11中的載荷曲線可以看出，不同葉素位置對于推力與轉(zhuǎn)速的貢獻(xiàn)是不同的，在靠近葉根的位置，由于轉(zhuǎn)速較低，且輪轂渦流系統(tǒng)所產(chǎn)生的渦旋尾流的影響，使得靠近輪轂處的葉片所受的軸向推力載荷較小，同時(shí)其動(dòng)力轉(zhuǎn)矩也較小。隨著葉素的半徑加大，軸向推力載荷與轉(zhuǎn)矩均有所變大，轉(zhuǎn)矩在速比 λ為3.5時(shí)達(dá)到最大。由于轉(zhuǎn)動(dòng)線速度隨著半徑的增大而變大，即半徑較大的葉素迎流速度較大，大半徑處葉素所受的流體動(dòng)力載荷也大于小半徑葉素。在最優(yōu)速比附近，提供葉片旋轉(zhuǎn)轉(zhuǎn)矩的部位主要分布在葉片展向上60%～90%處。但在低速比時(shí)，30%～60%小半徑處的葉素的轉(zhuǎn)矩貢獻(xiàn)更大，即小半徑處的葉素對于葉輪的啟動(dòng)性能起著至關(guān)重要的作用。因?yàn)樵谌~輪剛剛啟動(dòng)時(shí)大半徑處的葉素攻角很大，對啟動(dòng)轉(zhuǎn)矩貢獻(xiàn)較小，而小半徑處的葉素安裝角度較大，迎流攻角較小，因此對啟動(dòng)轉(zhuǎn)矩的貢獻(xiàn)更有流體動(dòng)力的優(yōu)勢。

5 結(jié) 語

本文基于Glauert渦流設(shè)計(jì)理論進(jìn)行水平軸潮流能葉輪的設(shè)計(jì)，并采用BEM葉素動(dòng)量理論與CFD數(shù)值模擬2種方法，對所設(shè)計(jì)的葉輪模型進(jìn)行載荷與性能的預(yù)報(bào)，證明此水輪機(jī)模型的工作性能達(dá)到了設(shè)計(jì)要求。由研究結(jié)果可知：

1）葉片的形狀參數(shù)和入流角只與工作速比有關(guān)，與翼型氣動(dòng)性能無關(guān)；

2）對于水平軸葉輪水動(dòng)力性能的預(yù)報(bào)，BEM方法和CFD方法均有較高的精度；

3）葉輪在最優(yōu)速比時(shí)轉(zhuǎn)矩并不是最大的，葉輪轉(zhuǎn)矩的最大值提前了0.5個(gè)速比；

4）工作速比下，提供葉片旋轉(zhuǎn)轉(zhuǎn)矩的部位主要分布在葉片展向上60%～90%處；

5）在低速比時(shí)，30%～60%小半徑處葉素的轉(zhuǎn)矩貢獻(xiàn)更大，即小半徑處的葉素對葉輪的自啟動(dòng)性能起著至關(guān)重要的作用。