水平軸潮流能水輪機設計與數值模擬

盛其虎，趙東亞，張亮

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

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

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

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

1 葉片設計的葉素動量理論

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

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

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

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

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

定義葉尖速比λ:

式中:R為水輪機半徑。

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

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

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

沿葉片展長方向取不同的截面，以產生的功率最大為設計優化目標。在各界面處，由單位葉素做功可得目標函數:

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

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

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

2 水輪機的設計參數

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

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

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

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

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

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

表1 葉片參數Table 1 Parameters of blade

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

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

3 數值計算

3.1 求解區域網格劃分與設置

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

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

表2 網格數據Table 2 Mesh data

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

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

3.2 數值方法的驗證

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

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

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

4 計算結果與分析

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

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

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

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

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的范圍內輸出功率始終保持較高的值，但在1.7 m＜r＜2 m的范圍內出現急劇下降，在葉尖處降為負值，這一現象和葉尖處流場的三維效應有關。

4.3 三維流動特性分析

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

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

定義壓力系數Cpre:

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

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

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

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

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

4.4 葉片表面流場分析

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

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

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

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

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

5 結論

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

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

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

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

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

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

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

［2］李志川.垂直軸潮流能水輪機水動力特性數值模擬和實驗研究［D］.哈爾濱:哈爾濱工程大學，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］韓璐.水平軸風力機葉片設計及氣動性能研究［D］.南京:南京航空航天大學，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.風能原理［M］.武鑫，譯.北京:科學出版社，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.