新型雙體波浪能轉換裝置多自由度耦合水動力特性研究

顧興遠, 牛玉博, 鄭 陽, 何培杰, 陳啟卷

(1.武漢大學 動力與機械學院,武漢 430072;2.武漢大學 水力機械過渡過程教育部重點實驗室,武漢 430072)

在大多數情況下,波浪能的能量密度在離岸較遠的位置往往比靠近海岸的位置更高[1-2],但與此同時,離岸較遠也使得波浪能轉換裝置(wave energy converter, WEC)生產的電能難以輸送到陸上電網,且安裝成本也會因此而提高。如果將WEC與其它海上設備結合,使WEC生產的電能為海上設備供電,可以有效解決這兩個問題。因此本文將一種可主動共振的縱搖式WEC[3-4]與海洋浮標相結合,設計出一種新型雙體WEC。

目前大多數雙體WEC都是雙體垂蕩系統,即利用兩個物體之間的垂蕩運動來轉換波浪能[5]。Al等[6]通過CFD獲取淹沒體的黏性阻力系數,再通過動力學建模,研究垂蕩雙體WEC中水下物體的不同形狀的特性,但僅考慮了垂蕩運動,且鑒于雙體距離較遠而忽略了雙體之間的輻射力影響。Ma等[7]結合Fortran編程與AQWA軟件,分析了垂蕩雙體WEC的頻域水動力性能和時域運動特性,但忽略了其它自由度的運動。周亞輝等[8]對穿孔式雙浮體裝置開展頻域仿真,研究垂蕩雙體的特性,同樣忽略了垂蕩以外的運動。對垂蕩雙體WEC的建模研究,往往僅考慮垂蕩這單一自由度,而忽略其它自由度的影響。Dai等研究一種新型雙體WEC,主體為由繩索連接的兩個球體,同時對雙體的垂蕩、縱蕩、縱搖運動開展研究,但考慮到裝置為球體,忽略了黏性影響,此外還忽略了兩個物體運動產生的輻射力對其它自由度的影響,即未計及各自由度之間的耦合。Tran等[9]同時考慮了黏性以及各自由度之間耦合,對一種圓柱體式的單體WEC開展了水動力性能研究。對多自由度的雙體WEC的建模研究,往往通過簡化裝置外形等方式忽略黏性以及各自由度之間耦合的影響,從而避免模型過于復雜。

本文提出的新型雙體縱搖WEC的外形較為復雜,因此對該裝置的建模仿真中同時考慮了雙體各自由度耦合以及黏性影響。首先通過AQWA開展頻域仿真,計算裝置的水動力參數;通過理論推導,在同時考慮雙體多自由度耦合以及黏性影響的情況下建立時域模型,并將AQWA仿真獲取的水動力參數代入時域模型,分別開展規則波、不規則波的仿真,分析質心位置、PTO阻尼大小對裝置發電性能的影響,根據真實海況劃分若干不規則波并分別開展仿真,研究裝置在真實海況下的發電性能。

1 裝置設計

本文提出的新型波能裝置由兩部分組成:一部分是漂浮在水面上的圓盤形浮標,另一部分是重力與浮力相等的可調節質心的擺體,如圖1所示。

圖1 新型雙體WEC

支架上端與海洋浮標固定,下端的橫軸穿過擺體,使得擺體可繞橫軸轉動。擺體內部密封有配重、配重調節機構、PTO機構。配重位置調節機構可控制配重升降,從而調節擺體的質心位置。PTO機構可將擺體與橫軸之間相對轉動的機械能轉化為電能。裝置在工作時,浮標漂浮在水面,而經過設計的擺體受到的重力與浮力相等,完全淹沒在水下,因此可減少風浪沖擊和能量輻射,同時減小自由表面黏性帶來的影響。相比于將PTO機構安裝在海床上的單體WEC[10],新型雙體WEC更容易維護。

2 水動力建模

此前研究人員對雙體WEC的建模尚未見同時考慮雙體各自由度耦合以及黏性影響的情況,因此本節將詳細推導新型雙體縱搖WEC的水動力模型。

2.1 勢流理論

以垂蕩運動為例,波浪與物體相互作用的控制方程可表示為

(1)

式中:m為物體質量;g為重力加速度;p為作用在物體上的流體壓強;S為物體的濕表面;z為物體在垂蕩方向的位移,正方向向上;Fother為包括PTO力在內的其它作用力的合力。

為進一步表達等式右側各項,引入線性勢流理論,入射流體作用在裝置表面的總壓強pa可由Bernoulli方程計算為

(2)

式中:ρw為海水密度;zs為濕表面面元位置;φt為入射流體的總速度勢。

根據線性勢流理論,入射流體的總速度勢φt可表示為三部分

φt=φI+φdif+φrad

(3)

式中,φI,φdif和φrad分別為裝置不存在時入射波的速度勢、由于裝置存在導致的衍射波的速度勢和裝置運動導致的輻射波的速度勢。

將式(3)代入式(2),應用線性勢流理論忽略二次項,將各面元壓強沿濕表面積分,可得波浪-物體相互作用的控制方程

(4)

式中:Fe為波浪激勵力;Fr為輻射力;Fres為回復力;FPTO為PTO力。

2.2 運動方程

如圖2所示,設A為浮標與擺體的連接軸中心,B為浮標質心,C為擺體殼質心,D為配重質心,P為擺體整體的質心。

圖2 裝置簡圖

忽略橫蕩、艏搖、橫搖運動,只考慮裝置的縱蕩、垂蕩、縱搖運動,將這三種運動的方向分別記為Z、Y、RX方向,于是裝置的運動共四個自由度,分別是裝置整體的縱蕩、裝置整體的垂蕩、擺體的縱搖、浮標的縱搖。根據Lagrange力學理論,將yA、zA、θ1和θ2作為廣義坐標,其中yA、zA為連接軸中心位置,θ1為擺體的轉角,θ2為浮標轉角,以順時針為正方向。于是A點在慣性系中的坐標為

A=(yA,zA)

(5)

則P、B點在慣性系中的坐標可表示為

r1=p=(yA+L1sinθ1,zA+L1cosθ1)

(6)

r2=B=(yA+L2sinθ2,zA+L2cosθ2)

(7)

式中,L1、L2分別為AP、AB長度。根據Lagrange方程,系統的時域運動模型可表示為

(8)

式中:qj為廣義坐標;T為動能;Qj為廣義力。廣義力的定義為

(9)

式中:Fv為外力;rv為外力作用點的位置矢量。將外力作用點移至A點,合外力在Y方向、Z方向的分量分別為Fy、Fz,合外力作用在擺體、浮標上的扭矩分別為M1、M2。將式(9)代入式(8),對式(6)、式(7)求導得到速度,再根據動能定理,可得裝置四個自由度的運動方程

(10)

式中:m1,m2分別為擺體、浮標的質量;L1,L2分別為擺體、浮標的質心到轉動軸的距離;v1,v2分別為擺體、浮標質心的速度;J1,J2分別為擺體、浮標繞質心的轉動慣量;ω1,ω2分別為擺體、浮標的角速度。根據勢流理論,Fy、Fz、M1、M2分別可表示為

(11)

式中,MPTO為PTO力矩,各下標的含義為:1為力作用于擺體,2為力作用于浮標,y為力作用于Y方向,z為作用于Z方向,r為輻射力,e為波浪激勵力,b為浮力,d為黏性力。

波浪激勵力可表示為

(12)

式中:下標i(i=1, 2) 為受力物體;下標j(j=y,z) 為力的方向;k為波浪激勵力系數;β為波浪激勵力相角,即波浪激勵力與入射波的相位差。

浮力可表示為

(13)

式中:S為浮標底面積;ρw為海水密度;Lb1為擺體浮心與A點的距離;kb2為浮標的水靜力剛度。

根據Cummins方程,擺體受到的輻射力可表示為卷積形式

(14)

式中:K(t)為K(ω)的脈沖響應函數,K(ω)由附加質量和輻射阻尼計算得到;下標i(i=y,z) 為力的方向;下標11、21分別為擺體的運動對擺體自身的影響、浮標的運動對擺體的影響,以K11yi為例,表示擺體自身的y方向運動對i(i=y,z)方向造成的輻射力影響的脈沖響應函數。浮標受到的輻射力計算式與擺體類似,僅將下標11、21改為22、12即可得到。

2.3 輻射力卷積近似

時域中的輻射力計算式是以卷積的形式表示的,不便于計算,因此這里采用Realization方法[11],對卷積項進行近似,將卷積項轉化為便于計算的狀態空間模型

(15)

式中:xr∈n×1為輻射力的狀態向量,Ar∈n×n、Br∈n×1、Cr∈1×n為參數矩陣,n為輻射力的狀態空間模型的階數,階數越高則狀態空間模型對卷積項的近似效果越好,但是需要更多的計算資源。

2.4 黏性力模型

WEC一般尺寸較小,黏性損失較明顯,因此有必要將黏性力引入到模型中做修正,本文采用Morison方程計算黏性力

(16)

式中:Sw為裝置的特征面積;Cd為裝置的黏性阻力系數;v為參考點處流體質點的速度,參考點選擇為物體豎直中線上各自高度的一半位置處,因此浮標、擺體的參考點分別為水面下0.5 m、4.5 m。根據Airy波浪理論,對于規則波

ζ=Acosα

(17)

參考點處流體的流速v可表示為

(18)

式中:ω為波浪頻率;ht為參考點深度;k為波數,滿足色散關系式

(19)

對于不規則波,參考點處流體的流速為各頻率子波在該點的流速之和。

3 仿真與分析

3.1 仿真對象

以典型的10 m直徑海洋浮標為載體,設計裝置的具體尺寸、質量等參數如表1所示。

表1 裝置設計參數

如圖2所示,擺體外輪廓由上下兩個60°的圓弧與兩側的曲線構成,R1為上部圓弧的半徑,R2為下部圓弧的半徑[3],W為擺體沿橫軸方向的寬度。D為浮標的直徑,H2為浮標的厚度,Hw為浮標淹沒深度,m1為擺體含配重的總質量,m2為浮標的質量,mt為配重的質量,J1s為擺體無配重時的轉動慣量,J2為浮標的轉動慣量,Jt為配重的轉動慣量,轉動慣量均為縱搖方向的值,參考軸在物體各自質心處。配重位置可調節范圍為0～0.8 m。

3.2 水動力參數計算

本節通過ANSYS AQWA軟件開展頻域仿真,計算WEC的輻射阻尼、附加質量和波浪激勵力系數等水動力參數。

為避免網格密度對計算結果的影響,本文對于浮標-擺體結構建立三個密度不同的網格進行仿真,進行網格無關性檢驗,結果表明,網格節點數為1 015時的仿真結果與網格節點數為3 355的結果略有差別,而當網格節點數增大到11 823時,仿真結果與網格節點數為3 355時的基本相同,可認為此時網格節點數繼續增大也基本不會對計算結果的準確度造成影響,因此最終選用11 823個網格節點的仿真結果。

圖3為水動力參數的計算結果,其中kr為輻射阻尼。由于與輻射力相關的水動力參數(輻射阻尼、附加質量)包含大量雙體各自由度之間的耦合項,在圖3中只展示其中一部分。可以看出,擺體運動對浮標所受的輻射力的影響較小,明顯小于浮標對擺體的影響。

圖3 雙體縱蕩運動在RX方向產生的輻射阻尼

部分輻射力參數(輻射阻尼、附加質量)的值基本為0,具體包括:縱蕩產生的輻射力參數在Z方向的分量,垂蕩產生的輻射力參數在Y、RX方向的分量,以及縱搖產生的輻射力參數在Z方向的分量。根據Falnes的理論,對于對稱物體,根據其對稱平面的特征,有部分的輻射力參數基本為0[12],而擺體與浮標均有x=0以及y=0這兩個對稱面,因此存在大量等于0的輻射力參數。忽略部分參數,最終保留20個輻射力參數分別進行卷積的狀態空間模型近似。

3.3 黏性力參數計算

對于本文中浮標外形的扁圓柱體,縱蕩、垂蕩、縱搖方向的黏性阻力系數分別取為CdY2=1,CdZ2=1,CdRX2=0.2。由于擺體在縱蕩和垂蕩方向上近似橫置的扁圓柱體,同樣分別取CdY1=1,CdZ1=1。對于擺體的縱搖方向的黏性阻力系數,可以通過CFD的自由衰減振蕩仿真獲得[13]。

本文通過Fluent對擺體單體的自由衰減振蕩過程開展二維的CFD仿真,最終得到擺體縱搖方向的黏性阻力系數為CdRX2=0.213 3。圖4為三種模型(CFD、加入黏性項的非線性模型、無黏性的線性模型)的仿真結果,設定擺體的初始釋放轉角為50°。結果表明,加入黏性項后與CFD仿真結果擬合較好,明顯優于無黏性模型。同時也驗證了模型的可靠性,能夠以此模型開展后續仿真研究。

圖4 擺體自由衰減振蕩

3.4 規則波仿真

3.4.1 固有周期

為了研究在規則波作用下,裝置的運動與發電性能,將規則波代入前述建立的四自由度耦合非線性模型中開展仿真。

已有研究表明,對于以單個擺體作為俘能裝置的縱搖式WEC,通過改變其內部配重的質量分布,可以顯著影響其固有周期,通過控制固有周期與波浪周期一致,能夠使WEC與波浪共振,從而有效提高WEC的發電能力。對于本文的浮標與擺體結合的雙體WEC,有必要研究調節擺體內部配重的質量分布對裝置整體的發電性能所帶來的影響。

通過改變配重位置,可以改變擺體的質心位置、轉動慣量,從而改變擺體的固有頻率,擺體縱搖方向固有頻率的具體表達式為

(20)

式中:ωin,Keq,Jeq分別為擺體的固有頻率、等效剛度和等效慣量;Fb1為擺體的浮力;lb1為擺體浮心到連接軸中心的長度;Jadd1為擺體的附加質量,其值與擺體的固有頻率有關,可以通過迭代運算得到固有頻率。擺體的固有周期Tin可由固有頻率計算得到。擺體固有周期4～7 s對應配重位置范圍為0.070～0.595 m。

3.4.2 參數響應

圖5 規則波(T=5 s, H=1 m)下配重位置對裝置的影響

圖5的結果表明,隨著配重位置的上移,θ1m、θdm會先增后減,而θ2m變化較小。裝置發電功率的變化趨勢與θdm相似,二者均在配重位置接近0.314 m時達到峰值,可見調節擺體內的配重位置對裝置整體的發電功率影響顯著。

PTO阻尼對于振蕩體式WEC的發電能力的影響同樣顯著,對于單體、單自由度的WEC裝置,已有理論可以計算其最優PTO阻尼kopt,具體公式為

(21)

式中,kr為輻射阻尼;J為物體轉動慣量;Ja為附加質量;cPTO為PTO剛度;Fb為浮力;Lb為浮心到轉動軸的距離;G為重力;L為質心到轉動軸的距離;ω為波浪角頻率。根據式(21)計算可得,單自由度的擺體在5 s周期的波浪下的最優PTO阻尼為2.83×104N·m·s/rad。

而對于雙體WEC,目前未見有關最優PTO阻尼的計算方法,因此下面研究在規則波作用下,不同PTO阻尼對于裝置發電能力的影響。分別設定不同的PTO阻尼開展仿真,結果如圖6所示,其中kPTO表示PTO阻尼。

圖6 規則波(T=5 s, H=1 m)下PTO阻尼對裝置的影響

隨著PTO阻尼的增大,θ1m先快速減小,然后趨近于θ2m;而θ2m的變化幅度較小,基本穩定為8°;θdm則快速減小,逐漸趨近于2°。當PTO阻尼約為1×105N·m·s/rad時裝置發電功率達到峰值,明顯大于單自由度的擺體在5 s周期的波浪下的最優PTO阻尼。

3.5 不規則波仿真

真實海況的波浪并不規則波,而是由頻率、波高不同的子波組成的不規則波,因此有必要研究在不規則波作用下,裝置的運動與發電性能。

3.5.1 不規則波模型

本文采用標準JONSWAP譜[14]建立實際不規則波的模型,標準JONSWAP譜的表達式為

(22)

式中:Sω為波浪譜密度;Tp為峰值周期;Hs為有義波高,參數σ可表示為

(23)

根據波浪譜,即可以計算任意角頻率ω對應子波的振幅A(ω)

(24)

3.5.2 參數響應

由于不規則波與規則波的特性差別較大,有必要進一步研究在不規則波作用下配重位置對裝置發電性能的影響。將Hs=1 m,Tp=5 s 的標準JONSWAP譜的不規則波代入四自由度耦合非線性模型中,在初始配重位置不同的情況下,分別開展仿真,結果如圖7所示。豎直虛線為擺體固有周期與波浪周期一致時的配重位置(Ltr=0.314 m)。

圖7 不規則波(Tp=5 s, He=1 m)下配重位置對裝置影響

從圖7可以看出,裝置在不規則波的作用下,一些特性仍與在規則波的作用下相近:隨著配重位置的上移,θ1m、θdm會先增后減,而θ2m變化較小。裝置發電功率的變化趨勢與θdm相似,二者在相近的位置達到峰值,可見在不規則波的作用下調節擺體內的配重位置對裝置整體發電功率的影響依舊顯著,為了使裝置達到最大的發電功率,對擺體內配重位置的調節可以參考擺體固有周期與波浪周期一致時的配重位置。

而與規則波的仿真結果不同的是:θ1m、θdm在峰值附近變化較為平緩。浮標的轉角平均幅值較小,約為5°,而θ1m、θdm均較大,可見在不規則波作用下,浮標十分穩定,不會出現大幅縱搖。

同樣地,為研究PTO阻尼變化的影響,分別設定不同的PTO阻尼開展仿真,結果如圖8所示。

圖8 不規則波(Tp=5 s, He=1 m)下PTO阻尼對平均功率影響

與規則波仿真結果相似的是,隨著PTO阻尼的增大,θ1m先快速減小,然后減小幅度逐漸平緩,而θ2m基本穩定在6°。

與規則波仿真結果不同的是,θ1m與θdm相近,當PTO阻尼約為2×105N·m·s/rad 時,裝置發電功率達到峰值,且發電功率曲線在峰值附近的變化較為平緩。

3.5.3 真實海況仿真

進一步地,為研究裝置在實際海況的發電性能,選取中國東海、南海開展真實海況仿真[15-17],研究裝置在不同海域中的性能。

表2 中國東海海況下的最大時均發電功率

表3 中國南海海況下的最大時均發電功率

俘獲寬度比(capture width ratio, CWR)常用于衡量WEC的發電能力,其計算式為

(25)

(26)

式中:ρw為海水密度,取1 025 kg/m3; g為重力加速度;Te為不規則波的能量周期,對于本文采用的JONSWAP標準不規則波譜,基本滿足Te=0.9Tp的關系。

計算不同波況下裝置共振時的俘獲寬度比,如圖9所示。結果表明,裝置共振時的CWR普遍較高,能夠較好地俘獲波浪能。波高較小時,裝置的CWR較高,這表明裝置在小波高情況下可以有效利用波浪能量,從而保持一定的發電能力。

圖9 俘獲寬度比

4 結 論

本文提出了一種將擺體與海洋浮標結合的新型波浪能轉換裝置,通過AQWA建立頻域模型,計算裝置的水動力參數,并通過理論推導,建立非線性的時域模型,分別開展規則波、不規則波的仿真,對該裝置進行了研究。

通過AQWA開展頻域仿真,獲得裝置的水動力參數,對裝置建立四個自由度耦合的時域模型,并在模型中加入非線性黏性項,通過CFD仿真擺體自由衰減振蕩,獲得了相應的黏性阻力系數,并驗證了模型的準確性。

通過多自由度耦合的非線性模型開展時域仿真,研究裝置在規則波、不規則波作用下運行時,配重位置以及PTO阻尼對裝置的運動狀態以及發電功率的影響。結果表明:在不規則波作用下,裝置的最優PTO阻尼大于在規則波作用下的最優PTO阻尼,且明顯大于單個擺體的最優PTO阻尼;無論是在規則波還是不規則波作用下,調節配重位置都可以使裝置發電功率達到峰值,且調節時可以參考擺體固有周期與波浪周期一致時對應的配重位置。

最后,研究裝置在真實海況中運行的性能,分別將中國東海、中國南海的海況劃分為有義波高、峰值周期不同的不規則波矩陣,并分別優化配重位置,得到裝置在各海況的最大功率矩陣。結果表明:裝置在周期為4～7 s的波況下發電能力較好,波高為1 m時的平均CWR可達0.54;且裝置在小波高情況下也可以有效利用波浪能量,從而保持一定的發電能力。因此裝置較為適用于中國東海、中國南海海域。

本文對新型波浪能轉換裝置的研究可以為海洋浮標等海上設備的供電方式提供參考,并且本文建立的雙體四自由度耦合的非線性水動力模型可以為研究人員研究復雜外形的雙體WEC時提供參考。在將來的研究中還可以進一步探索系泊系統對裝置發電性能的影響。