潮流能水平軸葉輪縱搖運動水動力分析

張亮，王樹齊，馬勇，張學(xué)偉，國威

（哈爾濱工程大學(xué)深海工程技術(shù)研究中心，黑龍江哈爾濱150001）

潮流能水平軸葉輪縱搖運動水動力分析

張亮，王樹齊，馬勇，張學(xué)偉，國威

（哈爾濱工程大學(xué)深海工程技術(shù)研究中心，黑龍江哈爾濱150001）

漂浮式潮流能裝置載體的波浪運動響應(yīng)使水平軸葉輪的水動力特性發(fā)生變化。采用滑移網(wǎng)格及動網(wǎng)格技術(shù)對無界均勻來流中旋轉(zhuǎn)葉輪強迫縱搖時的三維水動力特性進行分析，研究不同縱搖頻率、縱搖幅值、速比等參數(shù)對葉輪水動力的影響規(guī)律，以及根據(jù)葉輪縱搖運動時的水動力時歷曲線，擬合得到葉輪的縱搖阻尼系數(shù)。研究表明：葉輪軸向載荷和能量利用率瞬時值產(chǎn)生波動，波動幅值隨縱搖頻率、縱搖幅值及速比的增大而增大；葉片表面壓力分布時刻發(fā)生變化，縱搖頻率、縱搖幅值及速比越大，變化越明顯；縱搖阻尼系數(shù)與縱蕩頻率、振幅無關(guān)，而與葉輪旋轉(zhuǎn)速度相關(guān)，旋轉(zhuǎn)速度越大，阻尼系數(shù)越大。

潮流能；水平軸葉輪；強迫縱搖；水動力載荷；阻尼系數(shù)；計算流體力學(xué)

21世紀以來，面對常規(guī)能源日益枯竭、環(huán)境惡化的局面，越來越多的國家開始關(guān)注清潔的可再生能源。其中潮流能［1］沒有明顯的豐枯年以及季節(jié)的大幅度能量變化，具有可持續(xù)、能量密度大、可預(yù)測、同時不占有陸地面積等優(yōu)勢越來越受到人們的青睞。潮流能水輪機［2?4］是潮流能電站的關(guān)鍵部件，它的水動力性能直接關(guān)系到電站的效益。目前，在水平軸水輪機［5?11］水動力性能研究中，基本都是水輪機葉輪只繞固定軸做旋轉(zhuǎn)運動。而在實際漂浮式［12］潮流電站中，漂浮式載體由于波浪的影響會產(chǎn)生運動，從而導(dǎo)致水輪機跟隨漂浮式載體一起運動。例如，當(dāng)漂浮式載體發(fā)生縱搖時，水平軸葉輪也會跟著載體發(fā)生縱搖運動，改變了葉輪旋轉(zhuǎn)平面處的速度分布，進而影響葉輪的水動力載荷。因此研究葉輪在縱搖運動條件下的水動力載荷顯得尤為重要。

本文基于CFX軟件，研究縱搖頻率、幅值、速比等參數(shù)對葉輪水動力載荷的影響規(guī)律；以及對葉輪縱搖運動時的水動力載荷時歷曲線進行擬合來獲得阻尼系數(shù)。為研究漂浮式載體的波浪運動響應(yīng)、水輪機的結(jié)構(gòu)設(shè)計及電能控制等提供相關(guān)數(shù)據(jù)。

1 數(shù)值模擬

1.1 理論基礎(chǔ)

水平軸葉輪旋轉(zhuǎn)中心處的縱搖速度為u，來流速度為U，水輪機直徑為D，葉片數(shù)為N，旋轉(zhuǎn)角速度為ωT，為了方便分析，定義無量綱參數(shù)如下：

式中：λ為葉尖速比，R為葉輪半徑，u-為無因次運動速度，a-為無因次運動加速度，Cz為葉輪軸向載荷系數(shù)，F(xiàn)z為葉輪軸向載荷，Cp為能量利用率，T為葉輪軸向轉(zhuǎn)矩，ρ為介質(zhì)密度。Cz對u-、u-2、a-和a-2取變分可得

式中：nzz為阻尼系數(shù)，mzz為附加值量系數(shù)。考慮到葉片數(shù)N的影響，當(dāng)水平軸水輪機以頻率ωT旋轉(zhuǎn)時，水輪機整體水動力以N倍頻變化，因此式（1）中的各項水動力系數(shù)都應(yīng)是以為基頻的周期函數(shù)。

假設(shè)水動力系數(shù)可寫成級數(shù)形式：

由于動態(tài)失速以及尾流場的存在，水平軸水輪機的流場非常復(fù)雜，很難采用直接的方法來求解相關(guān)水動力系數(shù)。本文基于最小二乘法原理，對葉輪縱搖運動時的水動力載荷時歷曲線進行擬合來獲得各項縱搖系數(shù)。

1.2 CFX數(shù)值模擬

1.2.1 建立計算域

計算模型的葉片由哈爾濱工程大學(xué)海洋可再生能源研究所研發(fā)，除葉片根部外，展向各位置處的截面翼型為S809。葉輪直徑0.7 m，葉片數(shù)為2，輪轂直徑為0.1D，如圖1所示。在CFX數(shù)值模擬時，需要同時模擬水輪機的旋轉(zhuǎn)運動和縱搖運動（a＝Asin（ωt），A為縱搖幅值，ω為縱搖頻率），因此將整個計算域分為2個區(qū)域：靜止域和旋轉(zhuǎn)域，旋轉(zhuǎn)域為圓柱體，旋轉(zhuǎn)域又分為縱搖旋轉(zhuǎn)域和葉輪旋轉(zhuǎn)域，葉輪旋轉(zhuǎn)域在縱搖旋轉(zhuǎn)域的內(nèi)部，葉輪旋轉(zhuǎn)軸距縱搖中心0.4 m，計算域如圖1所示。

圖1 模型及網(wǎng)格示意圖Fig.1 Schematic of blade and mesh models

1.2.2 劃分網(wǎng)格

網(wǎng)格化分是數(shù)值模擬技術(shù)的一個非常重要的組成部分，網(wǎng)格質(zhì)量的好壞直接影響到數(shù)值計算的精度和效率。整個計算域的網(wǎng)格劃分采用結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格結(jié)合的方式，葉輪旋轉(zhuǎn)域由于葉片外形復(fù)雜，采用非結(jié)構(gòu)化網(wǎng)格形式，縱搖旋轉(zhuǎn)域和靜止域采用結(jié)構(gòu)化網(wǎng)格形式。葉輪表面網(wǎng)格的第1層高度為0.000 5 m，網(wǎng)格總數(shù)為201×104個。如圖1所示。

1.2.3 邊界條件設(shè)置

設(shè)定大氣壓為參考氣壓；將進口邊界設(shè)為速度進口，給定均勻來流速度U和湍流參數(shù)；出口邊界采用壓力出口邊界條件；流體計算域的側(cè)面設(shè)定為自由滑動墻面；葉輪表面為無滑移壁面；旋轉(zhuǎn)域給定縱搖角速度（滑移網(wǎng)格形式），葉輪旋轉(zhuǎn)域又以子域的形式給定旋轉(zhuǎn)運動（動網(wǎng)格形式），即葉輪旋轉(zhuǎn)域在旋轉(zhuǎn)運動的同時又和縱搖旋轉(zhuǎn)域一起產(chǎn)生縱搖運動，網(wǎng)格不發(fā)生變形；旋轉(zhuǎn)域和靜止域及葉輪旋轉(zhuǎn)域和縱搖旋轉(zhuǎn)域之間通過交界面連接。計算中湍流模型采用的是SST模型，求解器為瞬態(tài)求解器，時間步長為葉輪旋轉(zhuǎn)3°所用的時間。

2 葉輪載荷分析

在同樣的網(wǎng)格模型下，水平軸葉輪在無界均勻來流中僅做旋轉(zhuǎn)運動時的能量利用率及軸向載荷系數(shù)曲線如圖2所示。圖中實驗值是在哈爾濱工程大學(xué)船舶拖曳水池中進行測試的，由圖2可知，計算值與實驗值吻合很好，驗證了此網(wǎng)格模型模擬水平軸葉輪的準確性，下面對縱搖進行分析。

圖2 計算值與實驗值對比曲線Fig.2 Comparison between calculated and experimen?tal values

2.1 縱搖頻率對葉輪載荷的影響

選取縱搖頻率在常見波浪頻率范圍內(nèi)，即取縱搖頻率為0.4、0.8、1.2、1.6、2.0 rad／s。固定葉輪轉(zhuǎn)速230 r／min，縱搖幅值取15°，入口速度1.5 m／s。圖3給出了不同縱搖頻率下葉輪軸向載荷系數(shù)及能量利用率隨時間的變化曲線，ω＝0表示無縱搖的情況。由圖3可知，當(dāng)不發(fā)生縱搖時，葉輪的軸向載荷系數(shù)及能量利用率不隨時間推移而改變。而當(dāng)葉輪發(fā)生縱搖運動時，葉輪軸向載荷系數(shù)及能量利用率隨著葉輪的縱搖而產(chǎn)生明顯波動，縱搖頻率越大，波動幅值越大。當(dāng)縱搖頻率為2.0 rad／s時，葉輪軸向載荷系數(shù)和能量利用率波動幅度分別達到葉輪不發(fā)生縱搖時的45.7%和103.2%左右。原因是縱搖改變了葉輪旋轉(zhuǎn)平面處的速度分布，相當(dāng)于來流速度的大小及方向時刻發(fā)生變化，從而使葉輪載荷產(chǎn)生波動。從圖中還可以看出，在葉輪旋轉(zhuǎn)過程中，葉輪的軸向載荷系數(shù)及能量利用率會發(fā)生微幅波動，這是由于葉輪旋轉(zhuǎn)平面不是正對來流，最大存在15°的夾角。因此，在實際海洋環(huán)境中，當(dāng)遭遇高頻波浪時，葉輪載荷將發(fā)生高頻率、高幅值波動，這對水輪機的電能輸出的穩(wěn)定性及結(jié)構(gòu)安全可靠性來說是非常不利的。

圖3 不同縱搖頻率瞬時值對比Fig.3 Instantaneous value of different oscillation fre?quency

2.2 縱搖幅值對葉輪載荷的影響

取水平軸葉輪轉(zhuǎn)速230 r／min，水輪機縱搖頻率為1.2 rad／s，取3種不同縱搖幅值A(chǔ)（10°、15°和20°）進行計算，軸向載荷系數(shù)及能量利用率隨時間的變化曲線如圖4所示，從圖4可以明顯看出，縱搖幅值越大，軸向載荷系數(shù)及能量利用率瞬時值的波動幅值越大，這是由相對來流方向的改變引起的。在實際潮流電站中，水輪機運動與浮體平臺的運動是相互耦合的，因此在潮流電站的設(shè)計中應(yīng)盡可能減小浮體平臺的運動幅度，以提高水輪機的結(jié)構(gòu)安全性。

圖4 不同縱搖振幅瞬時值對比Fig.4 Instantaneous values of different oscillation am?plitude

2.3 速比對葉輪載荷的影響

考慮葉輪在縱搖條件下不同速比時葉輪載荷的變化規(guī)律。取縱搖頻率為1.2 rad／s，縱搖幅值為15°，來流速度1.5 m／s，取7種不同轉(zhuǎn)速進行計算，軸向載荷系數(shù)及能量利用率瞬時值隨時間的變化曲線如圖5所示，從圖5可以看出，在同一縱搖頻率及幅值條件下，隨著葉輪轉(zhuǎn)速的增加，葉輪軸向載荷系數(shù)及能量利用率波動頻率保持不變，但波動幅值增大，這樣不利于電能的穩(wěn)定輸出及水輪機的結(jié)構(gòu)安全。

圖5 不同速比瞬時值對比Fig.5 Instantaneous values of different tip speed ratios

2.4 翼型表面壓力分析

為了研究縱搖運動時載荷發(fā)生波動的機理，對葉片表面壓力分布進行分析，對不同縱蕩幅值，分別提取2個縱搖位置（1＃位置是縱搖平衡位置，縱搖速度與來流速度相反，葉輪處的實際流速是來流速度與縱搖速度之和；2＃位置是縱搖平衡位置，縱搖速度與來流速度相同，葉輪處的實際流速是來流速度與縱搖速度之差）時的葉片表面壓力，圖6給出了葉片0.4R處時的翼型表面壓力分布。圖中l(wèi)是截面翼型弦長坐標歸一化，Cpre表示表面相對壓力系數(shù)，表達式為

式中：p為絕對壓力，p0為參考大氣壓。

圖6中，1＃、2＃分別表示1＃位置、2＃位置，A表示縱搖振幅，A＝0表示縱搖振幅為0，即無縱搖。從圖6中可以看出，當(dāng)葉輪在1＃位置時，隨著縱蕩幅值的增加，葉片壓力面的壓力逐漸增大，吸力面的壓力逐漸減小，從而使葉片表面壓力差隨著幅值的增大而增大，而當(dāng)葉輪在2＃位置時，隨著縱蕩幅值的增加，葉片壓力面的壓力逐漸減小，吸力面的壓力逐漸增大，從而使葉片表面壓力差隨著幅值的增大而減小。同理，當(dāng)葉輪在1＃位置時，葉片表面壓力差隨著縱搖頻率的增大而增大，而當(dāng)葉輪在2＃位置時，葉片表面壓力差隨著頻率的增大而減小。因此當(dāng)葉輪發(fā)生縱搖運動時，葉輪水動力載荷會發(fā)生波動，且波動幅值的大小隨頻率及幅值的增加而增加。

圖6 表面壓力系數(shù)分布Fig.6 Distribution of the airfoil surface pressure

3 縱搖阻尼系數(shù)擬合

3.1 縱搖頻率對阻尼系數(shù)的影響

對不同縱搖頻率下軸向載荷的時歷曲線進行最小二乘法擬合，得出Cz級數(shù)形式的各項系數(shù)，如表1，表中只給出了Cz展開中前3項的常數(shù)項和一次項。

表1 Cz展開系數(shù)Table1 Czexpansion coefficients

從表中可以看出，Cz展開的各項系數(shù)中，前2項的常數(shù)項遠大于一次項系數(shù)，且隨著縱蕩頻率的增加，常數(shù)項基本不發(fā)生變化，因此，縱搖頻率對阻尼系數(shù)沒有影響。

3.2 縱搖幅值對阻尼系數(shù)的影響

對不同縱搖振幅下軸向載荷的時歷曲線進行最小二乘法擬合，得出Cz級數(shù)形式的各項系數(shù)，如表2，表中只給出了Cz展開中前3項的常數(shù)項和一次項。

表2 Cz展開系數(shù)Table2 Czexpansion coefficients

從表中可以得到與3.1節(jié)相同的結(jié)論，即Cz展開的各項系數(shù)中，前2項的常數(shù)項遠大于一次項系數(shù)，隨著縱搖振幅的增加，常數(shù)項基本不發(fā)生變化，因此，縱搖振幅對阻尼系數(shù)沒有影響。

3.3 速比對阻尼系數(shù)的影響

對不同轉(zhuǎn)速（速比）下軸向載荷的時歷曲線進行最小二乘法擬合，得出Cz級數(shù)形式的各項系數(shù)，如表3，表中只給出了Cz展開中前3項的常數(shù)項和一次項。

表3 Cz展開系數(shù)Table3 Czexpansion coefficients

從表中可以看出，Cz展開的各項系數(shù)中，前兩項的常數(shù)項遠大于一次項系數(shù)，隨著轉(zhuǎn)速（速比）的增加，阻尼系數(shù)常數(shù)項（絕對值）逐漸增加，因此，轉(zhuǎn)速對阻尼系數(shù)影響明顯，即轉(zhuǎn)速越大，阻尼系數(shù)越大。

4 結(jié)論

采用CFX對三維水平軸葉輪在均勻來流中強迫縱搖時的水動力特性進行了計算分析，研究結(jié)果表明：

1）葉輪縱搖時，葉輪軸向載荷及能量利用率產(chǎn)生明顯波動，波動幅值隨著縱搖頻率、縱搖幅值及速比的增加而增加，這對葉輪的結(jié)構(gòu)安全、疲勞壽命及潮流電站電能的穩(wěn)定輸出產(chǎn)生影響，在結(jié)構(gòu)和電控系統(tǒng)設(shè)計中需要特別考慮。

2）葉輪縱搖時，葉輪旋轉(zhuǎn)平面處的速度大小及方向時刻發(fā)生變化，從而使葉片表面壓力分布發(fā)生明顯改變，這是導(dǎo)致葉輪載荷及能量利用率產(chǎn)生波動的原因所在。

3）根據(jù)最小二乘法原理，對葉輪縱搖時的軸向載荷進行數(shù)值擬合，得出縱搖時的阻尼系數(shù)，阻尼系數(shù)的大小與縱搖頻率、幅值無關(guān)，而與葉輪旋轉(zhuǎn)速度相關(guān)，旋轉(zhuǎn)速度越大，阻尼系數(shù)越大。這為研究漂浮式載體的運動響應(yīng)提供了重要的數(shù)據(jù)。

［1］呂忻，郭佩芳.我國潮流能資源開發(fā)評述［J］.海洋湖沼通報，2011（1）：26?30.LYU Xin，GUO Peifang.Review of China tidal current ener?gy resource development［J］.Transactions of Oceanology and Limnology，2011（1）：26?30.

［2］戴軍，單忠德，王西峰，等.潮流水輪機的研究進展［J］.可再生能源，2010，28（4）：130?133.DAI Jun，SHAN Zhongde，WANG Xifeng，et al.Current re?search progress of water turbine［J］.Renewable Energy.2010，28（4）：130?133.

［3］陳展，馬勇，張亮，等.矩形潮流能水輪機性能研究［J］.華中科技大學(xué)學(xué)報：自然科學(xué)版，2013，6（41）：128?132.CHEN Zhan，MA Yong，ZHANG Liang，et al.Study on properties of rectangular tidal turbine［J］.Journal of Hua?zhong University of Science and Technology：Natural Science Edition，2013，6（41）：128?132.

［4］LEE J H，PARK S，KIM D H，et al.Computational meth?ods for performance analysis of horizontal axis tidal stream turbines［J］.Applied Energy，2012，98：512?523.

［5］BATTEN W，BAHAJ A S，MOLLAND A F，et al.The pre?diction of the hydrodynamic performance of marine current turbines［J］.Renewable Energy，2008，33（5）：1085?1096.［6］LI Y，CAHSAL S M.A discrete vortex method for simulating a stand?alone tidal?current turbine：modeling and validation［J］.Journal of Offshore Mechanics and Arctic Engineering，2010，132（3）：1?9.

［7］CHOI H J，ZULLAH M A，ROH H W，et al.CFD validation of performance improvement of a 500 kW Francis turbine［J］.Renewable Energy，2013，54（SI）：111?123.

［8］BAHAJ A S，MOLLAND A F，CHAPLIN J R，et al.Power and thrust measurements of marine current turbines under va?rious hydrodynamic flow conditions in a cavitation tunnel and a towing tank［J］.Renewable Energy，2007，32（3）：407?426.

［9］MYERS L E，BAHAJ A S.Experimental analysis of the flow field around horizontal axis tidal turbines by use of scale mesh disk rotor simulators［J］.Ocean Engineering，2010，37：218?227.

［10］沈云，李龍，朱多彪.水平軸潮流能水輪機轉(zhuǎn)輪尾流特性數(shù)值分析［J］.水電能源科學(xué)，2013，31（10）：149?151.SHEN Yun，LI Long，ZHU Duobiao.Numerical analysis of wake performance of horizontal?axis tidal current turbine［J］.Water Resources and Power，2013，31（10）：149?151.

［11］王曉航.潮流能水平軸葉輪設(shè)計與性能預(yù)測方法研究［D］.哈爾濱：哈爾濱工程大學(xué)，2011：3?8.WANG Xiaohang.Study on the design and performance pre?diction of tidal current horizontal axis turbine［D］.Harbin：Harbin Engineering University，2011：3?8.

［12］馬勇，張亮，盛其虎，等.漂浮式潮流能發(fā)電裝置水動力特性試驗研究［J］.華中科技大學(xué)學(xué)報：自然科學(xué)版，2012，10（40）：123?127.MA Yong，ZHANG Liang，SHENG Qihu，et al.The test study on hydrodynamic characteristics of floating tidal power generation devices［J］.Journal of Huazhong University of Science and Technology：Natural Science Edition，2012，10（40）：123?127.

The pitch hydrodynamic analysis of tidal current energy horizontal axis impeller

ZHANG Liang，WANG Shuqi，MA Yong，ZHANG Xuewei，GUO Wei
（Deepwater Engineering and Technology Research Center，Harbin Engineering University，Harbin 150001，China）

Under the condition of actual sea state，hydrodynamic characteristic of floating horizontal axis impeller is related to wave characteristics and floating carrier wave motion response.This paper uses slipping mesh and dynam?ic mesh to analyze the hydrodynamic characteristic in uniform unbounded stream when the rotating impeller is forced to pitch，and studies the influence of different pitch frequencies，pitch amplitudes and tip speed ratios on turbine hydrodynamics.Pitch damping coefficient can be got by hydrodynamic time?history curve of pitching impeller.The results showed that the impeller axial load and energy utilization ratio will make a fluctuation in pitching and the wave amplitude of momentary value increases with the increasing of pitch frequency，pitch amplitude，and tip speed ratio.The blade surface pressure distribution changes when the impeller is pitching，and the larger the pitch fre?quency，pitch amplitude and tip speed ratio，the more obvious the changes.The frequency and amplitude of the pitch have little impact on damping coefficient，but this coefficient is positively proportional to the rotational speed of the impeller.

tidal current energy；horizontal axis impeller；forced pitching；hydrodynamic load；damping coeffi?cient；computational fluid dynamics（CFD）

10.3969／j.issn.1006?7043.201403015

http：／／www.cnki.net／kcms／detail／23.1390.U.20150109.1505.006.html

TK730；O352

A

1006?7043（2015）03?0307?05

2014?03?06.網(wǎng)絡(luò)出版時間：2015?01?09.

國家自然科學(xué)基金資助項目（51309069，51309068，51106034）；海洋能專項基金資助項目（GHME2010GC02）；哈爾濱市科技創(chuàng)新人才研究專項資金資助項目（RC2014QN001008，2012RFQXG086）；中國博士后科學(xué)基金面上資助項目（2014M561334）；黑龍江省博士后基金資助項目（LBH?Z14060）.

張亮（1959?），男，教授，博士生導(dǎo)師；馬勇（1980?），男，講師，博士.

馬勇，E?mail：mayong02＠hrbeu.edu.cn.