何廣華，汪 鵬，劉朝綱，欒政曉

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

振蕩浮子式波浪能發(fā)電裝置根據(jù)工作原理的不同可分為機(jī)械式[1]和液壓式[2-3]。機(jī)械式是通過齒輪齒條等機(jī)械傳動(dòng)裝置將浮子的運(yùn)動(dòng)轉(zhuǎn)化為旋轉(zhuǎn)運(yùn)動(dòng)，然后通過齒輪箱進(jìn)行增速，最后在發(fā)電環(huán)節(jié)中通過發(fā)電機(jī)將旋轉(zhuǎn)運(yùn)動(dòng)的機(jī)械能轉(zhuǎn)化為電能。液壓式則是利用液壓系統(tǒng)進(jìn)行發(fā)電。對于機(jī)械式的振蕩浮子式波浪能發(fā)電裝置，在其設(shè)計(jì)過程中，需要對某一確定的浮子進(jìn)行負(fù)載匹配，即確定最佳的傳動(dòng)比和發(fā)電機(jī)的參數(shù)，以達(dá)到最佳發(fā)電功率。當(dāng)使用水動(dòng)力軟件對浮子水動(dòng)力性能進(jìn)行計(jì)算時(shí)，只能在浮子上添加附加質(zhì)量和負(fù)載阻尼等類型的負(fù)載，但是在實(shí)際情況中，發(fā)電機(jī)是一種較為復(fù)雜的負(fù)載，其阻尼系數(shù)難以根據(jù)產(chǎn)品基本參數(shù)信息查得，因此難以直接使用水動(dòng)力軟件計(jì)算以發(fā)電機(jī)為負(fù)載的浮子水動(dòng)力響應(yīng)，除此之外還需要考慮齒輪箱傳動(dòng)比的影響。因此本文提出了一種考慮上述影響因素的負(fù)載匹配模型，能夠根據(jù)浮子快速確定最佳負(fù)載參數(shù)。

Ronald等[4-7]利用數(shù)學(xué)公式計(jì)算了漂浮圓柱體的水動(dòng)力系數(shù)和無量綱波浪力振幅。張萬超[8]等基于單浮體垂蕩動(dòng)力學(xué)模型得到了最優(yōu)PTO阻尼，并分析了PTO對垂向運(yùn)動(dòng)RAO和俘獲寬度比的影響，但是該動(dòng)力學(xué)模型并沒有考慮PTO的轉(zhuǎn)動(dòng)慣量的影響，并且PTO阻尼不能反映負(fù)載參數(shù)的影響。Felipe等[9-10]建立了雙體系統(tǒng)的動(dòng)力學(xué)模型，模型中包含了PTO的轉(zhuǎn)動(dòng)慣量并得到了PTO阻尼的表達(dá)式，但是該P(yáng)TO阻尼是根據(jù)與發(fā)電機(jī)相連的電阻所得到的，沒有充分考慮發(fā)電機(jī)自身性能參數(shù)的影響。李增亮等[11]研究了浮子式波能轉(zhuǎn)換裝置與發(fā)電機(jī)功率的最佳匹配，采用線性阻尼力代替直流發(fā)電機(jī)的特性，分析了波浪參數(shù)對匹配功率的影響規(guī)律以及變化負(fù)載對最大發(fā)電功率的影響。葉寅等[12-13]建立了波浪能裝置液壓轉(zhuǎn)換系統(tǒng)的仿真模型，該系統(tǒng)包括了各液壓元件與發(fā)電機(jī)，較好地模擬出了波浪輸入下液壓能量轉(zhuǎn)換系統(tǒng)的發(fā)電特性曲線。張?zhí)m勇等[14]對直流力矩電機(jī)進(jìn)行了理想建模，得到了直流電機(jī)的數(shù)學(xué)模型，同時(shí)也對電機(jī)的機(jī)械摩擦與內(nèi)部損耗等非線性因素進(jìn)行了分析。

本文利用Simulink搭建了機(jī)械式振蕩浮子式波浪能發(fā)電裝置的浮子垂蕩運(yùn)動(dòng)模型，模型中考慮了傳動(dòng)比與發(fā)電機(jī)的影響。根據(jù)發(fā)電機(jī)的性能參數(shù)，對不同直徑的浮子進(jìn)行負(fù)載匹配。除此之外，本文從Simulink模型中推導(dǎo)出能夠施加于水動(dòng)力軟件的附加質(zhì)量和負(fù)載阻尼，并將水動(dòng)力軟件的計(jì)算結(jié)果與Simulink的計(jì)算結(jié)果進(jìn)行了對比驗(yàn)證。所建模型可以應(yīng)用于機(jī)械式振蕩浮子式波浪能發(fā)電裝置的設(shè)計(jì)、機(jī)械結(jié)構(gòu)參數(shù)的確定和發(fā)電機(jī)的選型等，并對其他類型的波浪能發(fā)電裝置的設(shè)計(jì)提供借鑒。

1 Simulink模型推導(dǎo)

1.1 浮子垂蕩運(yùn)動(dòng)模型

采用機(jī)械傳動(dòng)的振蕩浮子式波浪能發(fā)電裝置的結(jié)構(gòu)如圖1所示，浮子與齒條固定在一起，齒條與齒輪嚙合傳動(dòng)，齒輪轉(zhuǎn)軸連接齒輪箱，齒輪箱連接發(fā)電機(jī)。發(fā)電機(jī)作為負(fù)載時(shí)，其具有阻尼，并且發(fā)電機(jī)轉(zhuǎn)子的轉(zhuǎn)動(dòng)慣量J會(huì)導(dǎo)致轉(zhuǎn)子轉(zhuǎn)動(dòng)時(shí)產(chǎn)生慣性轉(zhuǎn)矩。將發(fā)電機(jī)的阻尼系數(shù)與轉(zhuǎn)子的轉(zhuǎn)動(dòng)慣量折算到浮子上，分別得到負(fù)載阻尼CP和轉(zhuǎn)子等效質(zhì)量mJ，轉(zhuǎn)子等效質(zhì)量是添加在浮子上的附加質(zhì)量，由此得到浮子動(dòng)力學(xué)模型，如圖2所示，其中kw為靜水回復(fù)力剛度，μ33為垂蕩附加質(zhì)量，λ33為垂蕩輻射阻尼，F(xiàn)e為垂向波浪激勵(lì)力。

圖1 波浪能發(fā)電裝置結(jié)構(gòu)Fig.1 Structure diagram of wave energy converter

圖2 浮子動(dòng)力學(xué)模型Fig.2 Buoy dynamic model

本論文暫不考慮海流對機(jī)構(gòu)垂向運(yùn)動(dòng)的影響，因此根據(jù)牛頓第二定律可以得到浮子在入射波作用下的動(dòng)力學(xué)方程：

(1)

式中：m為浮子質(zhì)量；z為浮子垂向運(yùn)動(dòng)位移；Fr為垂向輻射力；Fe為垂向波浪激勵(lì)力；Fs為垂向靜水回復(fù)力；FP為負(fù)載力。

根據(jù)勢流理論，垂向輻射力Fr可分解為與浮子垂向運(yùn)動(dòng)加速度有關(guān)的附加質(zhì)量力和與浮子垂向運(yùn)動(dòng)速度有關(guān)的輻射阻尼力。只考慮一階作用力，垂向輻射力可表示為：

(2)

垂向波浪激勵(lì)力Fe可用無量綱波浪力振幅表示：

(3)

垂向靜水回復(fù)力Fs可表示為：

Fs=-kwz=-ρgAwz

(4)

式中：kw為靜水回復(fù)力剛度；Aw為水線面面積，對于圓柱型浮子Aw=πR2。

對于發(fā)電機(jī)系統(tǒng)，負(fù)載力FP可用負(fù)載阻尼CP和轉(zhuǎn)子等效質(zhì)量mJ表示：

(5)

將上述表達(dá)式整理可得浮子垂蕩運(yùn)動(dòng)方程：

(6)

對式(6)兩端進(jìn)行拉式變換，可得浮子垂蕩運(yùn)動(dòng)傳遞函數(shù)：

(7)

式中，M=m+μ33。

當(dāng)無負(fù)載時(shí)，mJ=CP=0，此時(shí)浮子做自由垂蕩運(yùn)動(dòng)，其傳遞函數(shù)為：

(8)

1.2 發(fā)電機(jī)負(fù)載方程

以直流電機(jī)作為發(fā)電環(huán)節(jié)中的發(fā)電機(jī)，電樞回路圖如圖3所示。轉(zhuǎn)子轉(zhuǎn)動(dòng)從而在電樞回路中產(chǎn)生電樞電流和反電動(dòng)勢。

圖3 電樞回路圖Fig.3 Armature circuit diagram

根據(jù)電路理論，列出電樞回路的電壓平衡方程：

(9)

電機(jī)感應(yīng)電勢方程為：

Ea=Ke·ωM

(10)

電磁轉(zhuǎn)矩基本方程為：

TF=Kt·Ia

(11)

式中：Ea為反電動(dòng)勢；Ia為電樞電流；La為電樞電感；Ra為電樞電阻；Ke為反電動(dòng)勢系數(shù)，ωM為轉(zhuǎn)子角速度；Kt為轉(zhuǎn)矩系數(shù)，TF為電磁轉(zhuǎn)矩。

將上述表達(dá)式進(jìn)行拉式變換并整理可得電磁轉(zhuǎn)矩傳遞函數(shù)：

(12)

由于機(jī)械摩擦、磁滯損耗以及渦流損耗等原因，在發(fā)電機(jī)中會(huì)產(chǎn)生阻轉(zhuǎn)矩Tf，但是阻轉(zhuǎn)矩Tf的影響較小，本論文中將此項(xiàng)忽略。由于轉(zhuǎn)子的轉(zhuǎn)動(dòng)慣量J，轉(zhuǎn)子轉(zhuǎn)動(dòng)時(shí)會(huì)產(chǎn)生慣性轉(zhuǎn)矩TJ。因此負(fù)載轉(zhuǎn)矩TP表示為：

TP=TF+TJ

(13)

1.3 傳動(dòng)機(jī)構(gòu)方程與轉(zhuǎn)子等效質(zhì)量

對于采用機(jī)械傳動(dòng)的振蕩浮子式波浪能發(fā)電裝置，其傳動(dòng)機(jī)構(gòu)為齒輪齒條和齒輪箱。齒輪軸與齒輪箱輸入軸相連，齒輪箱輸出軸與發(fā)電機(jī)轉(zhuǎn)軸相連。此傳動(dòng)機(jī)構(gòu)的方程為：

(14)

將轉(zhuǎn)子的慣性轉(zhuǎn)矩折算到浮子上，并與轉(zhuǎn)子等效質(zhì)量產(chǎn)生的慣性力相等，由此可得：

(15)

聯(lián)立式(14)和式(15)，可得轉(zhuǎn)子等效質(zhì)量為：

(16)

1.4 Simulink模型

結(jié)合1.1、1.2和1.3的內(nèi)容，可以建立浮子垂蕩運(yùn)動(dòng)的Simulink模型，如圖4所示。該模型中的“l(fā)ambda_33”為垂蕩輻射阻尼λ33；“gamma”為等效傳動(dòng)系數(shù)γ。

該模型在浮子自由垂蕩運(yùn)動(dòng)的基礎(chǔ)上，加入發(fā)電機(jī)負(fù)載與等效傳動(dòng)系數(shù)。浮子運(yùn)動(dòng)的速度經(jīng)過傳動(dòng)機(jī)構(gòu)后，轉(zhuǎn)換為發(fā)電機(jī)轉(zhuǎn)子轉(zhuǎn)動(dòng)的角速度。發(fā)電機(jī)的負(fù)載轉(zhuǎn)矩經(jīng)過傳動(dòng)機(jī)構(gòu)后，轉(zhuǎn)換為負(fù)載力作用在浮子上。輸入的垂向波浪激勵(lì)力需要克服負(fù)載力才能驅(qū)動(dòng)浮子進(jìn)行運(yùn)動(dòng)。反電動(dòng)勢Ea與感應(yīng)電流Ia的乘積為Simulink模型的理想瞬時(shí)發(fā)電功率。

圖4 仿真模型

1.5 系統(tǒng)傳遞函數(shù)的簡化

為了驗(yàn)證圖4中的Simulink模型計(jì)算結(jié)果，采用水動(dòng)力軟件計(jì)算相同負(fù)載下的浮子垂蕩運(yùn)動(dòng)響應(yīng)，因此需要根據(jù)Simulink模型來求解負(fù)載阻尼CP。

根據(jù)圖4可得系統(tǒng)閉環(huán)傳遞函數(shù)為：

(17)

式中：X=Laγ2

Y=Raγ2

E=JLa+MLaγ2

F=JRa+(MRa+λ33La)γ2

G=KeKt+(kwLa+λ33Ra)γ2

H=kwRaγ2

由式(17)可知，Simulink模型為三階系統(tǒng)，無法直接求出負(fù)載阻尼CP。式(7)為二階系統(tǒng)，此二階系統(tǒng)包含了負(fù)載阻尼CP。因此將三階系統(tǒng)化簡為二階系統(tǒng)。

不妨令：

(18)

對式(18)分母中的多項(xiàng)式乘積進(jìn)行展開，并令式(17)和式(18)分母中的各項(xiàng)系數(shù)對應(yīng)相等，可得：

XMg=E

(19)

CgX+MgY=F

(20)

KgX+CgY=G

(21)

KgY=H

(22)

由式(19)和式(22)可得：

(23)

(24)

由式(20)可得：

(25)

由式(21)可得：

(26)

因?yàn)槭?25)中只有垂蕩輻射阻尼，式(26)中除了垂蕩輻射阻尼，還包含有一項(xiàng)阻尼，所以選取式(26)為Cg的表達(dá)式，并定義：

(27)

綜上，得到二階系統(tǒng)表達(dá)式：

(28)

1.6 誤差分析

式(17)與式(28)相減，可得三階系統(tǒng)與二階系統(tǒng)的誤差傳遞函數(shù)E(s)：

(29)

對式(3)進(jìn)行拉式變換，得：

(30)

當(dāng)輸入為垂向波浪激勵(lì)力Fe(s)時(shí)，輸出為浮子垂向位移誤差：

ΔZ(s)=E(s)·Fe(s)

(31)

對上式進(jìn)行拉式逆變換，得到位移誤差響應(yīng)：

(32)

式中，σk為多項(xiàng)式Eσ3+Fσ2+Gσ+H的第k個(gè)根。

則相對位移誤差為

(33)

式中,z3為三階系統(tǒng)浮子垂蕩位移響應(yīng)。

2 結(jié)果與分析

2.1 浮子模型和仿真參數(shù)

本文中的計(jì)算模型為三種尺寸的圓柱型浮子，浮子參數(shù)如表1所示。在ANSYS AQWA中建立三種浮子的數(shù)值模型，網(wǎng)格最大單元尺寸分別為0.03 m、0.04 m、0.05 m，網(wǎng)格個(gè)數(shù)分別為15 388、23 036、27 754。2號浮子模型正視圖和網(wǎng)格示意圖如圖5。

表1 圓柱型浮子參數(shù)Tab.1 Parameters of cylindrical buoys

圖5 AQWA中2號浮子模型正視圖(左)和網(wǎng)格示意圖(右)Fig.5 Front view (left) and grid diagram (right) of buoy No.2 in AQWA

在AQWA中Hydrodynamic Diffraction模塊計(jì)算圓柱型浮子垂蕩運(yùn)動(dòng)的附加質(zhì)量、輻射阻尼和無量綱波浪力振幅，計(jì)算頻率為0.05～6.00 rad/s，頻率個(gè)數(shù)為98，水深為50 m。求出三種浮子在表2所示“威海國家淺海綜合試驗(yàn)場”海域平均海況下的水動(dòng)力系數(shù)，如表3所示。

表2 海況Tab.2 Sea state

表3 浮子水動(dòng)力系數(shù)Tab.3 Hydrodynamic coefficients of buoys

根據(jù)以上參數(shù)，計(jì)算得到浮子垂向波浪激勵(lì)力，如表4所示。

表4 浮子垂向波浪激勵(lì)力Tab.4 Vertical wave excitation force of buoys

由式(16)和式(27)可知，傳動(dòng)機(jī)構(gòu)中實(shí)際影響轉(zhuǎn)子等效質(zhì)量與負(fù)載阻尼的物理量是等效傳動(dòng)系數(shù)γ，為了便于研究與分析，取rg為0.1 m，通過改變傳動(dòng)比來改變等效傳動(dòng)系數(shù)。發(fā)電機(jī)的參數(shù)如表5所示。

表5 130LY51型直流電機(jī)性能參數(shù)Tab.5 Performance parameters of 130LY51 DC motor

2.2 二階系統(tǒng)輸出響應(yīng)誤差

在本文的仿真參數(shù)下，計(jì)算得到不同直徑浮子在不同傳動(dòng)比下的最大相對位移誤差。由表6可知，簡化后的二階系統(tǒng)的最大相對位移誤差均小于1%，誤差很小，因此可以將二階系統(tǒng)得出的負(fù)載阻尼與轉(zhuǎn)子等效質(zhì)量施加于AQWA中的垂蕩浮子上，用AQWA的計(jì)算結(jié)果驗(yàn)證Simulink模型。

表6 二階系統(tǒng)最大相對位移誤差Tab.6 Maximum relative displacement error of second order system

2.3 負(fù)載參數(shù)

根據(jù)式(16)和式(27)，計(jì)算得到不同傳動(dòng)比下的轉(zhuǎn)子等效質(zhì)量mJ與負(fù)載阻尼CP，如表7所示。在AQWA中Additional Added Mass和Additional Damping分別添加mJ和CP,只解放浮子的垂蕩自由度，在表2的海況下進(jìn)行時(shí)域計(jì)算，得到浮子的垂蕩運(yùn)動(dòng)響應(yīng)，并與Simulink模型浮子的運(yùn)動(dòng)響應(yīng)做對比。

表7 轉(zhuǎn)子等效質(zhì)量與負(fù)載阻尼Tab.7 Equivalent mass of rotor and load damping

2.4 驗(yàn)證與分析

經(jīng)Simulink模型計(jì)算得到三種浮子在不同傳動(dòng)比處的垂蕩振幅AB，如表8所示。

表8 浮子在不同傳動(dòng)比處的垂蕩振幅Tab.8 Heave amplitude of buoys at different transmission ratios

用式(34)對浮子垂蕩運(yùn)動(dòng)的振幅進(jìn)行無量綱化，用式(35)計(jì)算AQWA浮子瞬時(shí)發(fā)電功率，用式(36)計(jì)算Simulink浮子瞬時(shí)發(fā)電功率。用式(37)計(jì)算平均發(fā)電功率PB，平均發(fā)電功率為瞬時(shí)發(fā)電功率PI在一個(gè)波浪周期內(nèi)的平均值。

(34)

式中，AB為浮子垂蕩運(yùn)動(dòng)的振幅。

AQWA瞬時(shí)發(fā)電功率計(jì)算式為：

(35)

Simulink瞬時(shí)發(fā)電功率計(jì)算式為：

(36)

平均發(fā)電功率計(jì)算式為：

(37)

式中，T為入射波周期。

1號、2號和3號浮子在不同傳動(dòng)比處的無量綱垂蕩振幅和平均發(fā)電功率如圖6、7和8所示。從上述對比圖中可以看出，Simulink模型的計(jì)算結(jié)果與AQWA的計(jì)算結(jié)果吻合很好，隨著傳動(dòng)比倒數(shù)的增大，無量綱垂蕩振幅近似呈現(xiàn)線性衰減，平均發(fā)電功率先增大后減小，出現(xiàn)了最大平均發(fā)電功率，并且同一浮子的兩種計(jì)算結(jié)果在相同的傳動(dòng)比處達(dá)到最大平均發(fā)電功率，因此Simulink模型反映出了無量綱垂蕩振幅和平均發(fā)電功率隨傳動(dòng)比變化的趨勢，并能準(zhǔn)確預(yù)測出最佳發(fā)電功率點(diǎn)。

表9為Simulink模型與AQWA的計(jì)算誤差。浮子垂蕩振幅的計(jì)算誤差不超過6 mm，對應(yīng)的相對誤差小于3%，平均發(fā)電功率幅值的相對誤差小于4%。

(a) 浮子無量綱垂蕩振幅

(a) 浮子無量綱垂蕩振幅

(a) 浮子無量綱垂蕩振幅

表9 計(jì)算誤差Tab.9 Calculation errors

圖9至圖14為不同浮子在最佳傳動(dòng)比(最大發(fā)電功率處的傳動(dòng)比)處的垂蕩位移和瞬時(shí)發(fā)電功率。從圖中可以看出Simulink模型很好地反映出了浮子的垂蕩運(yùn)動(dòng)。Simulink瞬時(shí)發(fā)電功率的計(jì)算結(jié)果是：在一個(gè)波浪周期內(nèi)，出現(xiàn)兩個(gè)峰值，并且兩個(gè)峰值大小相等。但是AQWA瞬時(shí)發(fā)電功率的計(jì)算結(jié)果是：在一個(gè)波浪周期內(nèi)出現(xiàn)的兩個(gè)峰值大小不等，但是二者差值較小。這是因?yàn)樵赟imulink模型中，浮子垂蕩運(yùn)動(dòng)的上行和下行最大垂蕩速度相等；而在AQWA中浮子垂蕩運(yùn)動(dòng)的上行和下行最大垂蕩速度有微小差異，由于負(fù)載阻尼的數(shù)量級遠(yuǎn)大于浮子垂蕩速度的數(shù)量級，因此這種微小差異在計(jì)算瞬時(shí)發(fā)電功率時(shí)被放大，產(chǎn)生了AQWA中一個(gè)波浪周期內(nèi)出現(xiàn)的兩個(gè)峰值大小不等的現(xiàn)象。

3 結(jié) 論

(1) Simulink與AQWA的計(jì)算結(jié)果吻合很好，充分反映了浮子垂蕩運(yùn)動(dòng)振幅和平均發(fā)電功率幅值隨負(fù)載變化的趨勢，相對誤差不超過4%。

(2) 本文中發(fā)電機(jī)的負(fù)載模型是根據(jù)發(fā)電機(jī)的基本性能參數(shù)和電磁理論推導(dǎo)計(jì)算得出，而不是單純的阻尼參數(shù)，更貼合實(shí)際。

(3) Simulink與AQWA的計(jì)算結(jié)果均在相同的最佳負(fù)載阻尼下達(dá)到最大平均發(fā)電功率。在本文中，負(fù)載阻尼是根據(jù)傳動(dòng)比和發(fā)電機(jī)負(fù)載模型計(jì)算得到，能夠用于確定最佳傳動(dòng)比。

(4) 最佳傳動(dòng)比與齒輪半徑的乘積即為最佳等效傳動(dòng)系數(shù)，因此可以在保證最佳等效傳動(dòng)系數(shù)不變的前提下，根據(jù)設(shè)計(jì)需求靈活選擇不同的齒輪半徑rg，并得到對應(yīng)的最佳傳動(dòng)比。

(5) 所建模型可以應(yīng)用于振蕩浮子式波浪能發(fā)電裝置的設(shè)計(jì)、機(jī)械結(jié)構(gòu)參數(shù)的確定和發(fā)電機(jī)的選型。