基于MLP方法的船舶大幅橫搖近似解析解*

劉亞沖 胡安康，2 韓鳳磊 汪春輝

(1.哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱150001;2.中集船舶海洋工程設(shè)計(jì)研究院，上海201206)

在長(zhǎng)期的研究過程中研究者發(fā)現(xiàn)利用線性模型來描述船舶的橫搖運(yùn)動(dòng)與實(shí)際情況相差甚遠(yuǎn)，且線性模型無法得到非線性所具有的一些諸如跳躍、分岔、混沌等現(xiàn)象［1-4］.由于船體幾何外形原因，在船舶發(fā)生大角度橫傾時(shí)船舶橫搖的恢復(fù)力矩本身即具有非線性，而由于真實(shí)流體具有粘性，橫搖阻尼也具有非線性特征.因此研究船舶的非線性橫搖運(yùn)動(dòng)，對(duì)非線性橫搖運(yùn)動(dòng)進(jìn)行解析求解具有重要的實(shí)際意義.

對(duì)船舶非線性橫搖阻尼的處理，目前有線性項(xiàng)加平方項(xiàng)(LPQD)方式即和線性項(xiàng)加立方項(xiàng)(LPCD)方式即.Bikdash 等［5］對(duì)不同的阻尼模式進(jìn)行了研究分析.從計(jì)算分析角度來看，LPCD 模式無疑方便處理，因?yàn)樵撟枘崮Ｊ奖苊饬艘蚪^對(duì)值的存在給數(shù)學(xué)分析帶來的困難.船舶恢復(fù)力矩是橫搖角φ 的函數(shù)，考慮到船舶關(guān)于中線面的對(duì)稱性，可以用奇數(shù)次的高階多項(xiàng)式函數(shù)來擬合，那么非線性恢復(fù)力矩可以表示為M(φ)= C1φ + C3φ3+ C5φ5+ …，通常取前兩項(xiàng)即可很好地模擬靜穩(wěn)性曲線.

傳統(tǒng)的L-P 方法對(duì)于弱非線性系統(tǒng)的求解是可行的，若將其應(yīng)用于強(qiáng)非線性系統(tǒng)則會(huì)帶來較大誤差.文中采用改進(jìn)的L-P 法(即MLP 方法)對(duì)船舶的非線性橫搖問題進(jìn)行求解，分別得到靜水中和規(guī)則波中的近似解析解，并用數(shù)值計(jì)算來驗(yàn)證該方法的正確性.

1 船舶非線性橫搖運(yùn)動(dòng)方程的建立

考慮阻尼力矩和恢復(fù)力矩的非線性，船舶的非線性橫搖運(yùn)動(dòng)方程可以表述為:

式中:Jφφ、ΔJφφ為轉(zhuǎn)動(dòng)慣量和附加轉(zhuǎn)動(dòng)慣量;D1為線性阻尼力矩系數(shù)，可以通過勢(shì)流理論求得;D3為非線性阻尼力矩系數(shù)，可以通過橫搖衰減試驗(yàn)［6-7］或CFD 仿真［8-9］得到;C1為線性恢復(fù)力矩系數(shù)，C1=D，D 為船舶排水量，為船舶初穩(wěn)心高;C3為非線性恢復(fù)力矩系數(shù);C1和C3可以由船舶靜穩(wěn)性曲線得到.式(1)右端表示波浪激勵(lì)項(xiàng)［10］，F(xiàn) = D·GMkξ，k 為波數(shù)，ξ 為波高，ω 為波浪圓頻率.式(1)兩邊同除以(Jφφ+ ΔJφφ)后，可將式(1)寫為

式中:2vφφ為橫搖衰減系數(shù)，2vφφ= D1/(Jφφ+ ΔJφφ);ω0為船舶橫搖近似固有圓頻率，ω20= C1/(Jφφ+ΔJφφ);d3= D3/(Jφφ+ΔJφφ);c3= C3/(Jφφ+ΔJφφ);d =D/(Jφφ+ΔJφφ)，引入?yún)?shù)ε，為便于表示，設(shè)，將式(2)表達(dá)為擬線性系統(tǒng)的形式:

2 改進(jìn)的L-P(MLP)法

對(duì)于非線性動(dòng)力方程(3)的求解，攝動(dòng)方法是進(jìn)行近似解析求解的行之有效的方法.常規(guī)的攝動(dòng)方法如直接攝動(dòng)法、LP 法、平均法、KBM 法和多尺度法等只能對(duì)弱非線性的動(dòng)力方程進(jìn)行近似解析［11］.當(dāng)方程(3)中的參數(shù)不能滿足ε ?1 時(shí)，針對(duì)弱非線性系統(tǒng)的攝動(dòng)法會(huì)帶來求解誤差，鑒于此，文中采用改進(jìn)的L-P 方法(Modified L-P 法，簡(jiǎn)稱MLP方法［12］)對(duì)式(2)進(jìn)行求解.

2.1 自由橫搖的二階近似解

首先求解船舶自由橫搖的二階近似解析解，令方程(2)右端激勵(lì)項(xiàng)為零即可，此時(shí)橫搖方程為

方程(4)中的ε 不再限定為小參數(shù).引入新的自變量 = ωt，ω 代表待求的原系統(tǒng)的非線性頻率，于是式(4)變?yōu)?/p>

式中，φ'表示φ 對(duì) 的一階微商，φ″表示φ 對(duì) 的二階微商.將ω2在ω0附近展開為ε 的冪級(jí)數(shù)，即

然后，引入一個(gè)參數(shù)變換

于是總有0 ＜α ＜1 成立，并且

將式(8)和(9)代入方程(5)得

式中，α、δ2為待定的未知常數(shù).將φ 展開成新參數(shù)α 的冪級(jí)數(shù)，對(duì)自由橫搖取前兩階近似解有

將式(11)代入式(10)，令兩端α 的同次冪相等，得:

考慮初始化條件φ0(0)= a0，φi(0)= 0 (i =1，2)，φj'(0)= 0 (j = 0，1，2)，式(12)的解為φ0()= a0cos ，將φ0()代入式(13)，消去久期項(xiàng)得

于是，新參數(shù)為

將式(15)、(16)代入方程(14)，根據(jù)久期項(xiàng)為0 的條件可以求出δ2和φ2().

最終，求得精確至o(α3)的解為

2.2 規(guī)則波中的一階近似解

規(guī)則波中船舶非線性橫搖方程可寫成

同樣，引入新的自變量 = ωt，式(19)成為

將ω2在ω0附近展開為ε 的冪級(jí)數(shù)，并將φ 展開成α的冪級(jí)數(shù).將ω2和φ 代入式(19)，令兩端α 的同次冪相等得:

此時(shí)，初始化條件為φ0(0)=a0;φ1(0)=φ1'(0)=0，式(21)的解為φ0()= a0cos +Csin ，C 為待定常數(shù)，可在下面的式(23)中求出.將φ0()代入式(22)，得到久期項(xiàng)為0 的條件為

根據(jù)式(23)，可求出ω1和C，見式(24)和(25).

3 數(shù)值算法

在文中，數(shù)值算法選擇高精度單步算法四階Runge-Kutta 法，對(duì)于高階微分方程式(19)，首先要將其降階化為一階微分方程組.不妨選取φ = φ1，，于是有

若采用向量的記號(hào)，記φ = (φ1，φ2)T，f =(f1，f2)T，并考慮初值條件，則有式(27)成立.

求解這一初值問題的四階Runge-Kutta 公式為

式中，

根據(jù)式(28)和(29)，考慮到初始條件后就可以編制相應(yīng)的計(jì)算程序進(jìn)行時(shí)間步進(jìn)求解.

4 算例分析

為了便于分析比較，文中選取Gaul 拖船［13］作為算例.分別采用MLP 法和數(shù)值算法對(duì)算例進(jìn)行計(jì)算，以驗(yàn)證MLP 法的有效性.該拖船的主要特征參數(shù)見表1.

表1 算例參數(shù)Table1 Parameters of example

首先選取靜水工況，不妨選取10°的初始橫傾角，即φ(0)= 0.175 rad，采用數(shù)值算法和采用MLP得到的橫搖運(yùn)動(dòng)的時(shí)歷圖與Poincare 截面見圖1.

圖1表明采用MLP 法獲得的二階近似解與采用Runge-Kutta 法得到的數(shù)值解吻合較好.由于MLP法獲得的橫搖近似解只精確至二階，與更高階有關(guān)的橫搖頻率的信息無法體現(xiàn)，因此近似解析解與數(shù)值解之間有一些細(xì)微的差別，在時(shí)歷圖中則表現(xiàn)為沿時(shí)間軸方向有一些偏移.

對(duì)于規(guī)則波中情形，為使選取的波浪參數(shù)不失一般性，文中在Gaul 船的橫搖固有頻率附近選擇3個(gè)不同的波浪激勵(lì)頻率，根據(jù)深水色散關(guān)系可以求出波浪參數(shù)的波長(zhǎng)，根據(jù)波陡概率密度分布函數(shù)［14-15］，可以得到3 個(gè)波浪頻率對(duì)應(yīng)的波高.由此得到的3 個(gè)波浪參數(shù)見表2.

圖1 自由橫搖的時(shí)歷圖與Poincare 截面Fig.1 Free-roll time-domain diagram and Poincare section

表2 波浪參數(shù)Table2 Wave parameters

分別對(duì)表2表示的3 種波浪激勵(lì)進(jìn)行數(shù)值求解，并根據(jù)MLP 進(jìn)行攝動(dòng)求解，得到的結(jié)果見圖2.其中圖2(a)、2(c)和2(d)為兩種方法得到的時(shí)歷圖，圖2(b)為相圖.

圖2 規(guī)則波作用下的時(shí)歷圖與相圖Fig.2 Time-domain diagram under regular wave

從以上3 個(gè)時(shí)歷圖中可以看出數(shù)值解在前200 s的時(shí)間內(nèi)有波動(dòng)現(xiàn)象，然后逐步穩(wěn)定到穩(wěn)態(tài)的周期運(yùn)動(dòng);從圖2(b)相圖也可以看出，對(duì)于第1 組波浪參數(shù)，從第40 個(gè)周期開始橫搖就已經(jīng)是穩(wěn)定的周期運(yùn)動(dòng)，而MLP 求得的是系統(tǒng)穩(wěn)定橫搖階段的解.從波浪的激勵(lì)頻率角度來看，當(dāng)激勵(lì)頻率與橫搖派生系統(tǒng)的固有頻率接近時(shí)，橫搖的運(yùn)動(dòng)幅度達(dá)到最大，此時(shí)對(duì)應(yīng)系統(tǒng)的主共振情形.

5 結(jié)論

文中采用MLP 法對(duì)船舶非線性橫搖運(yùn)動(dòng)進(jìn)行近似解析求解，分別得到了靜水中自由衰減運(yùn)動(dòng)的二階近似解析解和在不同波浪參數(shù)規(guī)則波作用下的一階近似解析解.通過將解析解與Runge-Kutta 法得到的數(shù)值解進(jìn)行對(duì)比分析，驗(yàn)證了MLP 法求得的近似解的正確性.由于近似解析解略掉了與高階項(xiàng)有關(guān)的固有頻率，因而與數(shù)值解有細(xì)微的差別.

在MLP 方法中，通過引入?yún)?shù)變換，將ω2在展開可以使得對(duì)于任何的εω1，總能使新參數(shù)α 滿足0 ＜α ＜1，這樣就可以將原來相對(duì)于ε是強(qiáng)非線性的系統(tǒng)轉(zhuǎn)化為相對(duì)α 是弱非線性系統(tǒng)，這對(duì)于研究船舶的大幅橫搖運(yùn)動(dòng)和進(jìn)行橫搖穩(wěn)性評(píng)估具有重要意義.

［1］Nayfeh A H，Sanchez N E.Stability and complicated rolling responses of ship in regular beam seas［J］.International Shipbuilding Progress，1990，37(410):177-198.

［2］Taylan M.Static and dynamic aspects of a capsize phenomenon［J］.Ocean Engineering，2003，30(3):331-350.

［3］Francescutto A.Intact ship stability:the way ahead［J］.Marine Technology，2004，41(1):31-37.

［4］李遠(yuǎn)林.船舶非線性大幅橫搖導(dǎo)致傾覆的研究［J］.華南理工大學(xué)學(xué)報(bào):自然科學(xué)版，2001，29(6):91-96.Li Yuan-lin.Research in nonlinear large-amplitude roll leading to capsize in wind-waves ［J］.Journal of South China University of Technology:Natural Science Edition，2001，29(6):91-96.

［5］Bikdash M，Balachandran B，Nafey A H.Melnikov analysis for a ship with a general roll-damping model［J］.Nonlinear Dynamics，1994，6(1):101-124.

［6］李紅霞，唐友剛，胡楠，等.船舶橫搖非線性阻尼系數(shù)的識(shí)別［J］.天津大學(xué)學(xué)報(bào)，2005，38(12)1042-1045.Li Hong-xia，Tang You-gang，Hu Nan，et al.Identification for Nonlinear Ship Roll DampingCoefficients［J］.Journal of Tianjin University，2005，38(12):1042-1045.

［7］李遠(yuǎn)林，伍曉榕.非線性橫搖阻尼的試驗(yàn)確定-數(shù)據(jù)處理方法［J］.華南理工大學(xué)學(xué)報(bào):自然科學(xué)版，2002，30(2):79-82.Li Yuan-lin，Wu Xiao-rong.Experimental determination of nonlinear roll damping:a technique for data processing［J］.Journal of South China University of Technology:Natural Science Edition，2002，30(2):79-82.

［8］朱仁傳，郭海強(qiáng)，繆國(guó)平，等.一種基于CFD 理論船舶附加質(zhì)量與阻尼的計(jì)算方法［J］.上海交通大學(xué)學(xué)報(bào)，2009，42(2):198-203.Zhu Ren-chuan，Guo Hai-qiang，Miao Guo-ping，et al.A computational method for evaluation of added mass and damping of ship based on CFD theory ［J］.Journal of Shanghai Jiao Tong University，2009，42(2):198-203.

［9］黃昊，郭海強(qiáng)，朱仁傳，等.粘性流中船舶橫搖阻尼計(jì)算［J］.船舶力學(xué)，2008，12(4):568-573.Huang Hao，Guo Hai-qiang，Zhu Ren-chuan，et al.Computations of ship roll damping in viscousflow［J］.Journal of Ship Mechanics，2008，12(4):568-573.

［10］李積德.船舶耐波性［M］.哈爾濱:哈爾濱工程大學(xué)出版社，2001.

［11］胡海巖.應(yīng)用非線性動(dòng)力學(xué)［M］.北京:航空工業(yè)出版社，2000.

［12］陳樹輝.強(qiáng)非線性振動(dòng)系統(tǒng)的定量分析方法［M］.北京:科學(xué)出版社，2007.

［13］Morrall A.The gaul disaster:a investigation into the loss of a large stern trawler［J］.Naval Architect，1981(5):391-440.

［14］高志一，文凡，李潔.波群內(nèi)單個(gè)波的波陡分布與破碎波［J］.海洋科學(xué)，2011，35(9):96-106.Gao Zhi-yi，Wen Fan，Li Jie.The steepness distributions and break of individual waves in wavegroups［J］.Marine Science，2011，35(9):96-106.

［15］裴玉華，鄭桂珍，叢培秀.海浪的波陡分布［J］.中國(guó)海洋大學(xué)學(xué)報(bào)，2007，37(7):73-77.Pei Yu-hua，Zheng Gui-zhen，Cong Pei-xiu.The probability distribution function of wave steepness ［J］.Periodical of Ocean University of China，2007，37(7):73-77.