垂蕩式波浪能裝置運動模型數值分析

丁 松，韓端鋒

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

垂蕩式波浪能裝置運動模型數值分析

丁 松，韓端鋒

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

垂蕩式波浪能裝置在海洋可再生能源開發中被廣泛應用，通過浮子與搖桿在垂蕩方向的相對運動吸收波浪能。在以往相關的運動預報數值分析中，通常基于微幅波假設，僅考慮浮子與搖桿在垂蕩方向的運動，忽略搖桿其他自由度運動。建立并求解了垂蕩式波浪能裝置的非線性聯合運動方程組，分析垂蕩式波浪能裝置的波浪載荷、浮子與平臺連接處的受力情況。數值計算系統的運動響應，并將計算結果與已有的試驗數據進行比較驗證，結果表明數值模擬的垂蕩式波能裝置的運動響應與試驗結果相符合。最后，應用本計算方法分析PTO(power take-off)參數對波能裝置發電性能的影響。

波能浪裝置；非線性運動方程；波能吸收系統；數值模擬

Abstract: Heaving Buoys Wave Energy Converter (WEC) is widely used in ocean energy exploitation. Most of the existing numerical simulations are based on the small amplitude assumption and the heave motion of the platform is ignored. In this paper, non-linear equations for wave loading and motion are developed and solved. The force between floater and platform is analyzed. A limited validation is carried out with experimental data. Then this numerical method is used to study the effect of WEC’s power take-off system on its power performance.

Keywords: wave energy converter; non-linear motion equation; power take-off system; numerical simulation

波浪能具有可再生、清潔的特點，相較于風能(0.4～0.6 kW/m2)和太陽能(0.1～0.2 kW/m2)，波浪能具有更大的能量密度(2～3 kW/m2)。很多年以前，就有人開始對波浪能裝置做過研究，經過長時間的模型和理論研究，到1975年左右，關于波浪能裝置的研究開始逐步進入實用階段，例如利用波浪能轉換裝置將波浪能量轉換為裝備的機械能，再通過能量轉換裝置將裝備的機械能轉換為電能等能量。目前，各種類型的波浪能裝置以達百種，可根據波浪能裝置的工作原理、設備尺寸以及作業海域對其進行分類[1-2]。

如圖1中(a)、(b)、(c)所示，可將波浪能裝置按照工作原理分成三類。圖1(a)為振蕩水柱式波浪能裝置，該裝置以氣室作為能量轉換機構，氣室的下部在水下開口，氣室內的海水與周邊海水聯通，上部的開口通過氣輪發電機與外部大氣聯通，通過氣室內氣壓的變化推動汽輪機轉動進而發電；圖1(b)為振蕩浮體式波浪能裝置，通過波浪與浮體相互作用吸收波浪能，再將自身機械能轉換為電能；圖1(c)為收縮波道式波浪能裝置，通過將波浪能轉換成水體的勢能，再通過水輪發電機發電。

如圖1中(d)、(e)、(f)所示，可將波浪能裝置按照設備尺寸及作業方向分成三類。圖1(d)為筏式波浪能裝置，一般由鉸接的筏體和液壓系統組成，筏式裝置順浪向布置，波浪作用下，筏體隨之運動，將波浪能轉換為筏體運動的機械能，再通過液壓系統轉化為電能；圖1(e)為橫式波浪能裝置，與筏式波浪能裝置類似，其布置方向垂直于來浪方向；圖1(f)為點吸式波浪能裝置，該裝置尺寸相對較小，可吸收多個自由度運動的機械能。

圖1 各類型波浪能裝置示意Fig. 1 Subcategories of Wave Energy Converters

如圖1(g)所示，可將波浪能裝置按照作業水域的位置分成三類。岸邊式波浪能裝置只能吸收岸邊區域的波浪能，這部分波浪能由于海底摩擦等原因波浪能量密度較低；近岸式波浪能裝置處于近岸區域，這一區域的波浪能量密度相對岸邊區域豐富；離岸式波浪能裝置吸收遠海的波浪能，遠海區域波浪能資源豐富，但是裝備容易受到極度高能波的影響。

圖2 垂蕩式波浪能裝置簡圖Fig. 2 Simplified sketch for Heaving Buoys

在上述的眾多波浪能裝置中，目前沒有任何一種波浪能裝置能夠達到商業用途并且性能明顯優于其他種類的波浪能裝置。不過Lopez等[1]提到，點吸式波浪能裝置具有結構簡單、性價比高的特點。在近期研制的波浪能裝置中，點吸式波浪能裝置的開發數量明顯高于其他類型的波浪能裝置。本文的研究對象是垂蕩式波浪能裝置，如圖2所示，在波浪的作用下，物體做垂蕩運動，把做垂蕩運動的物體稱為浮子，浮子沿搖桿軸線做垂蕩運動。浮子將波浪能轉換成浮子自身的機械能，再通過能量轉換裝置將浮子的機械能轉換為電能。

垂蕩式波浪能裝置的實例較多[3-6]，在以往相關的運動預報數值分析中，通常基于微幅波假設，僅考慮浮子與搖桿的垂蕩運動，忽略搖桿其他自由度運動(例如縱搖)，本文建立并求解了垂蕩式波能裝置的非線性聯合運動方程組，對系統的波浪載荷、浮子與平臺連接處的受力進行分析，計算系統的運動響應，并將數值計算的結果與已有的試驗數據進行比較驗證，結果表明數值模擬的垂蕩式波能裝置的運動響應與試驗結果相符合。最后，應用該數值方法探討了PTO阻尼系數BPTO對浮子發電性能的影響。

1 聯合運動方程

在船舶與海洋工程計算當中，通常將線性化的運動方程應用于小幅度的運動[7]。然而，當物體面臨陡波，運動響應較大，此時線性化的運動方程并不能滿足計算精度的需求。在計算垂蕩式波浪能裝置的運動時，將浮子和搖桿假設為剛體，不考慮受力后的變形。因而，系統的運動可用剛體的完全非線性運動方程來描述。

圖3 描述垂蕩式波浪能裝置運動的坐標系Fig. 3 Three coordinate systems used in motion simulation

1.1 坐標系建立

在描述垂蕩式波浪能裝置運動的過程中，共建立三個坐標系，如圖3所示。

絕對坐標系OXYZ：其中O為絕對坐標系的原點，可取空間中的任意一固定點，OZ軸垂直向上，OX軸指向右側，OX，OY，OZ符合右手定則；

隨體平動坐標系O′X′Y′Z′：它是一個輔助坐標系，只隨搖桿做平移運動而不做旋轉運動；

固體坐標系oxyz：整個坐標系被固定在搖桿上，其坐標原點o取在搖桿的重心，oz軸沿搖桿軸向中心線垂直向上。

在垂蕩式波浪能裝置的運動中，浮子僅能沿搖桿軸向運動，也就是說浮子與搖桿具有相同的旋轉運動。通過隨體平動坐標系O′X′Y′Z′與固體坐標系oxyz之間的歐拉角(α,β,γ)來描述旋轉運動。因此，絕對坐標系OXYZ與固體坐標系oxyz之間的坐標轉換關系：

式中：T為隨體平動坐標系O′X′Y′Z′與固體坐標系oxyz之間的坐標轉換關系矩陣。

其中：(x,y,z)為固體坐標系oxyz中的坐標值；(X0,Y0,Z0)為隨體平動坐標系O′X′Y′Z′原點O′在絕對坐標系OXYZ中的坐標值；(X,Y,Z)為絕對坐標系OXYZ中的坐標值。在之后計算中，將涉及浮子的相對垂蕩位移、速度、加速度由固體坐標系oxyz到絕對坐標系OXYZ的轉換，在這里給出任意向量L的時間導數的表達式：

式中：Ω×xb為牽連速度，固體坐標系oxyz繞其原點o以角速度Ω轉動并帶動浮子一起轉動而引起的速度。

1.2 聯合運動方程的推導

搖桿和浮子在計算中被假設為剛體，搖桿的運動姿態可根據其質心的位移X0和速度U0、隨體平動坐標系O′X′Y′Z′相對于固體坐標系oxyz之間的歐拉角(α,β,γ)以及角速度Ω來描述；由于浮子僅能沿搖桿軸線上下運動，浮子和搖桿具有相同的旋轉運動。因此，浮子的運動姿態僅需固體坐標系oxyz下浮子的位移xb、速度ub便可以確定。基于1.1小節中對浮子速度、加速度由固體坐標系oxyz到絕對坐標系OXYZ的推導，則一般形式的垂蕩式波浪能裝置的聯合運動方程：

式中：X0，U0為搖桿質心位移、速度；θ，Ω為搖桿歐拉角位移、角速度；xb，ub為浮子相對搖桿的垂蕩位移、速度；FS，NS為搖桿的合外力、合外力矩；FW為浮子的合外力；MS，IS為搖桿的質量矩陣、轉動慣量矩陣；MW為浮子的質量矩陣；B是由歐拉角(α,β,γ)組成的矩陣。

為便于計算，搖桿的平動運動方程基于絕對坐標系OXYZ，而搖桿的轉動運動方程、浮子的平動運動方程基于固體坐標系oxyz。因此，MS和IS中所有非對角線元素均為零，如質量分布為軸對稱，則Ixx=Iyy。另外，如系統無初始角速度，基于細長體理論的假設條件，可以得到Ωzb=0，因此Ω×ISΩ=0，此時搖桿的運動僅有5個自由度。

1.3 聯合運動方程的求解

垂蕩式波浪能裝置的運動方程是非線性方程，方程中各運動參數是相互耦合的，如歐拉角矩陣B，很難單獨求解其中的某個運動參數。此外，聯合方程右側的各項載荷是隨搖桿、浮子的位移、速度等運動參數而實時變化的。在通常情況下，常利用數值計算的方法求解垂蕩式波浪能裝置的運動響應。在此之前，需要將式(7)整理成一階微分方程組。下面將給出聯合運動方程分量表達式的推導。此處，先將FS，NS，FW分解。

式(9)～(11)基于絕對坐標系OXYZ，將式(9)～(11)代入聯合運動方程式(7)中的后三式，有：

由于在固體坐標系oxyz中，浮子僅有z軸方向運動。因此，僅提出式(14)分量表達式的第三行并進行整理：

將式(16)～(18)聯立，可表達成Ma=F的矩陣形式，其中M為7×7的“質量矩陣”，是聯合運動方程左側的各加速度項前的系數；a為7×1向量，其分量分別為垂蕩式波浪能裝置的7個自由度上的運動加速度(其中前6項為搖桿在絕對坐標系OXYZ中的六自由度運動、最后一項為浮子在固體坐標系oxyz中z軸垂蕩運動)；F為7×1向量，為各方程右側與加速度無關的各項組成的“合外力”。

至此，已將一般形式的垂蕩式波浪能裝置聯合運動方程式(7)整理成Ma=F的矩陣公式，如已知系統t時刻各自由度運動的位移、速度，則該時刻對應的M、F為已知，利用數值解法便可以求得該時刻的a，再利用龍格庫塔法等數值方法進行求解，得到下一時刻系統各自由度運動的位移、速度。

2 垂蕩式波浪能裝置波浪載荷分析

2.1搖桿波浪力FWAVE-S的計算

搖桿受到的水動力可分為非黏性力和黏性力兩部分，非黏性力應用Rainey提出的完全細長體公式進行計算[8-9]，將非黏性水動力分成三部分：單位浸沒長度力F1、底端力F2和自由表面力F3；黏性力Fdrag則應用莫里森公式中的阻力項進行計算。

2.1.1 單位浸沒長度力F1

式(22)為單位浸沒長度力，將單位長度浸沒力沿搖桿水下垂直方向積分得到該項波浪力作用于搖桿的合力。

2.1.2 底端力F2

底端力F2作用于搖桿底部，其中壓力項p=ph+ps，為水動壓力與靜水壓力之和，靜水壓力項psSl與式(11)中ρSgn一起，提供搖桿軸向的浮力和靜水回復力。式(25)中后兩項與Munk力矩相關，與波幅的平方成比例。

2.1.3 自由表面力F3

自由表面力F3作用于搖桿與波浪瞬時濕表面的相交處：

式中：α為搖桿中心線與波浪自由表面法向量之間的夾角，它的取值由搖桿的運動參數和波浪瞬時自由表面的空間位置共同決定;t為垂直于搖桿中心線的單位切向量，方向指向流體外側。從式(26)中不難看出，自由表面力與波幅的三次方成比例。

2.1.4 黏性水動力Fdrag

單位長度搖桿的黏性力Fdrag采用莫里森公式中的拖曳力項計算，CD為拖曳力系數，wn為相對速度w在與搖桿軸線垂直的平面內分量，R為搖桿截面半徑。

2.2 浮子波浪力FWAVE-W的計算

在各項波浪力的表達上，浮子所受波浪力FWAVE-W與搖桿所受波浪力FWAVE-S的計算方法是一致的，此處將不再贅述。值得注意的是，雖然FWAVE-W與FWAVE-S中各項波浪力的表達式相同，但由于浮子在固體坐標系中有相對運動，而搖桿與固體坐標系之間沒有相對運動，因此浮子的運動速度、加速度的表達式是和搖桿不同的，這點區別可以在展開的分量表達式中體現出來。此處特別再次給出搖桿和浮子上任意一點的速度、加速度計算式，以便對比：

2.3 附加質量系數CM和拖曳系數CD的取值

應用Sesam軟件[10]，分別對搖桿、浮子的水下部分建模。在本文數值計算中，搖桿高度(120 m)遠大于浮子(4 m)，浮子的存在對搖桿整體水動力的影響不大；而在計算浮子的水動力模型中，為防止浮子內部(月亮井)水面共振，在內部水面處加入浮板，在一定程度上考慮了搖桿的存在對浮子水動力性能的影響。

如圖4和圖5所示。在得到附加質量曲線后，根據入射波波浪周期，即可得到該工況下的各附加質量系數。例如，在之后的對比計算中，入射波波浪周期為12 s，計算得到搖桿水平方向附加質量系數為0.8，豎直方向附加質量系數為0.06；浮子水平方向附加質量系數為0.8，豎直方向附加質量系數為1.5。拖曳力系數的選取參考[8]中的取值，搖桿的拖曳系數為0.6，浮子的拖曳系數為1.0。

3 垂蕩式波浪能裝置連接阻尼分析

在垂蕩式波浪能裝置中，搖桿與浮子之間的連接裝置是非常重要的部分，它的作用是使浮子能且僅能沿搖桿軸向自由運動，而其他自由度的運動保持一致，顯然FC-S，FC-W是兩個等值反向的矢量。

在下面的表述中以浮子受到的連接力FC-W為主，計算當中，將FC-W分成兩個部分：

其中，FT-W為搖桿與浮子之間的正壓力，表達式見(15)；FPTO-W為能量吸收裝置在浮子工作期間對其運動產生的阻尼力：

式中：BPTO、KPTO分別為能量吸收裝置的阻尼系數和回復剛度系數，xb、ub分別為浮子垂蕩位移、速度。根據FPTO-W和ub可以得到一段時間內的浮子平均輸出功率，也就是反應了瞬時輸出功率曲線的均值：

4 聯合運動方程的數值解法及驗證

4.1 數值計算步驟

本文的數值計算程序在MATLAB中編寫，數值計算步驟如下：

1)給出系統初始時刻的運動參數，即搖桿位移X0、速度U0、歐拉角(α0,β0,γ0)、角速度Ω；浮子在固體坐標系oxyz中位移xb、速度ub。本文計算中假定它們在初始時刻為靜止狀態。

2)由搖桿、浮子的運動姿態，計算它們的瞬時浸沒長度、相對流速以及導纜孔在絕對坐標系中的位置。

3)根據搖桿、浮子的瞬時浸沒長度等參數計算搖桿、浮子的合外力、合外力矩。

4)根據之前的推導，將計算好的各合外力、合外力矩代入聯合運動方程式(8)中，并整理為一階非線性常系數微分方程組，采用自適應階的四階龍格庫塔法求解搖桿、浮子的聯合運動方程，在每一個時間階內，系統的各項運動參數可被求解，通過迭代法達到所要求的誤差控制限，輸出并記錄搖桿、浮子的運動參數：搖桿位移X0、速度U0、歐拉角(α0,β0,γ0)、角速度Ω；浮子在固體坐標系oxyz中位移xb、速度ub。

5)返回步驟2)，進行下一時刻的計算。計算結束于最終計算時刻。

運動模型的數值計算流程圖如圖6所示。

圖6 運動模型數值計算流程圖Fig. 6 Flow diagram of numerical simulation

4.2 數值計算結果與試驗數據對比

本文的對比驗證針對STC型風浪能混合利用系統中的垂蕩式波浪能裝置[11]。需要指出的是，本數值計算中不考慮風載荷，且波浪輸入為微幅波，并不涉及極限海況下系統載荷的計算問題。該物理模型試驗在挪威MARINTEK拖曳水池進行[12]，Froude縮尺比為1∶50，該模型中浮子可沿搖桿軸向上下自由滑動，能量吸收裝置PTO不工作(阻尼力為零)。對比工況的波浪輸入為規則波，波浪周期1.697 s，波高0.08 m。各自由度運動時歷曲線的對比如圖7所示，對比數據均為無量綱化數據。從圖中可見，在該工況下，各自由度運動的數值計算結果與試驗數據相符合，驗證了本文關于垂蕩式波浪能裝置的模型建立、數值計算方法的正確性和適用性。

圖7 數值計算結果與試驗數據對比Fig. 7 The comparison between numerical calculation and experimental result

5 PTO阻尼系數BPTO參數分析

由式(35)(36)可知，根據FPTO-W和ub可以得到浮子的瞬時輸出功率曲線，也就得到了一段時間TC內的浮子平均輸出功率：

在數值計算中KPTO=0，浮子輸出功率的來源是第一項，因此，阻尼系數BPTO的取值對輸出功率具有很大的影響。本節應用提出的數值計算方法，對能量吸收系統PTO阻尼系數BPTO進行參數分析，計算不同PTO阻尼系數(1 000 kNs/m、2 000 kNs/m、4 000 kNs/m、8 000 kNs/m)下波浪能裝置的俘獲寬度比[13]。測試工況中，波浪頻率范圍0.3～2 rad/s，波高分別為2 m和4 m。數值計算中的設備主要參數如表1所示。

表1 數值計算中設備主要參數Tab. 1 Model dimensions used in numerical simulation

圖8 各PTO參數下波能裝置的輸出功率曲線Fig. 8 WEC’s power performance with different PTO coefficients

經過數值計算得到各工況下浮子俘獲寬度比曲線如圖8所示，通過分析圖8可以得到如下結論：

由于數值計算中的波浪均為一階線性波，當BPTO一定時，波高對俘獲寬度比幾乎沒有影響。例如，當BPTO為4 000 kNs/m時，波高2 m和波高4 m的俘獲寬度比曲線峰值均出現在波浪頻率為0.8 rad/s左右，且吸波效率均為30.6%。

當BPTO增大時，俘獲寬度比曲線的峰值增大，例如，波高2 m時，俘獲寬度比曲線的峰值由BPTO為1 000 kNs/m時的35.9%增長到PTO阻尼系數為8 000 kNs/m時的57.6%；而隨著BPTO的增大，峰值所對應的俘獲寬度比減小，峰值寬度也隨之減小，例如，當BPTO為1 000 kNs/m時，俘獲寬度比曲線峰值對應的入射波頻率為1.3 rad/s，設備對0.9～1.7 rad/s波浪頻率區間的入射波都有良好的吸收效率(俘獲寬度比大于15%)，而當BPTO增長到8 000 kNs/m時，俘獲寬度比曲線峰值對應的入射波頻率減小到0.7 rad/s，此時設備僅對0.7～0.9 rad/s波浪頻率區間的入射波有良好的吸收效率。

6 結 語

對垂蕩式波浪能裝置的運動響應進行了數值分析，在剛體假設條件下建立包含搖桿、浮子在內的垂蕩式波浪能裝置的非線性聯合運動方程。分析垂蕩式波浪能裝置的波浪載荷、浮子與平臺連接處的受力情況。編寫MATLAB程序，運用自適應階的四階龍格庫塔法求解聯合運動方程組，模擬垂蕩式波浪能裝置的運動響應。將系統運動的數值計算結果與已有的試驗數據進行驗證，數值計算結果與試驗數據相符合，驗證了本文關于垂蕩式波浪能裝置的模型建立、數值計算方法的正確性和適用性。在對比工況中，搖桿縱搖運動RAO約為1.5 rad-1，在本文提出的聯合運動方程中，搖桿縱搖運動與浮子垂蕩運動是存在相互耦合作用的，因此在垂蕩式波浪能裝置的數學模型中，有必要考慮搖桿在縱搖等其他自由度的運動。

最后，應用該數值方法探討了PTO阻尼系數BPTO對吸波浮子發電性能的影響：在一定范圍內，增加BPTO將使得波能俘獲寬度比曲線峰值“變尖”，且峰值對應的入射波頻率向低頻區轉移。通過分析阻尼系數BPTO的取值對俘獲寬度比的影響，為之后PTO系統的設計提供依據。

[1] LóPEZ I, ANDREU J, CEBALLOS S, et al. Review of wave energy technologies and the necessary power-equipment[J]. Renewable and Sustainable Energy Reviews, 2013, 27: 413-434.

[2] ANTONIO F O. Wave energy utilization: A review of the technologies[J]. Renewable and Sustainable Energy Reviews, 2010, 14(3): 899-918.

[3] ELWOOD D, SCHACHER A, RHINEFRANK K, et al. Numerical modeling and ocean testing of a direct-drive wave energy device utilizing a permanent magnet linear generator for power take-off[C]//Proceedings of the 28th International Conference on Ocean, Offshore and Arctic Engineering. American Society of Mechanical Engineers, 2009: 817-824.

[4] FERDINANDE V, VANTORRE M. The concept of a bipartite point absorber[M]//Hydrodynamics of Ocean Wave-Energy Utilization. Springer Berlin Heidelberg, 1986: 217-226.

[5] NOREN S A. Plant for utilizing kinetic energy: U.S. Patent 4,277,690[P]. 1981-07-07.

[6] WEBER J, MOUWEN F, PARISH A, et al. Wavebob—research & development network and tools in the context of systems engineering[C]//Proceedings of the Eighth European Wave and Tidal Energy Conference. 2009.

[7] MEI C C. The applied dynamics of ocean surface waves[M]. World Scientific, 1989.

[8] RAINEY R C T. Slender-body expressions for the wave load on offshore structures[C]//Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, 1995, 450(1939): 391-416.

[9] MA Q W, PATEL M H. On the non-linear forces acting on a floating spar platform in ocean waves[J]. Applied Ocean Research, 2001, 23(1): 29-40.

[10] Norske Veritas D. SESAM-User Manual[M]. Hovik, Norway, 1998.

[11] MULIAWAN M J, KARIMIRAD M, MOAN T. Dynamic response and power performance of a combined spar-type floating wind turbine and coaxial floating wave energy converter[J]. Renewable Energy, 2013, 50: 47-57.

[12] WAN L, GAO Z, MOAN T. Model test of the STC concept in survival modes[C]//Proceedings of the 33rd International Conference on Ocean, Offshore and Arctic Engineering. American Society of Mechanical Engineers, 2014: V09AT09A010-V09AT09A010.

[13] FALNES J. Ocean waves and oscillating systems: linear interactions including wave-energy extraction[M]. Cambridge:Cambridge University Press, 2002.

Numerical modeling of heaving buoys wave energy converter

DING Song， HAN Duanfeng

(College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China)

1005-9865(2016)04-0107-11

P743.2

A

10.16483/j.issn.1005-9865.2016.04.015

2015-11-03

國家自然科學基金資助項目(51379051)

丁 松(1989-), 男, 博士研究生，主要從事海洋結構物水動力性能分析研究。E-mail: dingsong2014@hotmail.com