王驥源，賀冠翔，趙書玲

（江西理工大學(xué)電氣工程與自動化學(xué)院，江西 贛州 341000）

波浪能是一種重要的海洋可再生能源，分布廣泛，儲量豐富，具有可觀的應(yīng)用前景，如何將波浪能最大化地轉(zhuǎn)化利用成為新能源領(lǐng)域的熱門話題。波浪能轉(zhuǎn)換裝置有液壓式、氣壓式、機械式3 種形式，李燕等[1]研究表明，液壓式波浪能轉(zhuǎn)換裝置由于轉(zhuǎn)換效率高、輸出穩(wěn)定，為首選波浪能轉(zhuǎn)換裝置。液壓式裝置又可細(xì)分為點頭鴨式、搖擺式、閥式、振浮式4 種，其中振浮式應(yīng)用研究最廣，振浮式波浪能裝置主要通過垂蕩和縱搖獲取并轉(zhuǎn)換波浪能[2]。對于此類裝置輸出功率的優(yōu)化，國內(nèi)做了不少研究，已有的相關(guān)研究中，研究人員通常借助軟件仿真、算法優(yōu)化等方式找出相應(yīng)的最大輸出功率。潘海鵬等[3]基于遺傳算法對以浮子為驅(qū)動裝置的直驅(qū)式海浪發(fā)電系統(tǒng)的最大功率進行跟蹤控制，此方法可大幅度提升發(fā)電系統(tǒng)對波浪能的捕獲效率；黃俊豪等[4]借助WFT 法獲取波浪能發(fā)電系統(tǒng)最大功率的捕獲條件，并結(jié)合Simulink 仿真，最終提升了使得系統(tǒng)功率輸出的優(yōu)化效果；楊俊華等[5]根據(jù)商用軟件求解浮子的非線性參數(shù)，結(jié)合改進型的鯨魚算法，優(yōu)化了相應(yīng)負(fù)載參數(shù)，進而提高了波浪能捕獲效率。

本文根據(jù)振浮式波浪能轉(zhuǎn)換裝置垂蕩和縱搖這2 種運動方式，并建立該裝置運動和功率模型，在此基礎(chǔ)上進行數(shù)值求解和算法優(yōu)化，設(shè)計出了實現(xiàn)波浪能裝置輸出功率最大化的最優(yōu)系統(tǒng)參數(shù)，為實際應(yīng)用提供參考。

1 問題的提出及理論依據(jù)

圖1 為振浮式波浪能裝置示意圖，由浮子、振子、中軸及能量輸出系統(tǒng)（簡稱“PTO 系統(tǒng)”，包括彈簧、直線阻尼器、扭轉(zhuǎn)彈簧、旋轉(zhuǎn)阻尼器和轉(zhuǎn)軸）等構(gòu)成。在波浪激勵力和波浪激勵力矩作用下，浮子同時進行垂蕩和縱搖，并通過PTO 系統(tǒng)使振子沿著中軸進行垂蕩和縱搖。除此之外，浮子在線性周期微幅波作用下[6]，還受到附加慣性力和附加慣性力矩、興波阻尼力和興波阻尼力矩、靜水恢復(fù)力和靜水恢復(fù)力矩的影響。振子沿中軸垂蕩時，直線阻尼器的阻尼力與浮子和振子的相對速度成正比，比例系數(shù)為直線阻尼器的直線阻尼系數(shù)；振子沿中軸縱搖時，旋轉(zhuǎn)阻尼器的扭矩與浮子和振子的相對角速度成正比，比例系數(shù)為旋轉(zhuǎn)阻尼器的旋轉(zhuǎn)阻尼系數(shù)，2 個阻尼系數(shù)在[0，100 000]區(qū)間取值。浮子和振子通過相對垂蕩運動和相對縱搖運動驅(qū)動直線阻尼器和旋轉(zhuǎn)阻尼器做功，并將所做的功作為該波浪能裝置的能量輸出。結(jié)合相關(guān)條件，計算出波浪能裝置輸出的最大功率，并確定該功率下直線阻尼器的阻尼系數(shù)和旋轉(zhuǎn)阻尼器的阻尼系數(shù)。分析問題時，考慮海水是無粘無旋的，忽略中軸、PTO 的質(zhì)量和各種摩擦。

圖1 振浮式波浪能裝置示意圖

遺傳算法最早由美國的JOHO 于20 世紀(jì)70 年代提出，該算法起源于對生物系統(tǒng)所進行的計算機模擬研究，是一種通過模擬自然進化過程搜索最優(yōu)解的方法。該算法通過數(shù)學(xué)的方式，利用計算機仿真運算將問題的求解過程轉(zhuǎn)換成生物進化中染色體基因的交叉、變異等過程。它的大致實現(xiàn)過程為從任意初始種群出發(fā)，通過隨機選擇、交叉和變異操作[7]，產(chǎn)生一群更適合環(huán)境的個體，使群體進化到搜索空間中越來越好的區(qū)域，這樣多代繁衍進化，最后收斂到一群最適應(yīng)環(huán)境的個體，求得問題最優(yōu)解。

2 數(shù)值模型

2.1 垂蕩運動狀態(tài)分析

裝置受力情況如下。

振子進行垂蕩時，只受PTO 系統(tǒng)力，而PTO 系統(tǒng)力包含彈簧力和直線阻尼力，具體受力情況如圖2 所示。

圖2 振子受力圖

振子的受力方程如下：

式中：F1為振子所受合力的數(shù)值；FPTO為PTO 系統(tǒng)力的數(shù)值；Fc為彈簧力的數(shù)值；Fx為旋轉(zhuǎn)阻尼器的阻尼力的數(shù)值。

結(jié)合上文可知，浮子受力情況如圖3 所示。

圖3 浮子受力圖

波浪激勵力可表示為：

式中：f為波浪激勵力振幅的數(shù)值；ω為入射波浪頻率的數(shù)值。

結(jié)合浮子的受力情況，并根據(jù)牛頓第二定律，可以列出浮子的垂蕩運動方程為：

式中：mb為浮子質(zhì)量的數(shù)值；m′為垂蕩附加質(zhì)量的數(shù)值；為浮子垂蕩運動加速度的數(shù)值；Fb為興波阻尼力的數(shù)值；Fl為靜水恢復(fù)力的數(shù)值。

同理，振子的垂蕩運動方程為：

式中：Cz為直線阻尼系數(shù)；為浮子垂蕩運動速度的數(shù)值；為振子垂蕩運動速度的數(shù)值；Kc為彈簧剛度的數(shù)值；xb（t）為浮子垂蕩位移的數(shù)值；xv（t）為振子垂蕩位移的數(shù)值。

2.2 縱搖運動狀態(tài)分析

裝置所受力矩如下。

類比垂蕩，振子進行縱搖時，受到的力矩如下：

式中：M1為振子所受力矩的數(shù)值；MPTO為PTO 系統(tǒng)力矩的數(shù)值；Mc為彈簧扭矩的數(shù)值；Mx為旋轉(zhuǎn)阻尼器的阻尼扭矩的數(shù)值。

結(jié)合上文可知浮子受到的力矩，并根據(jù)剛性力學(xué)中的轉(zhuǎn)動定律，可以列出浮子的縱搖運動方程為：

式中：Ib為浮子轉(zhuǎn)動慣量的數(shù)值；I′ 為縱搖附加轉(zhuǎn)動慣量的數(shù)值；為浮子角加速度的數(shù)值；Mw為興波阻尼力矩的數(shù)值；Mj為靜水恢復(fù)力矩的數(shù)值。

振子的縱搖運動方程為：

式中：Iv為振子轉(zhuǎn)動慣量的數(shù)值；為振子角加速度的數(shù)值；Kx為旋轉(zhuǎn)彈簧剛度的數(shù)值；θb（t）為浮子角位移的數(shù)值；θv（t）為振子角位移的數(shù)值；Cx為旋轉(zhuǎn)阻尼系數(shù)的數(shù)值；為浮子角速度的數(shù)值；為振子角速度的數(shù)值。

2.3 波浪能輸出功率

該裝置垂蕩時，由直線阻尼器和彈簧組成的PTO系統(tǒng)對外做功，因為該裝置并非是受到恒力運動，結(jié)合微元法和積分法，給出相應(yīng)的功率公式如下：

該裝置縱搖時，由旋轉(zhuǎn)阻尼器和扭轉(zhuǎn)彈簧組成的PTO 系統(tǒng)對外做功，該裝置縱搖時受到的力矩亦非平衡，同理，輸出功率公式如下：

結(jié)合式（1）和式（2）可知，該裝置對外輸出的總功率如下：

3 數(shù)值計算和結(jié)果分析

3.1 轉(zhuǎn)動慣量的確定

在求解最大功率前，為確定裝置縱搖運動對應(yīng)的常量參數(shù)，需先求出浮子和振子的轉(zhuǎn)動慣量，求解過程如下。

先建立浮子的直角坐標(biāo)系，如圖4 所示。

再確定浮子質(zhì)心。分別確定圓柱和圓錐對應(yīng)質(zhì)心坐標(biāo)Z1和Z2，結(jié)合質(zhì)心運動定理，可得浮子的質(zhì)心坐標(biāo)：

然后確定浮子面密度：

式中：ρb為浮子面密度的數(shù)值；mb為浮子質(zhì)量的數(shù)值；L圓柱為圓柱半徑的數(shù)值；L圓錐為圓錐半徑的數(shù)值。

圓柱側(cè)面和圓錐側(cè)面轉(zhuǎn)動慣量計算如下：

式中：I1為圓柱側(cè)面轉(zhuǎn)動慣量的數(shù)值；I2為圓錐側(cè)面轉(zhuǎn)動慣量的數(shù)值。

最后，結(jié)合平行軸定理[8]，可求得圓柱側(cè)面和圓錐側(cè)面對于質(zhì)心軸的轉(zhuǎn)動慣量之和，即浮子的轉(zhuǎn)動慣量，公式如下：

式中：d1為圓柱質(zhì)心到浮子質(zhì)心距離的數(shù)值；d2為圓柱質(zhì)心到浮子質(zhì)心距離的數(shù)值。

同理，可求得振子轉(zhuǎn)動慣量：

3.2 最大輸出功率的計算

求得轉(zhuǎn)動慣量后，結(jié)合浮子振子的運動模型和輸出功率的模型，引入波浪能裝置的常量參數(shù)（如表1所示）和直線阻尼系數(shù)、旋轉(zhuǎn)阻尼系數(shù)這2 個變量參數(shù)，編寫優(yōu)化算法，求解該波浪能裝置的最大輸出功率。

表1 波浪能裝置的常量參數(shù)

算法步驟如下：①給出求解浮子振子運動狀態(tài)的算法框架1。將浮子和振子垂蕩和縱搖2 個運動對應(yīng)的微分方程模型導(dǎo)入Matlab 的ode45 函數(shù)中，設(shè)定浮子和振子的初始位移和速度均為0，運動總時間為184.2 s，步長為0.2 s，最終可求得浮子和振子各自在0～184.2 s 內(nèi)每隔0.2 s 的速度，由于尚未代入變量值，暫時無法求解。②在框架1 的基礎(chǔ)上，結(jié)合式（3）和式（4），先離散化計算每隔0.2 s 波浪能裝置對外所做的功，再通過數(shù)值積分計算出0～184.2 s 所做出的總功，并除以總周期，最終給出求解輸出功率所需要的算法框架2。③設(shè)定GA 工具箱的初始參數(shù)，具體設(shè)定為lb=[0，0]，表示2 個變量下限；ub=[100 000，100 000]，表示2 個變量上限；PopulationSize=10，表示種群數(shù)量；Tolfun=1e-3，表示允許誤差；CrossoverFraction=0.75，表示交叉率；EliteCount=2，表示每次遺傳中的篩選個數(shù)；Generation=20，表示迭代遺傳次數(shù)。④將步驟②中編寫好的算法框架2 導(dǎo)入設(shè)置好參數(shù)的GA 工具箱，并開始運行，最終求得結(jié)果。

3.3 結(jié)果分析

GA 收斂圖如圖5 所示。根據(jù)圖5 可以看到GA 算法的收斂效果較好。

圖5 GA 收斂圖

波浪能裝置輸出的最大功率最終收斂于315.77 W，相應(yīng)的最優(yōu)直線阻尼系數(shù)為589 00 s/m，最優(yōu)旋轉(zhuǎn)阻尼系數(shù)為95 541 s/m。

為了驗證結(jié)果的可行性，以時間為自變量、位移為因變量，給出浮子振子在最優(yōu)直線阻尼系數(shù)下的垂蕩運動狀態(tài)圖，如圖6 所示。

圖6 垂蕩運動圖

由圖6 不難看出浮子和振子垂蕩狀態(tài)類似波浪進行簡諧運動的狀態(tài)。除了垂蕩，再以時間為自變量、角位移為因變量，給出浮子振子在最優(yōu)旋轉(zhuǎn)阻尼系數(shù)下的縱搖狀態(tài)圖，如圖7 所示。

圖7 縱搖運動圖

縱搖運動也類似波浪的簡諧運動。結(jié)合遺傳算法收斂圖和裝置運動狀態(tài)圖，最終可以較好地驗證該最優(yōu)系統(tǒng)參數(shù)的可行性，并與實際情況相符合。

4 結(jié)束語

本研究結(jié)合運動學(xué)和能量轉(zhuǎn)換等理論建立振浮式波浪能裝置輸出功率的數(shù)學(xué)模型，利用Matlab 編寫智能優(yōu)化算法，對實現(xiàn)裝置輸出功率最大化的最優(yōu)變量參數(shù)進行搜索。經(jīng)Matlab 的可視化檢驗，最終得出切合實際情況的最大輸出功率和2 個最優(yōu)阻尼系數(shù)，為實際推廣應(yīng)用提供一定參考。

但此研究仍有不足之處，如在實際情況中，波高與風(fēng)力的影響不能忽略不計，尤其是在縱搖過程中，風(fēng)力矩帶來的影響不可忽略。其次，垂蕩與縱搖具有強耦合性，即垂蕩必然伴隨著縱搖，發(fā)生縱搖的同時一定會發(fā)生垂蕩。本文并未分析垂蕩縱搖的強耦合性對運動的影響，所以在今后研究中有待進一步改進。