何廣華,楊凱博,欒政曉,張志剛,劉朝綱,荊芃霖
(1.哈爾濱工業(yè)大學(威海)海洋工程學院,山東威海264209;2.山東船舶技術(shù)研究院,山東威海264209;3.哈爾濱工業(yè)大學機械工程學院,哈爾濱150001)
由于能源消耗和環(huán)境污染等問題,海洋可再生能源的開發(fā)和利用逐漸受到各國的關注[1]。其中波浪能以儲量大、無污染、可重復開發(fā)利用等優(yōu)點,成為國內(nèi)外海洋能開發(fā)利用研究的熱點。波浪能發(fā)電裝置根據(jù)其工作原理可分為振蕩體式、振蕩水柱式、聚波越浪式等[2]。其中,振蕩體式又可分為振蕩浮子式(點吸收式)、筏式和擺式。近年來,振蕩浮子式WEC 因其體積小、制造安裝方便、不受波浪方向影響等優(yōu)勢受到了廣泛的關注[3]。常見的振蕩浮子式WEC 的PTO 系統(tǒng)可分為機械式和液壓式,其中液壓式PTO 利用液壓缸和液壓泵等元件將液壓能轉(zhuǎn)化成電能,機械式PTO 利用齒輪齒條、滾珠絲杠[4-5]等機械裝置將機械能轉(zhuǎn)化成電能。
波浪能發(fā)電裝置的PTO 系統(tǒng)是波浪能發(fā)電的關鍵技術(shù)之一,不同研究者對PTO 系統(tǒng)進行了不同程度的簡化研究。Rico 等[6-8]利用數(shù)學公式計算了搖臂圓柱形浮子的運動學方程,分析了線性PTO 阻尼對浮子俘獲功率的影響;張亞群等[9]利用物模試驗分析了振蕩浮子的水動力性能,實驗顯示在最優(yōu)PTO 負載下浮子可獲得最大俘獲寬度比;宋文杰等[10]利用實驗分析了雙行程工作形式PTO 的輸出特性,但未對單行程工作形式進行對比分析;馮輝等[11]基于Simulink和AMESim 搭建了聯(lián)合仿真模型,該模型包括了各液壓元件與發(fā)電機,較好地模擬出了能量轉(zhuǎn)換系統(tǒng)的發(fā)電特性曲線,但其未根據(jù)仿真結(jié)果對PTO 元件進行優(yōu)化選型;Li等[12]利用CFD 工具研究了波物相互作用與機械元件對浮子的影響,發(fā)現(xiàn)機械元件產(chǎn)生的力和力矩對單個浮子的運動響應影響巨大;劉延偉等[13]對超越離合器進行了理想建模,并分析了在切換瞬間超越離合器的動態(tài)特性。
然而,對于采用機械式PTO 的振蕩浮子波浪能發(fā)電裝置,PTO 系統(tǒng)的傳動方式、發(fā)電機選型及電機的電子負載等參數(shù)都會對波浪能發(fā)電裝置的獲能產(chǎn)生影響,導致將PTO 負載簡化為線性阻尼的方法不夠精確。針對該問題,本文采用水動力模型與PTO系統(tǒng)聯(lián)合仿真的方法,使所添加的負載更接近真實情況,其中,在STAR-CCM+中建立波浪能發(fā)電裝置數(shù)值模型,在AMESim 中建立PTO 系統(tǒng)仿真模型,通過二者的聯(lián)合仿真進行水動力與PTO系統(tǒng)的耦合研究;研究PTO負載的不同行程工作模式與不同傳動比對浮子運動響應和輸出功率的影響,并針對目標海況對PTO 參數(shù)進行優(yōu)化。本文使用的聯(lián)合仿真方法可應用于機械式PTO 的振蕩浮子式波能裝備設計,為不同海況與各類振蕩浮子式波浪能發(fā)電裝置的機械參數(shù)優(yōu)化與選型工作提供參考。
將水視為粘性不可壓縮流體,其控制方程為連續(xù)性方程和納維-斯托克斯(N-S)方程,可表示為
式中,vx、vy和vz分別表示x、y和z方向的速度分量,fx、fy和fz分別表示x、y和z方向受到的質(zhì)量力,?2為拉普拉斯算子,t表示時間,ρ表示流體密度,p表示流體壓力,ν表示流體的運動粘度。
由于流域包括氣、液兩相,因此本文以流體體積分數(shù)法(VOF)對自由液面進行追蹤和捕捉。其中,相的分布和液面的位置由相體積分數(shù)描述,相體積分數(shù)αi定義為
式中,Vi為網(wǎng)格單元中第i相的體積,V為網(wǎng)格單元總體積。網(wǎng)格單元中所有相的體積分數(shù)滿足下式:
式中,N為網(wǎng)格單元中相的總個數(shù),文中N=2。
采用的數(shù)值水池如圖1 所示,流域左側(cè)為速度入口邊界,右側(cè)為壓力出口邊界,兩側(cè)和底部均為無滑移壁面邊界,頂部為速度入口邊界。其中上游和下游均設置3倍波長,兩側(cè)設置2.5倍波長,四周均設置消波區(qū),沿計算域邊界設置了力消波區(qū)域[14]。

圖1 數(shù)值水池模型示意圖Fig.1 Schematic diagram of numerical pool model
采用切割體網(wǎng)格進行網(wǎng)格劃分。由于涉及到浮子的運動,需采用重疊網(wǎng)格。為了更好地捕捉自由液面,對兩相交界處進行加密處理。網(wǎng)格的劃分如圖2所示。

圖2 搖臂浮子網(wǎng)格劃分示意圖Fig.2 Diagram of float-arm buoy mesh division
為驗證數(shù)值水池的可行性,對上述建立的三維波浪水池進行網(wǎng)格收斂性驗證和時間步收斂性驗證。首先選取了三種網(wǎng)格尺寸,具體網(wǎng)格數(shù)量與尺寸如表1所示。

表1 網(wǎng)格尺寸與數(shù)量Tab.1 Mesh size and quantity
由圖3 可以看出,與網(wǎng)格A 的垂蕩振幅對比,網(wǎng)格B 和網(wǎng)格C的相對誤差分別為1.43%和6.14%,h代表浮子垂蕩。綜合考慮數(shù)值精度和計算成本,最終采用網(wǎng)格B。

圖3 網(wǎng)格收斂性驗證Fig.3 Grid convergence verification
接下來根據(jù)網(wǎng)格B 尺寸的數(shù)值模型進行時間步驗證。由圖4可以看出,時間步Δt=0.02 s時,浮子垂蕩曲線明顯與Δt=0.01 s 和Δt=0.005 s 有偏移,而Δt=0.01 s 和Δt=0.005 s 的垂蕩振幅相對誤差為1.28%。故在后續(xù)計算中采用Δt=0.01 s的時間步長。

圖4 時間步收斂性驗證Fig.4 Time step convergence verification
本文設計的波浪能發(fā)電裝置可分為三部分,分別是搖臂浮子、機械傳動裝置和發(fā)電機,如圖5所示。

圖5 波能發(fā)電裝置組成示意圖(1-搖臂浮子,2-機械傳動裝置,3-發(fā)電機;俯視圖)Fig.5 Schematic diagram of WEC compo?sition(1-Buoy,2-Mechanical trans?mission,3-Motor;top view)
搖臂浮子作為發(fā)電裝置的波能捕獲單元,在波浪激勵力的作用下繞鉸接軸做縱搖運動,浮子的縱搖通過機械傳動裝置將運動傳遞給發(fā)電機。其具體的工作流程為:浮子的縱搖通過同步帶將運動傳遞給兩根帶輪軸,兩根帶輪軸上分別安裝有工作方向相反的超越離合器,單個入射波周期內(nèi),兩根帶輪軸將通過齒輪組交替向發(fā)電機的齒輪軸傳遞動力,經(jīng)齒輪箱增速后提高發(fā)電機的轉(zhuǎn)速。由于浮子的上升與下降行程分別對應兩個超越離合器的單程做功,當研究單個行程的做功特性時,也可以在本機械裝置上實現(xiàn)。
AMESim 是系統(tǒng)工程高級建模和仿真平臺,用以建立復雜的多學科多領域的系統(tǒng)模型,包括流體、機械、電磁和控制等,并在此基礎上進行仿真計算和分析,研究元件或系統(tǒng)的穩(wěn)態(tài)與動態(tài)性能。本文集成機械、電氣與信號元件搭建出PTO 系統(tǒng)的完整模型,并參照仿真結(jié)果對系統(tǒng)進行優(yōu)化。根據(jù)上述機械傳動裝置的工作原理,通過AMESim 搭建了波能發(fā)電裝置PTO 系統(tǒng)的仿真模型,如圖6所示。

圖6 波能發(fā)電裝置PTO系統(tǒng)仿真模型Fig.6 PTO system simulation model of WEC
上述AMESim 仿真模型中,機械部分包括超越離合器、齒輪系、帶輪系與齒輪箱。當超越離合器處于接合狀態(tài)時,將超越離合器視作一個無轉(zhuǎn)動慣量的剛體,根據(jù)工作方向傳遞轉(zhuǎn)矩;當超越離合器處于超越狀態(tài)時,動力傳遞中斷。以其主動端與從動端的轉(zhuǎn)角差作為模型中工作狀態(tài)的判斷依據(jù),則超越離合器傳遞轉(zhuǎn)矩M(t)的表達式如下:
式中,Kc為超越離合器扭轉(zhuǎn)剛度,μ為等效粘性摩擦系數(shù),θm(t)為主動端扭轉(zhuǎn)角度,θs(t)為從動端扭轉(zhuǎn)角度。
忽略超越離合器主從動件之間的間隙與變形因素的影響,則機械傳動裝置的運動方程可表示為
式中,ωb為浮子的角速度,ωM為電機轉(zhuǎn)子的角速度,γ為等效傳動比,ig為增速齒輪箱傳動比,rb為同步帶輪的外徑比,rg為齒輪組的分度圓半徑比。
改變AMESim 內(nèi)信號部分的觸發(fā)方式,分別建立了下降行程做功與上升行程做功兩種單行程的動力傳遞方式,波能發(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 電氣部分以直流電機作為PTO 系統(tǒng)的發(fā)電機,電阻作為PTO 系統(tǒng)的電子負載,用電阻功率來衡量PTO 系統(tǒng)的發(fā)電功率。上述三種仿真模型的電氣部分與機械部分的參數(shù)設置完全相同,只根據(jù)做功方式改變了信號部分的觸發(fā)方式,用以執(zhí)行聯(lián)合仿真。
本文開展了水動力模型與PTO 系統(tǒng)聯(lián)合仿真研究,計算出了在PTO 系統(tǒng)作用下,浮子角速度、電機功率與電機轉(zhuǎn)速等相關數(shù)據(jù)。其仿真過程如圖9所示。

圖9 STAR-CCM+與AMESim聯(lián)合仿真示意圖Fig.9 Co-simulation diagram of STAR-CCM+and AMESim
AMESim 將PTO 的負載轉(zhuǎn)矩傳遞到STAR-CCM+中的水動力計算模塊,同時考慮了電機轉(zhuǎn)子轉(zhuǎn)動慣量的慣性轉(zhuǎn)矩與電機阻尼的電磁轉(zhuǎn)矩,兩者共同組成了電機的負載轉(zhuǎn)矩。電機的負載轉(zhuǎn)矩通過超越離合器和增速齒輪箱等機械傳動系統(tǒng)轉(zhuǎn)變?yōu)樽饔迷赟TAR-CCM+模型上浮子的負載力矩。因此,本聯(lián)合仿真方法相較于直接添加阻尼系數(shù)的計算方法更加接近真實PTO系統(tǒng)的工作狀態(tài)。
本文中的計算模型為帶擺臂的圓柱形浮子,所用海況來源于威海國家淺海綜合試驗場,Star-CCM+內(nèi)的浮子參數(shù)、海況參數(shù)設置如表2所示。

表2 Star-CCM+計算模型參數(shù)Tab.2 Calculation parameters of the model in Star-CCM+
浮子后端的PTO 系統(tǒng)仿真模型通過AMESim 搭建,發(fā)電機選用Z4-100-1 型號直流電機,AMESim內(nèi)的電機參數(shù)、機械元件參數(shù)設置如表3所示。

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

表4 工況設置Tab.4 Setting of working conditions

圖10 單個波周期內(nèi)浮子的運動響應Fig.10 Motion response of the buoy in a single wave period
1 號、2 號和3 號浮子在不同工況下的浮子角速度和電機瞬時發(fā)電功率PI如圖11~15 所示。由于采用直流電機,故電機瞬時功率與電機轉(zhuǎn)子的角速度成正比。圖中,當角速度小于0 時,浮子逆時針轉(zhuǎn)動,對應上升行程做功;角速度大于0時,浮子順時針轉(zhuǎn)動,對應下降行程做功。

圖11 1:20傳動比的浮子角速度和電機瞬時功率Fig.11 Buoy angular velocity and instantaneous power under 1:20 transmission ratio

圖12 1:40傳動比的浮子角速度和電機瞬時功率Fig.12 Buoy angular velocity and instantaneous power under 1:40 transmission ratio
可以看出,在低傳動比的情況下(圖11~12),三種浮子的電機功率與浮子角速度曲線吻合較好,說明低傳動比情況下PTO 負載力矩小,動力傳遞狀態(tài)下PTO 系統(tǒng)對浮子運動的影響幾乎可以忽略。隨著傳動比的增大(圖13~15),2號浮子對應的1號浮子上升行程段與3號浮子對應的1號浮子下降行程段,角速度曲線前半部分吻合較差,后半部分吻合良好。其原因為:相較于雙行程的做功,浮子在單行程做功狀態(tài)下只有單行程內(nèi)存在負載力矩,故從自由運動切換至負載狀態(tài)時的負載力矩發(fā)生突變,導致前半段吻合較差,待后半段運動穩(wěn)定后,功率曲線與角速度曲線吻合良好。在高傳動比情況下,隨著PTO 系統(tǒng)的負載力矩增大,會顯著影響浮子的運動狀態(tài),工作行程內(nèi)浮子角速度的幅值將顯著降低。

圖13 1:60傳動比的浮子角速度和電機瞬時功率Fig.13 Buoy angular velocity and instantaneous power under 1:60 transmission ratio

圖14 1:80電機瞬時功率傳動比的浮子角速度和電機瞬時功率Fig.14 Buoy angular velocity and instantaneous power under 1:80 transmission ratio

圖15 1:100傳動比的浮子角速度電機瞬時功率Fig.15 Buoy angular velocity and instantaneous power under 1:100 transmission ratio
圖16 為1 號、2 號和3 號浮子在不同傳動比工況下一個入射波周期內(nèi)電機的平均功率圖。平均功率PT的計算公式如下:

圖16 單個波周期內(nèi)三種浮子的平均功率Fig.16 Average power of three buoys in a single wave period
從圖中可以看出,1 號浮子的輸出功率始終大于2 號與3 號浮子,這是因為1 號浮子在一個波周期內(nèi)對PTO 系統(tǒng)全程做功。在傳動比小于40 時,1 號浮子的平均功率是2號和3號浮子的近兩倍。但隨著傳動比的增大,1號浮子的平均功率衰減嚴重,2號與3號浮子的功率逐漸超過1號浮子的單程功率,在傳動比為80到100之間時,1號浮子的輸出功率僅為3號浮子的1.3倍。這是因為在高傳動比情況下,浮子受到的負載力矩過大,1號浮子在雙行程均受到較大的負載,導致其運動受阻,無法完全吸收波浪的能量。而2號和3號浮子,只有一段行程受到了較大的負載力,故單段輸出功率較1號浮子更高。
由于計算模型中考慮了擺臂的重力影響,且3號浮子在上升過程不受PTO 的負載力矩,在一個波周期內(nèi)吸收了更多波浪的能量,使得其輸出功率會略大于2號浮子。
此外,隨著傳動比的增大,三種浮子的平均功率均先上升后下降,說明PTO 系統(tǒng)存在最佳傳動比,在60 到80 優(yōu)化區(qū)間內(nèi)等差選取10 組增速比工況進行計算,其平均功率如表5 所示。據(jù)表可知,電機的最大平均功率出現(xiàn)在70 至75 增速比之間,故在波能裝置設計中其最佳傳動比應取70 至75 之間,并在此區(qū)間內(nèi)配合選擇齒輪箱的參數(shù)與齒輪帶輪的外徑。

表5 優(yōu)化區(qū)間1號浮子平均功率Tab.5 No.1 buoy average power in optimized interval
(1)本文基于發(fā)電機和機械裝置的工作原理,利用AMESim 搭建了PTO 負載模型,相比設置線性阻尼系數(shù)更加貼近真實情況。
(2)低傳動比狀態(tài)下,負載力矩對浮子的運動影響很小。隨著傳動比的增大,浮子的角速度將逐漸衰減。聯(lián)合仿真模型可通過改變傳動比來尋找PTO 系統(tǒng)輸出功率的最佳參數(shù)區(qū)間,進而確定元件的最佳選型,在本文海況中,1 m 的圓柱形浮子機械裝置的最佳傳動比在72.5 附近,最大平均功率為35.1 W。
(3)三種做功方式中,浮子雙行程做功的平均功率要大于單行程的平均功率。隨著傳動比的增大,其平均功率先增大再減小,傳動比增大到80之后,PTO系統(tǒng)對浮子的負載力矩過大,浮子的縱搖角速度衰減嚴重,浮子無法充分吸收波浪的能量,導致雙行程的平均功率降低為單行程的1.3倍。
(4)聯(lián)合仿真的計算結(jié)果能指導裝置最大平均發(fā)電功率的參數(shù)優(yōu)化工作,本文的聯(lián)合仿真模型可用于不同海況與不同浮子尺寸的振蕩浮子式波浪能發(fā)電裝置的設計、機械傳動系統(tǒng)的參數(shù)優(yōu)化和發(fā)電機的選型。
(5)聯(lián)合仿真模型涉及水動力學與機械等學科領域,可用于指導實際工程中振蕩浮子式波浪能發(fā)電裝置的功率優(yōu)化工作,確定最優(yōu)功率下的發(fā)電機型號、齒輪箱增速比、浮子尺寸等數(shù)據(jù),降低實際開發(fā)成本。在未來多浮子波浪能發(fā)電裝置的開發(fā)中,聯(lián)合仿真模型可以計算更復雜的PTO負載,從而簡化多浮子波能裝置中PTO 的機械設計與浮子的數(shù)量和布置等復雜的優(yōu)化問題,縮短實際工程中波能裝置的開發(fā)周期,為不同海況下的振蕩浮子式波浪能發(fā)電裝置的設計提供借鑒。