張錦嵐，王 凱，張萬超，周亞輝

(1. 武漢第二船舶設(shè)計(jì)研究所， 湖北 武漢 430064；2. 江蘇科技大學(xué) 船舶與海洋工程學(xué)院， 江蘇 鎮(zhèn)江 212003；3. 江蘇科技大學(xué) 能源與動(dòng)力學(xué)院， 江蘇 鎮(zhèn)江 212003)

海洋可再生能源儲(chǔ)量豐富，總儲(chǔ)量可達(dá)7.66億千瓦，其開發(fā)利用逐漸受到世界各國(guó)的重視[1-3]。其中潮流能因其分布相對(duì)集中、能流密度高、有規(guī)律可循和相對(duì)穩(wěn)定的屬性特點(diǎn)[4-6]，受到更多的重視。目前潮流能開發(fā)利用，已經(jīng)進(jìn)入復(fù)雜海況下的定常流研究階段。通過限定水輪機(jī)運(yùn)動(dòng)特性，探究復(fù)雜特殊工況下的水輪機(jī)水動(dòng)力特性及能量轉(zhuǎn)換是目前急需解決的任務(wù)。立軸葉輪強(qiáng)迫晃蕩耦合運(yùn)動(dòng)是解決上述問題的關(guān)鍵，但相關(guān)研究成果還比較有限[7-9]。

目前為止，利用計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)方法模擬葉輪的運(yùn)行工況已經(jīng)比較成熟。2008年，Nabavi[10]應(yīng)用FLUENT軟件模擬了固定偏角水輪機(jī)的性能，并采用Spalart-Ailmaras湍流模型和滑移網(wǎng)格技術(shù)模擬水輪機(jī)二維及三維運(yùn)動(dòng)，此外還分析了時(shí)間步、邊界距離、葉片處網(wǎng)格厚度對(duì)模擬結(jié)果的影響。2014年，Danao等[11]應(yīng)用FLUENT軟件分析了二維和三維垂直軸風(fēng)機(jī)在穩(wěn)態(tài)風(fēng)作用下的性能，湍流度、葉片表面節(jié)點(diǎn)數(shù)、邊界距離及最小時(shí)間步等因素被討論，算例采用滑移網(wǎng)格技術(shù)，提出時(shí)間步應(yīng)取風(fēng)機(jī)轉(zhuǎn)過0.5°所需時(shí)間。2015年，Alaimo等[12]應(yīng)用FLUENT軟件模擬了二維及三維垂直軸水輪機(jī)在不同速比下的運(yùn)動(dòng)過程，著重分析了模擬時(shí)間步、模型網(wǎng)格數(shù)、葉片處網(wǎng)格尺寸等因素對(duì)模擬結(jié)果的影響。模型采用滑移網(wǎng)格技術(shù)且各算例均采用k-ε湍流模型。這些學(xué)者研究了葉輪三維和二維模擬下的葉輪水動(dòng)力性能影響。對(duì)于立軸葉輪在浪流作用下發(fā)生振蕩運(yùn)動(dòng)的三維和二維水動(dòng)力性能差別還少有學(xué)者進(jìn)行研究。

近幾年，一些學(xué)者對(duì)浪流作用下的水平軸水輪機(jī)的水動(dòng)力性能做了一些研究[13-16]，這些學(xué)者認(rèn)為浪流作用下水平軸水輪機(jī)的水動(dòng)力會(huì)發(fā)生波動(dòng)，力和力矩的峰值會(huì)提高，但是平均值基本不變。哈爾濱工程大學(xué)張亮團(tuán)隊(duì)從潮流能水平軸水輪機(jī)實(shí)際運(yùn)行中也發(fā)現(xiàn)浪流作用對(duì)水輪機(jī)的水動(dòng)力性能有明顯影響。

1 數(shù)值模擬

1.1 理論基礎(chǔ)

表1為立軸水輪機(jī)振蕩運(yùn)動(dòng)下的各參數(shù)定義。表1除去常規(guī)的葉輪參數(shù)定義外，特別定義了葉輪振蕩運(yùn)動(dòng)下的振蕩幅值、振蕩頻率和振蕩速度。

表1 立軸水輪機(jī)葉輪參數(shù)

通過表1，定義無量綱參數(shù)如下：

推力系數(shù)

CFX=FX/(0.5ρ·V2·D·H)

(1)

側(cè)向力系數(shù)

CFY=FY/(0.5ρ·V2·D·H)

(2)

能量效率系數(shù)

Cp=Q·ω/(0.5ρ·V3·D·H)

(3)

葉尖速比

λ=R·ω/V

(4)

無因次振蕩速度

(5)

無因次振蕩加速度

(6)

為了簡(jiǎn)化描述，使用變量Φ來表示坐標(biāo)軸方向，其中Φ為X時(shí)表示X軸方向，Φ為Y時(shí)表示Y軸方向。葉輪的受力可以分解為三類，即均勻水動(dòng)力、阻尼力和附加質(zhì)量力，則葉輪受力表示為

(7)

對(duì)阻尼力和附加質(zhì)量力進(jìn)行分解，表示為

(8)

(9)

其中，ψ表示葉輪位置角和各項(xiàng)展開項(xiàng)的錯(cuò)位角。將式(8)和式(9)代入式(7)，則葉輪振蕩運(yùn)動(dòng)下的受力可表示為

(10)

利用ANSYS-CFX軟件對(duì)立軸葉輪振蕩運(yùn)動(dòng)進(jìn)行模擬仿真，并將得到的葉輪振蕩運(yùn)動(dòng)下的推力和側(cè)向力時(shí)歷曲線利用最小二乘法擬合，就可以得到式(10)中的各項(xiàng)系數(shù)，總結(jié)對(duì)比二維模擬和三維模擬各參數(shù)對(duì)立軸葉輪振蕩運(yùn)動(dòng)的影響規(guī)律。

1.2 立軸水輪機(jī)模型

選用雙葉片葉輪模型，其葉片采用NACA0018對(duì)稱翼型，葉輪直徑為0.8 m，弦長(zhǎng)為0.12 m，來流速度為1.0 m/s。

圖1和圖2分別為立軸振蕩葉輪數(shù)值計(jì)算的二維和三維網(wǎng)格模型。本文主要計(jì)算葉輪旋轉(zhuǎn)運(yùn)動(dòng)和振蕩運(yùn)動(dòng)相結(jié)合工況下的二維和三維結(jié)果，其中振蕩運(yùn)動(dòng)包括橫蕩運(yùn)動(dòng)和縱蕩運(yùn)動(dòng)。兩幅圖中分別給出了二維和三維工況下的各運(yùn)動(dòng)形式計(jì)算域，主要包括旋轉(zhuǎn)域、振蕩域和隨振蕩運(yùn)動(dòng)網(wǎng)格發(fā)生微幅變形的變形域。

圖3為葉輪振蕩運(yùn)動(dòng)模擬設(shè)置的邊界條件。邊界條件的設(shè)置和常規(guī)旋轉(zhuǎn)葉輪的邊界條件設(shè)置類似，所不同的是本文中將旋轉(zhuǎn)域以外的流場(chǎng)域分為振蕩域和振蕩變形域，其中振蕩域在計(jì)算中發(fā)生橫蕩或者縱蕩運(yùn)動(dòng)，運(yùn)動(dòng)過程中網(wǎng)格只存在滑移運(yùn)動(dòng)，網(wǎng)格不發(fā)生變形；振蕩變形域在計(jì)算中網(wǎng)格發(fā)生微幅拉伸和壓縮，但是對(duì)網(wǎng)格質(zhì)量幾乎沒有影響。

(a) 葉片(a) Blade (b) 旋轉(zhuǎn)域(b) Rotation domain

(c) 橫蕩域(c) Sway region (d) 橫蕩變形域(d) Sway deformation region

(e) 縱蕩域(e) Surge region (f) 縱蕩變形域(f) Surge deformation region 圖1 二維網(wǎng)格模型Fig.1 2D grid model

(a) 葉片(a) Blade (b) 旋轉(zhuǎn)域(b) Rotation domain

(c) 橫蕩域(c) Sway region (d) 橫蕩變形域(d) Sway deformation region

(e) 縱蕩域(e) Surge region (f) 縱蕩變形域(f) Surge deformation region圖2 三維網(wǎng)格模型Fig.2 3D grid model

(a) 二維橫蕩(a) Two-dimensional sway motion

(b) 二維縱蕩(b) Two-dimensional surge motion

(c) 三維橫蕩(c) Three-dimensional sway motion

(d) 三維縱蕩(d) Three-dimensional surge motion圖3 邊界條件Fig.3 Boundary conditions

2 結(jié)果分析

2.1 CFD方法的可行性證明

在計(jì)算模擬前，首先驗(yàn)證立軸葉輪仿真計(jì)算的有效性，分別從二維模擬和三維模擬方面對(duì)計(jì)算值和試驗(yàn)值進(jìn)行對(duì)比分析。

二維計(jì)算驗(yàn)證：采用本文的計(jì)算方法對(duì)美國(guó)Strikland教授做的關(guān)于立軸葉輪的二維試驗(yàn)[17]進(jìn)行模擬，并將計(jì)算得到的單個(gè)葉片切向力、法向力與試驗(yàn)結(jié)果對(duì)比。圖4為試驗(yàn)值和本文計(jì)算結(jié)果的對(duì)比，從兩幅圖中可以看出，計(jì)算結(jié)果和試驗(yàn)值相近。教授也對(duì)該葉輪進(jìn)行了數(shù)值模擬計(jì)算，本文結(jié)果和其計(jì)算結(jié)果類似，在90°～270°的上游盤面，計(jì)算值和試驗(yàn)值誤差較小，而270°～450°的下游盤面誤差較大。Strikland教授分析可能是：一方面試驗(yàn)時(shí)葉片安裝有問題，造成下游盤面試驗(yàn)值存在凸起現(xiàn)象；另一方面下游盤面水流受到上游盤面影響，水流復(fù)雜，導(dǎo)致試驗(yàn)值和計(jì)算值存在一定差異。

為了驗(yàn)證本文計(jì)算方法三維模擬的有效性，對(duì)比了哈爾濱工程大學(xué)“HEU”型兩葉片立軸水輪機(jī)的水動(dòng)力性能，結(jié)果如圖5所示。圖5(a)和圖5(b)分別為葉輪效率和推力的試驗(yàn)值和三維計(jì)算值對(duì)比。不論是效率還是推力，三維仿真計(jì)算結(jié)果都略大于試驗(yàn)值，這是因?yàn)榉抡嬗?jì)算中未考慮葉輪運(yùn)行過程中機(jī)械摩擦、壁效應(yīng)的影響，這些影響都會(huì)造成葉輪效率和推力的降低。綜合來看，立軸葉輪三維模擬計(jì)算精度較高，能夠模擬立軸葉輪真實(shí)運(yùn)行情況。

(a) 切向力(a) Tangential force

(b) 法向力(b) Normal force圖4 二維計(jì)算對(duì)比圖Fig.4 Comparisons between calculated and experimental values

(a) 效率對(duì)比(a) Comparison of efficiency

(b) 推力系數(shù)對(duì)比(b) Comparison of thrust coefficient圖5 三維計(jì)算對(duì)比Fig.5 Results comparison between the experimental and 3D numerical simulation

2.2 立軸葉輪振蕩運(yùn)動(dòng)效率

當(dāng)展弦比較小時(shí)，葉輪振蕩運(yùn)動(dòng)下二維計(jì)算和三維計(jì)算結(jié)果相差較大。為了尋找不同展弦比下二維結(jié)果和三維結(jié)果的變化規(guī)律，首先計(jì)算了不同展弦比的立軸葉輪振蕩運(yùn)動(dòng)，并將效率作為對(duì)比對(duì)象和二維結(jié)果進(jìn)行比較。如圖6所示，兩圖為橫蕩運(yùn)動(dòng)和縱蕩運(yùn)動(dòng)下不同展弦比三維葉輪效率和二維效率對(duì)比的結(jié)果。隨著展弦比增大，三維計(jì)算結(jié)果逐漸增大，當(dāng)展弦比達(dá)到10左右時(shí)，三維計(jì)算結(jié)果和二維計(jì)算結(jié)果相近。葉輪三維計(jì)算結(jié)果和二維計(jì)算結(jié)果存在差異的最主要原因是三維計(jì)算中存在三維效應(yīng)，葉片展長(zhǎng)兩側(cè)存在機(jī)翼繞流，葉梢處高壓區(qū)流體會(huì)流向低壓區(qū)，造成葉梢處尾渦，帶走一部分能量，因此葉輪三維效率計(jì)算結(jié)果小于二維計(jì)算結(jié)果。

(a) 橫蕩運(yùn)動(dòng)(a) Sway motion

(b) 縱蕩運(yùn)動(dòng)(b) Surge motion圖6 不同展弦比的水輪機(jī)效率曲線Fig.6 Comparison of power output curves of vertical axis turbine with different span-chord ratios

2.3 立軸葉輪振蕩運(yùn)動(dòng)推力和側(cè)向力

圖7所示為葉輪橫蕩和縱蕩運(yùn)動(dòng)下的不同展弦比的推力系數(shù)和側(cè)向力系數(shù)的對(duì)比。其中圖7(a)和圖7(c)為推力系數(shù)對(duì)比，圖7(b)和圖7(d)為側(cè)向力系數(shù)對(duì)比。四幅圖中分別對(duì)比了展長(zhǎng)為0.2 m、0.6 m和1.0 m的三維計(jì)算結(jié)果和二維計(jì)算結(jié)果。從四幅圖中可以看出，隨著展長(zhǎng)增大，相對(duì)應(yīng)的推力和側(cè)向力系數(shù)逐漸接近二維計(jì)算結(jié)果。推力和側(cè)向力的波峰和波谷都隨著展長(zhǎng)增大而增大。這是因?yàn)檎归L(zhǎng)越長(zhǎng)，葉輪的三維效應(yīng)影響越小，葉梢處尾渦帶走的能量越少，計(jì)算結(jié)果越接近二維計(jì)算值。

(a) 橫蕩推力(a) Thrust of sway

(b) 橫蕩側(cè)向力(b) Lateral force of sway

(c) 縱蕩推力(c) Thrust of surge

(d) 縱蕩側(cè)向力(d) Lateral force of surge圖7 不同展長(zhǎng)的葉輪振蕩狀態(tài)下的推力和側(cè)向力系數(shù)對(duì)比Fig.7 Comparison of thrust and lateral force coefficients under swaying motion and surging motion of different elongation

圖8為三個(gè)不同展長(zhǎng)的葉輪流場(chǎng)圖，從三幅圖中可以看出上游葉片葉梢兩端產(chǎn)生流體分離，帶走了一部分能量。圖8(a)所示展長(zhǎng)0.2 m的葉輪流體分離幾乎布滿整個(gè)葉片，明顯大于圖8(b)和圖8(c)。結(jié)合圖9所示的不同展長(zhǎng)的葉輪渦量場(chǎng)，在同一渦量指標(biāo)下，展長(zhǎng)0.2 m的葉輪兩端渦量強(qiáng)度明顯大于0.4 m和0.6 m展長(zhǎng)的渦量強(qiáng)度，其帶走的能量也越多，因此導(dǎo)致葉輪水動(dòng)力性能降低。

從圖8和圖9可以看出，展弦比越大，葉輪由于三維效應(yīng)在葉片兩端產(chǎn)生的尾渦越小，帶走的能量也越小，因此葉輪效率、推力等水動(dòng)力系數(shù)越高。結(jié)合2.2節(jié)，當(dāng)葉輪的展弦比達(dá)到10左右時(shí)，葉輪尾渦帶走的能量可以忽略不計(jì)，水動(dòng)力性能也與二維葉輪的水動(dòng)力性能相近。

(a) H=0.2 m (b) H=0.4 m (c) H=0.6 m圖8 不同展長(zhǎng)的葉輪流場(chǎng)流線圖Fig.8 Streamlines of the flow field under different elongations

(a) H=0.2 m (b) H=0.4 m (c) H=0.6 m圖9 不同展長(zhǎng)的葉輪渦量場(chǎng)Fig.9 Vorticity fields under different elongations

3 水動(dòng)力系數(shù)

根據(jù)式(10)將三維模擬計(jì)算得到的不同展長(zhǎng)的立軸葉輪推力和側(cè)向力系數(shù)進(jìn)行擬合分析，得到各受力系數(shù)的均勻水動(dòng)力、阻尼力和附加質(zhì)量力系數(shù)。如表2～5所示，對(duì)比了0.2 m、0.4 m、0.6 m、0.8 m四個(gè)展長(zhǎng)和2D模擬橫蕩與縱蕩運(yùn)動(dòng)下的各項(xiàng)系數(shù)。

表2 橫蕩運(yùn)動(dòng)CFX展開系數(shù)

表3 縱蕩運(yùn)動(dòng)CFX展開系數(shù)

表4 橫蕩運(yùn)動(dòng)CFY展開系數(shù)

表5 縱蕩運(yùn)動(dòng)CFY展開系數(shù)

表4和表5為振蕩運(yùn)動(dòng)下側(cè)向力的展開系數(shù)，規(guī)律和表2、表3類似。隨著展長(zhǎng)增大，振蕩葉輪的側(cè)向力展開項(xiàng)中附加質(zhì)量力所占比例較小，且另外兩項(xiàng)系數(shù)絕對(duì)值隨著展長(zhǎng)增大而增大，逐漸接近二維結(jié)果。

4 結(jié)論

利用流體計(jì)算軟件ANSYS-CFX計(jì)算了不同展弦比的葉輪三維振蕩運(yùn)動(dòng)，包括橫蕩和縱蕩運(yùn)動(dòng)，并將計(jì)算結(jié)果和二維模擬結(jié)果做了對(duì)比，結(jié)論如下：

1)立軸葉輪三維模擬結(jié)果中水動(dòng)力系數(shù)和二維結(jié)果存在一定差異，當(dāng)展弦比接近10時(shí)，三維結(jié)果和二維結(jié)果相差不大。

2)隨著展弦比增大，三維振蕩葉輪的輸出效率、推力和側(cè)向力系數(shù)逐漸增大。

3)三維振蕩葉輪的水動(dòng)力性能計(jì)算結(jié)果低于二維計(jì)算結(jié)果，最主要原因是三維葉片兩端流體從高壓區(qū)流向低壓區(qū)，造成能量損失。

4)三維振蕩葉輪的推力和側(cè)向力各項(xiàng)中，附加質(zhì)量力和均勻水動(dòng)力、阻尼力相比為小量。