水平軸潮流能水輪機設(shè)計與數(shù)值模擬

盛其虎，趙東亞，張亮

(哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江哈爾濱150001)

當(dāng)今社會能源短缺問題已變得日益嚴(yán)峻，作為一種清潔可再生能源，海洋潮流能具有儲量豐富、分布集中及可預(yù)測性強等特點。如何高效開發(fā)利用潮流能已受到國際社會的重視。

水輪機是一種常見的潮流能發(fā)電裝置，可分為垂直軸水輪機和水平軸水輪機。水平軸水輪機和水平軸風(fēng)力機原理相近，具有輸出功率穩(wěn)定的特點。1994年英國IT動力公司設(shè)計的風(fēng)車式發(fā)電裝置、2003年英國MCT公司研制的出首臺商業(yè)規(guī)模的潮流發(fā)電裝置均為水平軸形式的水輪機［1-2］。我國由哈爾濱工程大學(xué)設(shè)計的10 kW水平軸水輪機已于2012年成功投入運行。

葉片是潮流能水輪機能量轉(zhuǎn)換設(shè)備的關(guān)鍵部件，直接影響能量轉(zhuǎn)換的效率。葉片的性能主要取決于翼型形狀以及葉展形狀。前者根據(jù)翼型的氣動性能選擇，后者主要由沿展長方向各站面葉素的弦長和安裝角決定。風(fēng)力機葉片形狀的設(shè)計主要有基于圓盤理論的簡化設(shè)計模型，基于渦流模型的Schmitz設(shè)計模型、Glauert設(shè)計模型和Wilson設(shè)計模型。由于Wilson模型考慮了軸向、切向誘導(dǎo)系數(shù)，并且計入葉尖葉根損失，使得設(shè)計結(jié)果具有較高的精度［3-4］。本文選用Wilson模型進行水輪機葉片設(shè)計，并對設(shè)計的水輪機進行數(shù)值模擬以預(yù)測其水動力性能。

1 葉片設(shè)計的葉素動量理論

葉素理論的基本思想是:將葉片沿展長方向分成許多微段，這些微段稱為葉素，當(dāng)考慮軸向和徑向速度誘導(dǎo)因子a、b時，葉素周圍速度分布如圖1所示。

圖1 葉素周圍流速及受力Fig.1 Flow and forces on blade element

葉素在合速度為W的流體作用下，所受流體力可以分解為垂直于W的升力L和平行于W的阻力D，ω為葉輪旋轉(zhuǎn)角速度，V為來流速度，r為剖面半徑，且有:

式中:Cl、Cd分別為該葉素在攻角α下的升、阻力系數(shù)，c為葉素弦長。

葉素－動量定理(blade element momentum，BEM)的基本假定是:葉素所受力僅與通過葉素掃過的圓盤的流體的動量變化有關(guān)，即葉片受力等于流過流體的動量變化率［5-6］。由此可在軸向和切向建立如下等式:

定義葉尖速比λ:

式中:R為水輪機半徑。

對于不脫體的粘性繞流，阻力D是由于葉片表面的摩擦力產(chǎn)生的，在通過葉輪時并不影響壓降，因此在確定誘導(dǎo)因子時排除阻力因素Cd的影響。

考慮葉尖、葉根損失影響時需對BEM方程進行修正，本文采用Prandtl近似法確定的葉尖、葉根損失因數(shù)ftip、fhub分別為［7］

將總損失系數(shù)f=ftipfhub代入動量方程可得約束條件:

沿葉片展長方向取不同的截面，以產(chǎn)生的功率最大為設(shè)計優(yōu)化目標(biāo)。在各界面處，由單位葉素做功可得目標(biāo)函數(shù):

利用Matlab中的fmincon工具進行最優(yōu)化計算，便可得求各站面處的誘導(dǎo)因子a、b。進而可求得安裝角:

式中:α為升阻比最大時的攻角，與具體翼型有關(guān)。

將誘導(dǎo)因子代入式(3)或者式(4)，即得對應(yīng)截面的弦長值。

2 水輪機的設(shè)計參數(shù)

中國潮流能資源分布調(diào)查顯示:流速在1.28～2.04 m/s的潮流蘊含的能量高1.13 GW［8］。盡管某些水道的流速可以達到3 m/s，但總能量較低，開發(fā)價值不大。為使設(shè)計的水輪機能有較大的工作范圍，本文選取V=1.6 m/s為水輪機設(shè)計工作環(huán)境流速。

定義水輪機能量利用率Cp和推力系數(shù)Ct為

水輪機的直徑Z可由下式給出:

式中:Pexp=10 kW為水輪機設(shè)計功率，Cp=0.4為預(yù)測的水輪機能量利用率，η=0.95為裝置的機電效率，ρ=1.025 kg/m3為海水密度。

葉片翼型采用NACA63-8XX系列，該系列翼型具有較高的升力系數(shù)。為保證葉片的結(jié)構(gòu)強度，翼型的相對厚度沿展長方向遞減。輪轂直徑取為1/10葉輪直徑。取葉尖速比λ=6時的轉(zhuǎn)速為設(shè)計轉(zhuǎn)速。

將葉片沿展長劃分為17個站面，根據(jù)約束條件求解目標(biāo)函數(shù)可得各站面的參數(shù)如表1示。

表1 葉片參數(shù)Table 1 Parameters of blade

對于單位翼型，依據(jù)求得的弦長和安裝角進行坐標(biāo)變換得到各剖面的優(yōu)化翼型。為減小輪轂形狀對流場的影響，輪轂兩端采用圓弧面過渡，建立的三維模型如圖2示。

圖2 水輪機三維模型Fig.2 Turbine＇s 3D model

3 數(shù)值計算

3.1 求解區(qū)域網(wǎng)格劃分與設(shè)置

本節(jié)對設(shè)計的葉輪進行CFD模擬以預(yù)測其工作性能，驗證設(shè)計的合理性。水輪機的性能驗證采用商業(yè)軟件 CFX，利用 ICEM-CFD進行網(wǎng)格劃分［9-10］。

本次數(shù)值模擬在長方體形狀的流體域內(nèi)進行，內(nèi)部建立圓柱旋轉(zhuǎn)域以模擬水輪機轉(zhuǎn)動。其中葉輪中心距入口為4倍葉片展長，距出口12倍葉片展長，距四周邊界為3倍葉片展長。對于流體域采用結(jié)構(gòu)網(wǎng)格劃分，由于葉輪表面的形狀很不規(guī)則，旋轉(zhuǎn)域內(nèi)采用非結(jié)構(gòu)網(wǎng)格并在葉輪表面添加邊界層。生成的網(wǎng)格數(shù)據(jù)如表2示。

表2 網(wǎng)格數(shù)據(jù)Table 2 Mesh data

湍流模型采用 Shear Stress Transport模型。SST模型結(jié)合了κ-ε模型與κ-ω模型的優(yōu)點，在靠近壁面部分采用κ-ω模型，在遠處的自由剪切流動采用κ-ε模型。由于考慮了湍動剪切應(yīng)力，SST湍流模型能很好地預(yù)測并計算逆壓梯度下流動的分離。

對于計算域的邊界，入口選擇為速度入口;出口為壓力出口，相對壓力為0，參考壓力為1個大氣壓，且不考慮重力;其余4個流體域邊界設(shè)為墻壁;葉輪表面為無滑移壁面。

3.2 數(shù)值方法的驗證

CFD數(shù)值模擬方法的可靠性可以采用與實驗數(shù)據(jù)對比的方法進行驗證。本文根據(jù)南安普頓大學(xué)公布的三葉片水輪機模型的實驗數(shù)據(jù)進行模擬計算，該模型直徑為0.8 m，采用NACA 63-8XX系列翼型［11］，選取流速V=1.73 m/s時的計算工況，得到的數(shù)值模擬結(jié)果如圖3所示。

圖3 能量利用率和推力系數(shù)對比Fig.3 Comparison of power and thrust coefficients

可以看出:模擬結(jié)果與實驗值較為一致，尤其是在設(shè)計的最佳葉速比附近。最大誤差出現(xiàn)在較高或較低葉速比時，但最大誤差小于10%。由此可以認(rèn)為該方法能很好地模擬水輪機的水動力性能。

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

4.1 水輪機的功率與推力特性

設(shè)定流速為1.6 m/s時，在不同葉速比時水輪機的功率和推力特性計算結(jié)果如圖4所示。

圖4 水輪機功率和推力特性Fig.4 Power and thrust properties of blade

由圖4可知，在初始階段，水輪機的效率隨轉(zhuǎn)速增大而逐漸提高，并在葉速比λ=6附近時達到最大值。說明水輪機滿足了在設(shè)計葉速比下的功率要求。當(dāng)轉(zhuǎn)速高于設(shè)計轉(zhuǎn)速時，水輪機效率開始下降。但與低轉(zhuǎn)速時相比，水輪機在高轉(zhuǎn)速時具有較高的效率。水輪機的推力系數(shù)隨葉速比增大急劇增加。在低轉(zhuǎn)速比時，推力以線性趨勢上升;當(dāng)轉(zhuǎn)速比較高時，速率增長有所放緩。

4.2 葉展長方向單位長度輸出功率分布

由于葉片在展長方向不同截面處的弦長攻角差異，葉片各處受力和輸出的功率不均。圖5為水輪機葉片單位長度輸出功率隨葉片展長的變化。r＜0.4 m的部分為輪轂及輪轂和葉片間的過渡部分，不考慮其輸出功率。

圖5 展長方向單位長度輸出功率Fig.5 Power output along spanwise direction

單位長度輸出功率沿葉片展長大致呈拋物線分布。從0.4 m＜r＜0.9 m的部分單位長度輸出功率呈近乎線性的增長，并在r=1.3 m的附近達到最大值。在0.9 m＜r＜1.7 m的范圍內(nèi)輸出功率始終保持較高的值，但在1.7 m＜r＜2 m的范圍內(nèi)出現(xiàn)急劇下降，在葉尖處降為負值，這一現(xiàn)象和葉尖處流場的三維效應(yīng)有關(guān)。

4.3 三維流動特性分析

在葉片設(shè)計中，葉素理論截取葉片沿展長方向的剖面，將三維葉片簡化為二維剖面進行分析;葉素－動量理論中假定軸向誘導(dǎo)因子沿展長方向不變，即流體無徑向的流動。但在實際三維流場中，壓力的分布不均會導(dǎo)致流體的徑向流動。

為研究三維效應(yīng)對水輪機特性的影響，提取三維流場中半徑為0.4、1.3、1.9、1.99 m的葉片截面壓力值，同時對這些截面處的翼型進行對應(yīng)攻角和流速下的二維模擬計算。

定義壓力系數(shù)Cpre:

各截面的壓力系數(shù)分布如圖6示。可以看出在r=0.4 m處三維流場的壓力系數(shù)分布與二維數(shù)值模擬的結(jié)果差別較大。三維流場中該處壓力面與吸力面的壓力絕對值很低，上下表面的壓差很小，并將導(dǎo)致產(chǎn)生的力矩較低。這是因為靠近輪轂處流場流動復(fù)雜，三維效應(yīng)強烈，葉片上下表面的壓差會誘導(dǎo)出側(cè)向的繞流，從而平衡上下表面的壓力值。

在r=1.3 m截面處的壓力系數(shù)值和二維計算結(jié)果吻合很好。這說明葉片中部的繞流有很強的二維性。但導(dǎo)邊處的壓力值波動較大，這是因為由于該處存在著很大的壓力梯度，而且靠近前駐點的位置流速較低，使得流動易受到臨近截面的干擾而產(chǎn)生沿徑向的流動。

r=1.9 m處三維繞流的壓力分布和二維模擬結(jié)果開始出現(xiàn)差值。具體表現(xiàn)為壓力面的壓力下降，吸力面的壓力上升，導(dǎo)致壓差減小。到r=1.99 m處，壓力分布與r=0.4 m處的已十分相似。繞葉尖端部的繞流(圖7)使得葉片上下面的壓差被抵消，甚至在壓力面出現(xiàn)了負壓。

圖6 二維與三維模擬的壓力系數(shù)比較Fig.6 Pressure coefficient comparison between 2D and 3D simulation

圖7 葉尖繞流流線分布Fig.7 Streamline distribution on blade tip

4.4 葉片表面流場分析

鑒于流場沿葉片展長的變化，對葉片表面流場的流線分析可以給出更為直觀的流動情況。

圖8為葉片表面的流線圖，可以看出:葉展中部的流線較為均勻，且與葉片截面平行。葉根和葉尖部位的流線有偏移，說明該處流場有徑向方向的流動。同時葉片導(dǎo)邊處流場亦有明顯的徑向流動。

圖8 葉片吸力面流線分布Fig.8 Streamline distribution on suction surface

從圖8還可以看出:流線未在葉片表面終止或產(chǎn)生環(huán)形繞流，這說明盡管葉片吸力面有明顯的逆壓梯度，但翼型尾部表面流動并未將至0，附面層的厚度沒有積累起來，未發(fā)生邊界層分離。因此不會導(dǎo)致葉片的失速。從r為0.4、1.3、1.9 m處的截面流線分布(圖9)也可觀察到，葉片周圍繞流比較有規(guī)律。即便是在具有大攻角的葉根(r=0.4 m)處，除了少許沿展長方向的流動引起的在葉片表面的流線終結(jié)，葉片表面的流動表現(xiàn)出很好的層流性質(zhì)。

圖9 葉片截面流線分布Fig.9 Streamline distribution around blade profile

5 結(jié)論

本文基于葉素－動量理論進行了潮流能水輪機的設(shè)計，并用CFX進行數(shù)值模擬證明了水輪機的工作性能達到了設(shè)計要求。由研究結(jié)果可知:

1)葉片設(shè)計的葉素－動量理論忽略展長方向的流動，為準(zhǔn)確預(yù)報水輪機功率特性，在確定葉片弦長、安裝角時，必須考慮三維損失系數(shù)。

2)輸出功率主要由葉片中間部分產(chǎn)生。該段葉片受三維效應(yīng)影響較小。而且由于各截面弦長差別不大，產(chǎn)生功率的大小受截面處線速度影響明顯，靠近外側(cè)部分輸出功率較大。

3)水輪機的三維數(shù)值模擬結(jié)果和對應(yīng)工況下的二維數(shù)值模擬結(jié)果存在一定的誤差，尤其是在葉根、葉尖處誤差很大。因此只有葉片中部的流線與截面相對平行，可采用二維分析;葉尖與葉根處截面的二維分析則失真嚴(yán)重。

4)葉片較小的弦長和較大的線速度成功避免了邊界層分離現(xiàn)象的發(fā)生。使得輸出功率穩(wěn)定，不易產(chǎn)生振動。

本文的分析主要針對水輪機水動力性能和周圍流場的分析，并未考慮載荷對葉片結(jié)構(gòu)的影響和葉片變形對流場的反作用［12］，同時水輪機在真實的工作環(huán)境下還有可能發(fā)生空化現(xiàn)象。這些問題都有待進一步的研究。

［1］FRAENKEL P L.Power from marine currents［J］.Journal of Power Energy，2002，216:1-14.

［2］李志川.垂直軸潮流能水輪機水動力特性數(shù)值模擬和實驗研究［D］.哈爾濱:哈爾濱工程大學(xué)，2011:4-12.

LI Zhichuan.Numerical simulation and experimental study on hydrodynamic performances of vertical axis tidal turbine［D］.Harbin:Harbin Engineering University，2011:4-12.

［3］JO C-H，KIM D Y，RHO Y H，et al.FSI analysis of deformation along offshore pile structure for tidal current power［J］.Renewable Energy，2013，54:248-252.

［4］韓璐.水平軸風(fēng)力機葉片設(shè)計及氣動性能研究［D］.南京:南京航空航天大學(xué)，2008:28-32.

LU Han.Design and experimental research on aerodynamic performances of horizontal axial wind turbine blade［D］.Nanjing:Nanjing University of Aeronautics and Astronautics，2008:28-32.

［5］GAO Qiuxin，JIN Wei，VASSALOS D.The calculations of propeller induced velocity by RANS and momentum theory［J］.Journal of Marine Science and Application，2012，11 (2):164-168.

［6］YAVUZ T，KOC E.Performance analysis of double blade airfoil for hydrokinetic turbine applications［J］.Energy Conversion and Management，2012，63:95-100.

［7］BURTON T，SHARPE D，JENKINS N，et al.風(fēng)能原理［M］.武鑫，譯.北京:科學(xué)出版社，2007:72-75.

［8］WANG Chuankun，LU Wei.Analysis method of ocean energy resource and storage estimation［M］.Beijing:Ocean Publisher，2009:6-9.

［9］HE Miao，WANG Chao，CHANG Xin，et al.Analysis of a propeller wake flow field using viscous fluid mechanics［J］.Journal of Marine Science and Application，2012，11(3): 295-300.

［10］WANG Qiang，ZHOU Hu，WAN Decheng.Numerical simulation of wind turbine blade-tower interaction［J］.Journal of Marine Science and Application，2012，11(3):321-327.

［11］BAHAJ A S，MOLLAND A F，CHAPLIN J R，et al.Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank［J］.Renewable Energy，2007，32(3):407-426.

［12］BAZILEVS Y，HSU M C，SCOTT M A.Isogeometric fluidstructure interaction analysis with emphasis on non-matching discretizations，and with application to wind turbines［J］.Computer Methods in Application and Mechanics Engineering，2012(1):28-41.