丁 松，韓端鋒

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

垂蕩式波浪能裝置運(yùn)動(dòng)模型數(shù)值分析

丁 松，韓端鋒

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

垂蕩式波浪能裝置在海洋可再生能源開(kāi)發(fā)中被廣泛應(yīng)用，通過(guò)浮子與搖桿在垂蕩方向的相對(duì)運(yùn)動(dòng)吸收波浪能。在以往相關(guān)的運(yùn)動(dòng)預(yù)報(bào)數(shù)值分析中，通?；谖⒎僭O(shè)，僅考慮浮子與搖桿在垂蕩方向的運(yùn)動(dòng)，忽略搖桿其他自由度運(yùn)動(dòng)。建立并求解了垂蕩式波浪能裝置的非線性聯(lián)合運(yùn)動(dòng)方程組，分析垂蕩式波浪能裝置的波浪載荷、浮子與平臺(tái)連接處的受力情況。數(shù)值計(jì)算系統(tǒng)的運(yùn)動(dòng)響應(yīng)，并將計(jì)算結(jié)果與已有的試驗(yàn)數(shù)據(jù)進(jìn)行比較驗(yàn)證，結(jié)果表明數(shù)值模擬的垂蕩式波能裝置的運(yùn)動(dòng)響應(yīng)與試驗(yàn)結(jié)果相符合。最后，應(yīng)用本計(jì)算方法分析PTO(power take-off)參數(shù)對(duì)波能裝置發(fā)電性能的影響。

波能浪裝置；非線性運(yùn)動(dòng)方程；波能吸收系統(tǒng)；數(shù)值模擬

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

波浪能具有可再生、清潔的特點(diǎn)，相較于風(fēng)能(0.4～0.6 kW/m2)和太陽(yáng)能(0.1～0.2 kW/m2)，波浪能具有更大的能量密度(2～3 kW/m2)。很多年以前，就有人開(kāi)始對(duì)波浪能裝置做過(guò)研究，經(jīng)過(guò)長(zhǎng)時(shí)間的模型和理論研究，到1975年左右，關(guān)于波浪能裝置的研究開(kāi)始逐步進(jìn)入實(shí)用階段，例如利用波浪能轉(zhuǎn)換裝置將波浪能量轉(zhuǎn)換為裝備的機(jī)械能，再通過(guò)能量轉(zhuǎn)換裝置將裝備的機(jī)械能轉(zhuǎn)換為電能等能量。目前，各種類(lèi)型的波浪能裝置以達(dá)百種，可根據(jù)波浪能裝置的工作原理、設(shè)備尺寸以及作業(yè)海域?qū)ζ溥M(jìn)行分類(lèi)[1-2]。

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

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

圖1 各類(lèi)型波浪能裝置示意Fig. 1 Subcategories of Wave Energy Converters

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

圖2 垂蕩式波浪能裝置簡(jiǎn)圖Fig. 2 Simplified sketch for Heaving Buoys

在上述的眾多波浪能裝置中，目前沒(méi)有任何一種波浪能裝置能夠達(dá)到商業(yè)用途并且性能明顯優(yōu)于其他種類(lèi)的波浪能裝置。不過(guò)Lopez等[1]提到，點(diǎn)吸式波浪能裝置具有結(jié)構(gòu)簡(jiǎn)單、性?xún)r(jià)比高的特點(diǎn)。在近期研制的波浪能裝置中，點(diǎn)吸式波浪能裝置的開(kāi)發(fā)數(shù)量明顯高于其他類(lèi)型的波浪能裝置。本文的研究對(duì)象是垂蕩式波浪能裝置，如圖2所示，在波浪的作用下，物體做垂蕩運(yùn)動(dòng)，把做垂蕩運(yùn)動(dòng)的物體稱(chēng)為浮子，浮子沿?fù)u桿軸線做垂蕩運(yùn)動(dòng)。浮子將波浪能轉(zhuǎn)換成浮子自身的機(jī)械能，再通過(guò)能量轉(zhuǎn)換裝置將浮子的機(jī)械能轉(zhuǎn)換為電能。

垂蕩式波浪能裝置的實(shí)例較多[3-6]，在以往相關(guān)的運(yùn)動(dòng)預(yù)報(bào)數(shù)值分析中，通?；谖⒎僭O(shè)，僅考慮浮子與搖桿的垂蕩運(yùn)動(dòng)，忽略搖桿其他自由度運(yùn)動(dòng)(例如縱搖)，本文建立并求解了垂蕩式波能裝置的非線性聯(lián)合運(yùn)動(dòng)方程組，對(duì)系統(tǒng)的波浪載荷、浮子與平臺(tái)連接處的受力進(jìn)行分析，計(jì)算系統(tǒng)的運(yùn)動(dòng)響應(yīng)，并將數(shù)值計(jì)算的結(jié)果與已有的試驗(yàn)數(shù)據(jù)進(jìn)行比較驗(yàn)證，結(jié)果表明數(shù)值模擬的垂蕩式波能裝置的運(yùn)動(dòng)響應(yīng)與試驗(yàn)結(jié)果相符合。最后，應(yīng)用該數(shù)值方法探討了PTO阻尼系數(shù)BPTO對(duì)浮子發(fā)電性能的影響。

1 聯(lián)合運(yùn)動(dòng)方程

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

圖3 描述垂蕩式波浪能裝置運(yùn)動(dòng)的坐標(biāo)系Fig. 3 Three coordinate systems used in motion simulation

1.1 坐標(biāo)系建立

在描述垂蕩式波浪能裝置運(yùn)動(dòng)的過(guò)程中，共建立三個(gè)坐標(biāo)系，如圖3所示。

絕對(duì)坐標(biāo)系OXYZ：其中O為絕對(duì)坐標(biāo)系的原點(diǎn)，可取空間中的任意一固定點(diǎn)，OZ軸垂直向上，OX軸指向右側(cè)，OX，OY，OZ符合右手定則；

隨體平動(dòng)坐標(biāo)系O′X′Y′Z′：它是一個(gè)輔助坐標(biāo)系，只隨搖桿做平移運(yùn)動(dòng)而不做旋轉(zhuǎn)運(yùn)動(dòng)；

固體坐標(biāo)系oxyz：整個(gè)坐標(biāo)系被固定在搖桿上，其坐標(biāo)原點(diǎn)o取在搖桿的重心，oz軸沿?fù)u桿軸向中心線垂直向上。

在垂蕩式波浪能裝置的運(yùn)動(dòng)中，浮子僅能沿?fù)u桿軸向運(yùn)動(dòng)，也就是說(shuō)浮子與搖桿具有相同的旋轉(zhuǎn)運(yùn)動(dòng)。通過(guò)隨體平動(dòng)坐標(biāo)系O′X′Y′Z′與固體坐標(biāo)系oxyz之間的歐拉角(α,β,γ)來(lái)描述旋轉(zhuǎn)運(yùn)動(dòng)。因此，絕對(duì)坐標(biāo)系OXYZ與固體坐標(biāo)系oxyz之間的坐標(biāo)轉(zhuǎn)換關(guān)系：

式中：T為隨體平動(dòng)坐標(biāo)系O′X′Y′Z′與固體坐標(biāo)系oxyz之間的坐標(biāo)轉(zhuǎn)換關(guān)系矩陣。

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

式中：Ω×xb為牽連速度，固體坐標(biāo)系oxyz繞其原點(diǎn)o以角速度Ω轉(zhuǎn)動(dòng)并帶動(dòng)浮子一起轉(zhuǎn)動(dòng)而引起的速度。

1.2 聯(lián)合運(yùn)動(dòng)方程的推導(dǎo)

搖桿和浮子在計(jì)算中被假設(shè)為剛體，搖桿的運(yùn)動(dòng)姿態(tài)可根據(jù)其質(zhì)心的位移X0和速度U0、隨體平動(dòng)坐標(biāo)系O′X′Y′Z′相對(duì)于固體坐標(biāo)系oxyz之間的歐拉角(α,β,γ)以及角速度Ω來(lái)描述；由于浮子僅能沿?fù)u桿軸線上下運(yùn)動(dòng)，浮子和搖桿具有相同的旋轉(zhuǎn)運(yùn)動(dòng)。因此，浮子的運(yùn)動(dòng)姿態(tài)僅需固體坐標(biāo)系oxyz下浮子的位移xb、速度ub便可以確定?；?.1小節(jié)中對(duì)浮子速度、加速度由固體坐標(biāo)系oxyz到絕對(duì)坐標(biāo)系OXYZ的推導(dǎo)，則一般形式的垂蕩式波浪能裝置的聯(lián)合運(yùn)動(dòng)方程：

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

為便于計(jì)算，搖桿的平動(dòng)運(yùn)動(dòng)方程基于絕對(duì)坐標(biāo)系OXYZ，而搖桿的轉(zhuǎn)動(dòng)運(yùn)動(dòng)方程、浮子的平動(dòng)運(yùn)動(dòng)方程基于固體坐標(biāo)系oxyz。因此，MS和IS中所有非對(duì)角線元素均為零，如質(zhì)量分布為軸對(duì)稱(chēng)，則Ixx=Iyy。另外，如系統(tǒng)無(wú)初始角速度，基于細(xì)長(zhǎng)體理論的假設(shè)條件，可以得到Ωzb=0，因此Ω×ISΩ=0，此時(shí)搖桿的運(yùn)動(dòng)僅有5個(gè)自由度。

1.3 聯(lián)合運(yùn)動(dòng)方程的求解

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

式(9)～(11)基于絕對(duì)坐標(biāo)系OXYZ，將式(9)～(11)代入聯(lián)合運(yùn)動(dòng)方程式(7)中的后三式，有：

由于在固體坐標(biāo)系oxyz中，浮子僅有z軸方向運(yùn)動(dòng)。因此，僅提出式(14)分量表達(dá)式的第三行并進(jìn)行整理：

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

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

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

2.1搖桿波浪力FWAVE-S的計(jì)算

搖桿受到的水動(dòng)力可分為非黏性力和黏性力兩部分，非黏性力應(yīng)用Rainey提出的完全細(xì)長(zhǎng)體公式進(jìn)行計(jì)算[8-9]，將非黏性水動(dòng)力分成三部分：?jiǎn)挝唤](méi)長(zhǎng)度力F1、底端力F2和自由表面力F3；黏性力Fdrag則應(yīng)用莫里森公式中的阻力項(xiàng)進(jìn)行計(jì)算。

2.1.1 單位浸沒(méi)長(zhǎng)度力F1

式(22)為單位浸沒(méi)長(zhǎng)度力，將單位長(zhǎng)度浸沒(méi)力沿?fù)u桿水下垂直方向積分得到該項(xiàng)波浪力作用于搖桿的合力。

2.1.2 底端力F2

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

2.1.3 自由表面力F3

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

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

2.1.4 黏性水動(dòng)力Fdrag

單位長(zhǎng)度搖桿的黏性力Fdrag采用莫里森公式中的拖曳力項(xiàng)計(jì)算，CD為拖曳力系數(shù)，wn為相對(duì)速度w在與搖桿軸線垂直的平面內(nèi)分量，R為搖桿截面半徑。

2.2 浮子波浪力FWAVE-W的計(jì)算

在各項(xiàng)波浪力的表達(dá)上，浮子所受波浪力FWAVE-W與搖桿所受波浪力FWAVE-S的計(jì)算方法是一致的，此處將不再贅述。值得注意的是，雖然FWAVE-W與FWAVE-S中各項(xiàng)波浪力的表達(dá)式相同，但由于浮子在固體坐標(biāo)系中有相對(duì)運(yùn)動(dòng)，而搖桿與固體坐標(biāo)系之間沒(méi)有相對(duì)運(yùn)動(dòng)，因此浮子的運(yùn)動(dòng)速度、加速度的表達(dá)式是和搖桿不同的，這點(diǎn)區(qū)別可以在展開(kāi)的分量表達(dá)式中體現(xiàn)出來(lái)。此處特別再次給出搖桿和浮子上任意一點(diǎn)的速度、加速度計(jì)算式，以便對(duì)比：

2.3 附加質(zhì)量系數(shù)CM和拖曳系數(shù)CD的取值

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

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

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

在垂蕩式波浪能裝置中，搖桿與浮子之間的連接裝置是非常重要的部分，它的作用是使浮子能且僅能沿?fù)u桿軸向自由運(yùn)動(dòng)，而其他自由度的運(yùn)動(dòng)保持一致，顯然FC-S，F(xiàn)C-W是兩個(gè)等值反向的矢量。

在下面的表述中以浮子受到的連接力FC-W為主，計(jì)算當(dāng)中，將FC-W分成兩個(gè)部分：

其中，F(xiàn)T-W為搖桿與浮子之間的正壓力，表達(dá)式見(jiàn)(15)；FPTO-W為能量吸收裝置在浮子工作期間對(duì)其運(yùn)動(dòng)產(chǎn)生的阻尼力：

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

4 聯(lián)合運(yùn)動(dòng)方程的數(shù)值解法及驗(yàn)證

4.1 數(shù)值計(jì)算步驟

本文的數(shù)值計(jì)算程序在MATLAB中編寫(xiě)，數(shù)值計(jì)算步驟如下：

1)給出系統(tǒng)初始時(shí)刻的運(yùn)動(dòng)參數(shù)，即搖桿位移X0、速度U0、歐拉角(α0,β0,γ0)、角速度Ω；浮子在固體坐標(biāo)系oxyz中位移xb、速度ub。本文計(jì)算中假定它們?cè)诔跏紩r(shí)刻為靜止?fàn)顟B(tài)。

2)由搖桿、浮子的運(yùn)動(dòng)姿態(tài)，計(jì)算它們的瞬時(shí)浸沒(méi)長(zhǎng)度、相對(duì)流速以及導(dǎo)纜孔在絕對(duì)坐標(biāo)系中的位置。

3)根據(jù)搖桿、浮子的瞬時(shí)浸沒(méi)長(zhǎng)度等參數(shù)計(jì)算搖桿、浮子的合外力、合外力矩。

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

5)返回步驟2)，進(jìn)行下一時(shí)刻的計(jì)算。計(jì)算結(jié)束于最終計(jì)算時(shí)刻。

運(yùn)動(dòng)模型的數(shù)值計(jì)算流程圖如圖6所示。

圖6 運(yùn)動(dòng)模型數(shù)值計(jì)算流程圖Fig. 6 Flow diagram of numerical simulation

4.2 數(shù)值計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)對(duì)比

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

圖7 數(shù)值計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)對(duì)比Fig. 7 The comparison between numerical calculation and experimental result

5 PTO阻尼系數(shù)BPTO參數(shù)分析

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

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

表1 數(shù)值計(jì)算中設(shè)備主要參數(shù)Tab. 1 Model dimensions used in numerical simulation

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

經(jīng)過(guò)數(shù)值計(jì)算得到各工況下浮子俘獲寬度比曲線如圖8所示，通過(guò)分析圖8可以得到如下結(jié)論：

由于數(shù)值計(jì)算中的波浪均為一階線性波，當(dāng)BPTO一定時(shí)，波高對(duì)俘獲寬度比幾乎沒(méi)有影響。例如，當(dāng)BPTO為4 000 kNs/m時(shí)，波高2 m和波高4 m的俘獲寬度比曲線峰值均出現(xiàn)在波浪頻率為0.8 rad/s左右，且吸波效率均為30.6%。

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

6 結(jié) 語(yǔ)

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

最后，應(yīng)用該數(shù)值方法探討了PTO阻尼系數(shù)BPTO對(duì)吸波浮子發(fā)電性能的影響：在一定范圍內(nèi)，增加BPTO將使得波能俘獲寬度比曲線峰值“變尖”，且峰值對(duì)應(yīng)的入射波頻率向低頻區(qū)轉(zhuǎn)移。通過(guò)分析阻尼系數(shù)BPTO的取值對(duì)俘獲寬度比的影響，為之后PTO系統(tǒng)的設(shè)計(jì)提供依據(jù)。

[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

國(guó)家自然科學(xué)基金資助項(xiàng)目(51379051)

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