何廣華，楊凱博，欒政曉，張志剛，劉朝綱，荊芃霖

（1.哈爾濱工業(yè)大學(xué)（威海）海洋工程學(xué)院，山東威海264209；2.山東船舶技術(shù)研究院，山東威海264209；3.哈爾濱工業(yè)大學(xué)機(jī)械工程學(xué)院，哈爾濱150001）

0 引 言

由于能源消耗和環(huán)境污染等問(wèn)題，海洋可再生能源的開發(fā)和利用逐漸受到各國(guó)的關(guān)注[1]。其中波浪能以儲(chǔ)量大、無(wú)污染、可重復(fù)開發(fā)利用等優(yōu)點(diǎn)，成為國(guó)內(nèi)外海洋能開發(fā)利用研究的熱點(diǎn)。波浪能發(fā)電裝置根據(jù)其工作原理可分為振蕩體式、振蕩水柱式、聚波越浪式等[2]。其中，振蕩體式又可分為振蕩浮子式（點(diǎn)吸收式）、筏式和擺式。近年來(lái)，振蕩浮子式WEC 因其體積小、制造安裝方便、不受波浪方向影響等優(yōu)勢(shì)受到了廣泛的關(guān)注[3]。常見的振蕩浮子式WEC 的PTO 系統(tǒng)可分為機(jī)械式和液壓式，其中液壓式PTO 利用液壓缸和液壓泵等元件將液壓能轉(zhuǎn)化成電能，機(jī)械式PTO 利用齒輪齒條、滾珠絲杠[4-5]等機(jī)械裝置將機(jī)械能轉(zhuǎn)化成電能。

波浪能發(fā)電裝置的PTO 系統(tǒng)是波浪能發(fā)電的關(guān)鍵技術(shù)之一，不同研究者對(duì)PTO 系統(tǒng)進(jìn)行了不同程度的簡(jiǎn)化研究。Rico 等[6-8]利用數(shù)學(xué)公式計(jì)算了搖臂圓柱形浮子的運(yùn)動(dòng)學(xué)方程，分析了線性PTO 阻尼對(duì)浮子俘獲功率的影響；張亞群等[9]利用物模試驗(yàn)分析了振蕩浮子的水動(dòng)力性能，實(shí)驗(yàn)顯示在最優(yōu)PTO 負(fù)載下浮子可獲得最大俘獲寬度比；宋文杰等[10]利用實(shí)驗(yàn)分析了雙行程工作形式PTO 的輸出特性，但未對(duì)單行程工作形式進(jìn)行對(duì)比分析；馮輝等[11]基于Simulink和AMESim 搭建了聯(lián)合仿真模型，該模型包括了各液壓元件與發(fā)電機(jī)，較好地模擬出了能量轉(zhuǎn)換系統(tǒng)的發(fā)電特性曲線，但其未根據(jù)仿真結(jié)果對(duì)PTO 元件進(jìn)行優(yōu)化選型；Li等[12]利用CFD 工具研究了波物相互作用與機(jī)械元件對(duì)浮子的影響，發(fā)現(xiàn)機(jī)械元件產(chǎn)生的力和力矩對(duì)單個(gè)浮子的運(yùn)動(dòng)響應(yīng)影響巨大；劉延偉等[13]對(duì)超越離合器進(jìn)行了理想建模，并分析了在切換瞬間超越離合器的動(dòng)態(tài)特性。

然而，對(duì)于采用機(jī)械式PTO 的振蕩浮子波浪能發(fā)電裝置，PTO 系統(tǒng)的傳動(dòng)方式、發(fā)電機(jī)選型及電機(jī)的電子負(fù)載等參數(shù)都會(huì)對(duì)波浪能發(fā)電裝置的獲能產(chǎn)生影響，導(dǎo)致將PTO 負(fù)載簡(jiǎn)化為線性阻尼的方法不夠精確。針對(duì)該問(wèn)題，本文采用水動(dòng)力模型與PTO系統(tǒng)聯(lián)合仿真的方法，使所添加的負(fù)載更接近真實(shí)情況，其中，在STAR-CCM+中建立波浪能發(fā)電裝置數(shù)值模型，在AMESim 中建立PTO 系統(tǒng)仿真模型，通過(guò)二者的聯(lián)合仿真進(jìn)行水動(dòng)力與PTO系統(tǒng)的耦合研究；研究PTO負(fù)載的不同行程工作模式與不同傳動(dòng)比對(duì)浮子運(yùn)動(dòng)響應(yīng)和輸出功率的影響，并針對(duì)目標(biāo)海況對(duì)PTO 參數(shù)進(jìn)行優(yōu)化。本文使用的聯(lián)合仿真方法可應(yīng)用于機(jī)械式PTO 的振蕩浮子式波能裝備設(shè)計(jì)，為不同海況與各類振蕩浮子式波浪能發(fā)電裝置的機(jī)械參數(shù)優(yōu)化與選型工作提供參考。

1 浮子水動(dòng)力模型的建立

1.1 控制方程

將水視為粘性不可壓縮流體，其控制方程為連續(xù)性方程和納維-斯托克斯（N-S）方程，可表示為

式中，vx、vy和vz分別表示x、y和z方向的速度分量，fx、fy和fz分別表示x、y和z方向受到的質(zhì)量力，?2為拉普拉斯算子，t表示時(shí)間，ρ表示流體密度，p表示流體壓力，ν表示流體的運(yùn)動(dòng)粘度。

由于流域包括氣、液兩相，因此本文以流體體積分?jǐn)?shù)法（VOF）對(duì)自由液面進(jìn)行追蹤和捕捉。其中，相的分布和液面的位置由相體積分?jǐn)?shù)描述，相體積分?jǐn)?shù)αi定義為

式中，Vi為網(wǎng)格單元中第i相的體積，V為網(wǎng)格單元總體積。網(wǎng)格單元中所有相的體積分?jǐn)?shù)滿足下式：

式中，N為網(wǎng)格單元中相的總個(gè)數(shù)，文中N=2。

1.2 三維數(shù)值波浪水池的建立與驗(yàn)證

采用的數(shù)值水池如圖1 所示，流域左側(cè)為速度入口邊界，右側(cè)為壓力出口邊界，兩側(cè)和底部均為無(wú)滑移壁面邊界，頂部為速度入口邊界。其中上游和下游均設(shè)置3倍波長(zhǎng)，兩側(cè)設(shè)置2.5倍波長(zhǎng)，四周均設(shè)置消波區(qū)，沿計(jì)算域邊界設(shè)置了力消波區(qū)域[14]。

圖1 數(shù)值水池模型示意圖Fig.1 Schematic diagram of numerical pool model

采用切割體網(wǎng)格進(jìn)行網(wǎng)格劃分。由于涉及到浮子的運(yùn)動(dòng)，需采用重疊網(wǎng)格。為了更好地捕捉自由液面，對(duì)兩相交界處進(jìn)行加密處理。網(wǎng)格的劃分如圖2所示。

圖2 搖臂浮子網(wǎng)格劃分示意圖Fig.2 Diagram of float-arm buoy mesh division

為驗(yàn)證數(shù)值水池的可行性，對(duì)上述建立的三維波浪水池進(jìn)行網(wǎng)格收斂性驗(yàn)證和時(shí)間步收斂性驗(yàn)證。首先選取了三種網(wǎng)格尺寸，具體網(wǎng)格數(shù)量與尺寸如表1所示。

表1 網(wǎng)格尺寸與數(shù)量Tab.1 Mesh size and quantity

由圖3 可以看出，與網(wǎng)格A 的垂蕩振幅對(duì)比，網(wǎng)格B 和網(wǎng)格C的相對(duì)誤差分別為1.43%和6.14%，h代表浮子垂蕩。綜合考慮數(shù)值精度和計(jì)算成本，最終采用網(wǎng)格B。

圖3 網(wǎng)格收斂性驗(yàn)證Fig.3 Grid convergence verification

接下來(lái)根據(jù)網(wǎng)格B 尺寸的數(shù)值模型進(jìn)行時(shí)間步驗(yàn)證。由圖4可以看出，時(shí)間步Δt=0.02 s時(shí)，浮子垂蕩曲線明顯與Δt=0.01 s 和Δt=0.005 s 有偏移，而Δt=0.01 s 和Δt=0.005 s 的垂蕩振幅相對(duì)誤差為1.28%。故在后續(xù)計(jì)算中采用Δt=0.01 s的時(shí)間步長(zhǎng)。

圖4 時(shí)間步收斂性驗(yàn)證Fig.4 Time step convergence verification

2 PTO系統(tǒng)仿真模型

2.1 波浪能發(fā)電裝置運(yùn)行模式簡(jiǎn)介

本文設(shè)計(jì)的波浪能發(fā)電裝置可分為三部分，分別是搖臂浮子、機(jī)械傳動(dòng)裝置和發(fā)電機(jī)，如圖5所示。

圖5 波能發(fā)電裝置組成示意圖（1-搖臂浮子，2-機(jī)械傳動(dòng)裝置，3-發(fā)電機(jī)；俯視圖）Fig.5 Schematic diagram of WEC compo?sition(1-Buoy,2-Mechanical trans?mission,3-Motor;top view)

搖臂浮子作為發(fā)電裝置的波能捕獲單元，在波浪激勵(lì)力的作用下繞鉸接軸做縱搖運(yùn)動(dòng)，浮子的縱搖通過(guò)機(jī)械傳動(dòng)裝置將運(yùn)動(dòng)傳遞給發(fā)電機(jī)。其具體的工作流程為：浮子的縱搖通過(guò)同步帶將運(yùn)動(dòng)傳遞給兩根帶輪軸，兩根帶輪軸上分別安裝有工作方向相反的超越離合器，單個(gè)入射波周期內(nèi)，兩根帶輪軸將通過(guò)齒輪組交替向發(fā)電機(jī)的齒輪軸傳遞動(dòng)力，經(jīng)齒輪箱增速后提高發(fā)電機(jī)的轉(zhuǎn)速。由于浮子的上升與下降行程分別對(duì)應(yīng)兩個(gè)超越離合器的單程做功，當(dāng)研究單個(gè)行程的做功特性時(shí)，也可以在本機(jī)械裝置上實(shí)現(xiàn)。

2.2 PTO系統(tǒng)的仿真模型建立

AMESim 是系統(tǒng)工程高級(jí)建模和仿真平臺(tái)，用以建立復(fù)雜的多學(xué)科多領(lǐng)域的系統(tǒng)模型，包括流體、機(jī)械、電磁和控制等，并在此基礎(chǔ)上進(jìn)行仿真計(jì)算和分析，研究元件或系統(tǒng)的穩(wěn)態(tài)與動(dòng)態(tài)性能。本文集成機(jī)械、電氣與信號(hào)元件搭建出PTO 系統(tǒng)的完整模型，并參照仿真結(jié)果對(duì)系統(tǒng)進(jìn)行優(yōu)化。根據(jù)上述機(jī)械傳動(dòng)裝置的工作原理，通過(guò)AMESim 搭建了波能發(fā)電裝置PTO 系統(tǒng)的仿真模型，如圖6所示。

圖6 波能發(fā)電裝置PTO系統(tǒng)仿真模型Fig.6 PTO system simulation model of WEC

上述AMESim 仿真模型中，機(jī)械部分包括超越離合器、齒輪系、帶輪系與齒輪箱。當(dāng)超越離合器處于接合狀態(tài)時(shí)，將超越離合器視作一個(gè)無(wú)轉(zhuǎn)動(dòng)慣量的剛體，根據(jù)工作方向傳遞轉(zhuǎn)矩；當(dāng)超越離合器處于超越狀態(tài)時(shí)，動(dòng)力傳遞中斷。以其主動(dòng)端與從動(dòng)端的轉(zhuǎn)角差作為模型中工作狀態(tài)的判斷依據(jù)，則超越離合器傳遞轉(zhuǎn)矩M(t)的表達(dá)式如下：

式中，Kc為超越離合器扭轉(zhuǎn)剛度，μ為等效粘性摩擦系數(shù)，θm(t)為主動(dòng)端扭轉(zhuǎn)角度，θs(t)為從動(dòng)端扭轉(zhuǎn)角度。

忽略超越離合器主從動(dòng)件之間的間隙與變形因素的影響，則機(jī)械傳動(dòng)裝置的運(yùn)動(dòng)方程可表示為

式中，ωb為浮子的角速度，ωM為電機(jī)轉(zhuǎn)子的角速度，γ為等效傳動(dòng)比，ig為增速齒輪箱傳動(dòng)比，rb為同步帶輪的外徑比，rg為齒輪組的分度圓半徑比。

改變AMESim 內(nèi)信號(hào)部分的觸發(fā)方式，分別建立了下降行程做功與上升行程做功兩種單行程的動(dòng)力傳遞方式，波能發(fā)電裝置的仿真模型如圖7～8所示。

圖7 浮子下降行程做功的PTO系統(tǒng)仿真模型Fig.7 PTO system simulation model of downward working buoy

圖8 浮子上升行程做功的PTO系統(tǒng)仿真模型Fig.8 PTO system simulation model of upward working buoy

AMESim 電氣部分以直流電機(jī)作為PTO 系統(tǒng)的發(fā)電機(jī)，電阻作為PTO 系統(tǒng)的電子負(fù)載，用電阻功率來(lái)衡量PTO 系統(tǒng)的發(fā)電功率。上述三種仿真模型的電氣部分與機(jī)械部分的參數(shù)設(shè)置完全相同，只根據(jù)做功方式改變了信號(hào)部分的觸發(fā)方式，用以執(zhí)行聯(lián)合仿真。

3 水動(dòng)力-PTO聯(lián)合仿真設(shè)置

本文開展了水動(dòng)力模型與PTO 系統(tǒng)聯(lián)合仿真研究，計(jì)算出了在PTO 系統(tǒng)作用下，浮子角速度、電機(jī)功率與電機(jī)轉(zhuǎn)速等相關(guān)數(shù)據(jù)。其仿真過(guò)程如圖9所示。

圖9 STAR-CCM+與AMESim聯(lián)合仿真示意圖Fig.9 Co-simulation diagram of STAR-CCM+and AMESim

AMESim 將PTO 的負(fù)載轉(zhuǎn)矩傳遞到STAR-CCM+中的水動(dòng)力計(jì)算模塊，同時(shí)考慮了電機(jī)轉(zhuǎn)子轉(zhuǎn)動(dòng)慣量的慣性轉(zhuǎn)矩與電機(jī)阻尼的電磁轉(zhuǎn)矩，兩者共同組成了電機(jī)的負(fù)載轉(zhuǎn)矩。電機(jī)的負(fù)載轉(zhuǎn)矩通過(guò)超越離合器和增速齒輪箱等機(jī)械傳動(dòng)系統(tǒng)轉(zhuǎn)變?yōu)樽饔迷赟TAR-CCM+模型上浮子的負(fù)載力矩。因此，本聯(lián)合仿真方法相較于直接添加阻尼系數(shù)的計(jì)算方法更加接近真實(shí)PTO系統(tǒng)的工作狀態(tài)。

本文中的計(jì)算模型為帶擺臂的圓柱形浮子，所用海況來(lái)源于威海國(guó)家淺海綜合試驗(yàn)場(chǎng)，Star-CCM+內(nèi)的浮子參數(shù)、海況參數(shù)設(shè)置如表2所示。

表2 Star-CCM+計(jì)算模型參數(shù)Tab.2 Calculation parameters of the model in Star-CCM+

浮子后端的PTO 系統(tǒng)仿真模型通過(guò)AMESim 搭建，發(fā)電機(jī)選用Z4-100-1 型號(hào)直流電機(jī)，AMESim內(nèi)的電機(jī)參數(shù)、機(jī)械元件參數(shù)設(shè)置如表3所示。

表3 AMESim計(jì)算模型參數(shù)Tab.3 Calculation parameters of the model in AMESim

下文將往返雙行程做功的浮子稱為1 號(hào)浮子，上升行程做功的浮子稱為2 號(hào)浮子，下降行程做功的浮子稱為3號(hào)浮子。圖10為雙行程做功浮子在一個(gè)周期內(nèi)的運(yùn)動(dòng)響應(yīng)。為研究做功模式以及等效傳動(dòng)比對(duì)浮子的角速度與電機(jī)功率的影響，對(duì)工況進(jìn)行設(shè)置，如表4所示，并根據(jù)仿真結(jié)果，在特定區(qū)間完成PTO系統(tǒng)的負(fù)載匹配與優(yōu)化選型。

表4 工況設(shè)置Tab.4 Setting of working conditions

圖10 單個(gè)波周期內(nèi)浮子的運(yùn)動(dòng)響應(yīng)Fig.10 Motion response of the buoy in a single wave period

4 結(jié)果分析

1 號(hào)、2 號(hào)和3 號(hào)浮子在不同工況下的浮子角速度和電機(jī)瞬時(shí)發(fā)電功率PI如圖11～15 所示。由于采用直流電機(jī)，故電機(jī)瞬時(shí)功率與電機(jī)轉(zhuǎn)子的角速度成正比。圖中，當(dāng)角速度小于0 時(shí)，浮子逆時(shí)針轉(zhuǎn)動(dòng)，對(duì)應(yīng)上升行程做功；角速度大于0時(shí)，浮子順時(shí)針轉(zhuǎn)動(dòng)，對(duì)應(yīng)下降行程做功。

圖11 1:20傳動(dòng)比的浮子角速度和電機(jī)瞬時(shí)功率Fig.11 Buoy angular velocity and instantaneous power under 1:20 transmission ratio

圖12 1:40傳動(dòng)比的浮子角速度和電機(jī)瞬時(shí)功率Fig.12 Buoy angular velocity and instantaneous power under 1:40 transmission ratio

可以看出，在低傳動(dòng)比的情況下（圖11～12），三種浮子的電機(jī)功率與浮子角速度曲線吻合較好，說(shuō)明低傳動(dòng)比情況下PTO 負(fù)載力矩小，動(dòng)力傳遞狀態(tài)下PTO 系統(tǒng)對(duì)浮子運(yùn)動(dòng)的影響幾乎可以忽略。隨著傳動(dòng)比的增大（圖13～15），2號(hào)浮子對(duì)應(yīng)的1號(hào)浮子上升行程段與3號(hào)浮子對(duì)應(yīng)的1號(hào)浮子下降行程段，角速度曲線前半部分吻合較差，后半部分吻合良好。其原因?yàn)椋合噍^于雙行程的做功，浮子在單行程做功狀態(tài)下只有單行程內(nèi)存在負(fù)載力矩，故從自由運(yùn)動(dòng)切換至負(fù)載狀態(tài)時(shí)的負(fù)載力矩發(fā)生突變，導(dǎo)致前半段吻合較差，待后半段運(yùn)動(dòng)穩(wěn)定后，功率曲線與角速度曲線吻合良好。在高傳動(dòng)比情況下，隨著PTO 系統(tǒng)的負(fù)載力矩增大，會(huì)顯著影響浮子的運(yùn)動(dòng)狀態(tài)，工作行程內(nèi)浮子角速度的幅值將顯著降低。

圖13 1:60傳動(dòng)比的浮子角速度和電機(jī)瞬時(shí)功率Fig.13 Buoy angular velocity and instantaneous power under 1:60 transmission ratio

圖14 1:80電機(jī)瞬時(shí)功率傳動(dòng)比的浮子角速度和電機(jī)瞬時(shí)功率Fig.14 Buoy angular velocity and instantaneous power under 1:80 transmission ratio

圖15 1:100傳動(dòng)比的浮子角速度電機(jī)瞬時(shí)功率Fig.15 Buoy angular velocity and instantaneous power under 1:100 transmission ratio

圖16 為1 號(hào)、2 號(hào)和3 號(hào)浮子在不同傳動(dòng)比工況下一個(gè)入射波周期內(nèi)電機(jī)的平均功率圖。平均功率PT的計(jì)算公式如下：

圖16 單個(gè)波周期內(nèi)三種浮子的平均功率Fig.16 Average power of three buoys in a single wave period

從圖中可以看出，1 號(hào)浮子的輸出功率始終大于2 號(hào)與3 號(hào)浮子，這是因?yàn)? 號(hào)浮子在一個(gè)波周期內(nèi)對(duì)PTO 系統(tǒng)全程做功。在傳動(dòng)比小于40 時(shí)，1 號(hào)浮子的平均功率是2號(hào)和3號(hào)浮子的近兩倍。但隨著傳動(dòng)比的增大，1號(hào)浮子的平均功率衰減嚴(yán)重，2號(hào)與3號(hào)浮子的功率逐漸超過(guò)1號(hào)浮子的單程功率，在傳動(dòng)比為80到100之間時(shí)，1號(hào)浮子的輸出功率僅為3號(hào)浮子的1.3倍。這是因?yàn)樵诟邆鲃?dòng)比情況下，浮子受到的負(fù)載力矩過(guò)大，1號(hào)浮子在雙行程均受到較大的負(fù)載，導(dǎo)致其運(yùn)動(dòng)受阻，無(wú)法完全吸收波浪的能量。而2號(hào)和3號(hào)浮子，只有一段行程受到了較大的負(fù)載力，故單段輸出功率較1號(hào)浮子更高。

由于計(jì)算模型中考慮了擺臂的重力影響，且3號(hào)浮子在上升過(guò)程不受PTO 的負(fù)載力矩，在一個(gè)波周期內(nèi)吸收了更多波浪的能量，使得其輸出功率會(huì)略大于2號(hào)浮子。

此外，隨著傳動(dòng)比的增大，三種浮子的平均功率均先上升后下降，說(shuō)明PTO 系統(tǒng)存在最佳傳動(dòng)比，在60 到80 優(yōu)化區(qū)間內(nèi)等差選取10 組增速比工況進(jìn)行計(jì)算，其平均功率如表5 所示。據(jù)表可知，電機(jī)的最大平均功率出現(xiàn)在70 至75 增速比之間，故在波能裝置設(shè)計(jì)中其最佳傳動(dòng)比應(yīng)取70 至75 之間，并在此區(qū)間內(nèi)配合選擇齒輪箱的參數(shù)與齒輪帶輪的外徑。

表5 優(yōu)化區(qū)間1號(hào)浮子平均功率Tab.5 No.1 buoy average power in optimized interval

5 結(jié) 論

（1）本文基于發(fā)電機(jī)和機(jī)械裝置的工作原理，利用AMESim 搭建了PTO 負(fù)載模型，相比設(shè)置線性阻尼系數(shù)更加貼近真實(shí)情況。

（2）低傳動(dòng)比狀態(tài)下，負(fù)載力矩對(duì)浮子的運(yùn)動(dòng)影響很小。隨著傳動(dòng)比的增大，浮子的角速度將逐漸衰減。聯(lián)合仿真模型可通過(guò)改變傳動(dòng)比來(lái)尋找PTO 系統(tǒng)輸出功率的最佳參數(shù)區(qū)間，進(jìn)而確定元件的最佳選型，在本文海況中，1 m 的圓柱形浮子機(jī)械裝置的最佳傳動(dòng)比在72.5 附近，最大平均功率為35.1 W。

（3）三種做功方式中，浮子雙行程做功的平均功率要大于單行程的平均功率。隨著傳動(dòng)比的增大，其平均功率先增大再減小，傳動(dòng)比增大到80之后，PTO系統(tǒng)對(duì)浮子的負(fù)載力矩過(guò)大，浮子的縱搖角速度衰減嚴(yán)重，浮子無(wú)法充分吸收波浪的能量，導(dǎo)致雙行程的平均功率降低為單行程的1.3倍。

（4）聯(lián)合仿真的計(jì)算結(jié)果能指導(dǎo)裝置最大平均發(fā)電功率的參數(shù)優(yōu)化工作，本文的聯(lián)合仿真模型可用于不同海況與不同浮子尺寸的振蕩浮子式波浪能發(fā)電裝置的設(shè)計(jì)、機(jī)械傳動(dòng)系統(tǒng)的參數(shù)優(yōu)化和發(fā)電機(jī)的選型。

（5）聯(lián)合仿真模型涉及水動(dòng)力學(xué)與機(jī)械等學(xué)科領(lǐng)域，可用于指導(dǎo)實(shí)際工程中振蕩浮子式波浪能發(fā)電裝置的功率優(yōu)化工作，確定最優(yōu)功率下的發(fā)電機(jī)型號(hào)、齒輪箱增速比、浮子尺寸等數(shù)據(jù)，降低實(shí)際開發(fā)成本。在未來(lái)多浮子波浪能發(fā)電裝置的開發(fā)中，聯(lián)合仿真模型可以計(jì)算更復(fù)雜的PTO負(fù)載，從而簡(jiǎn)化多浮子波能裝置中PTO 的機(jī)械設(shè)計(jì)與浮子的數(shù)量和布置等復(fù)雜的優(yōu)化問(wèn)題，縮短實(shí)際工程中波能裝置的開發(fā)周期，為不同海況下的振蕩浮子式波浪能發(fā)電裝置的設(shè)計(jì)提供借鑒。