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

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

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

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

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

1 船舶非線性橫搖運動方程的建立

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

式中:Jφφ、ΔJφφ為轉(zhuǎn)動慣量和附加轉(zhuǎn)動慣量;D1為線性阻尼力矩系數(shù)，可以通過勢流理論求得;D3為非線性阻尼力矩系數(shù)，可以通過橫搖衰減試驗［6-7］或CFD 仿真［8-9］得到;C1為線性恢復(fù)力矩系數(shù)，C1=D，D 為船舶排水量，為船舶初穩(wěn)心高;C3為非線性恢復(fù)力矩系數(shù);C1和C3可以由船舶靜穩(wěn)性曲線得到.式(1)右端表示波浪激勵項［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)表達為擬線性系統(tǒng)的形式:

2 改進的L-P(MLP)法

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

2.1 自由橫搖的二階近似解

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

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

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

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

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

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

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

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

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

于是，新參數(shù)為

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

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

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

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

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

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

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

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

3 數(shù)值算法

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

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

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

式中，

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

4 算例分析

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

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

［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］李遠林.船舶非線性大幅橫搖導(dǎo)致傾覆的研究［J］.華南理工大學(xué)學(xué)報:自然科學(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ù)的識別［J］.天津大學(xué)學(xué)報，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］李遠林，伍曉榕.非線性橫搖阻尼的試驗確定-數(shù)據(jù)處理方法［J］.華南理工大學(xué)學(xué)報:自然科學(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］朱仁傳，郭海強，繆國平，等.一種基于CFD 理論船舶附加質(zhì)量與阻尼的計算方法［J］.上海交通大學(xué)學(xué)報，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］黃昊，郭海強，朱仁傳，等.粘性流中船舶橫搖阻尼計算［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)用非線性動力學(xué)［M］.北京:航空工業(yè)出版社，2000.

［12］陳樹輝.強非線性振動系統(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)單個波的波陡分布與破碎波［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］.中國海洋大學(xué)學(xué)報，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.