李焱，王麗元,3，封培元，劉震，劉崢，唐友剛

(1.天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300350；2.天津大學(xué) 天津大學(xué)建筑工程學(xué)院，天津 300350；3.江蘇科技大學(xué) 海洋裝備研究院，江蘇 鎮(zhèn)江 212003；4.上海市船舶工程重點(diǎn)實(shí)驗(yàn)室，上海 200011；5.中國船舶及海洋工程設(shè)計(jì)研究院，上海 200011)

船舶騎浪運(yùn)動(dòng)是船舶第二代完整穩(wěn)性中的5種失效模式之一[1]，其非線性動(dòng)力學(xué)特性最為顯著。船舶騎浪運(yùn)動(dòng)失穩(wěn)模式，是指船舶在特定的波浪中航行時(shí)，船舶推進(jìn)控制系統(tǒng)失去對(duì)船的控制，并使船舶隨著波浪前進(jìn)的現(xiàn)象。即船舶在隨浪條件下航行時(shí)會(huì)被前進(jìn)的波浪所捕獲，并與波浪一起運(yùn)動(dòng)，使船舶的航速與波速相同。通常使船舶發(fā)生騎浪運(yùn)動(dòng)的波長約為船長的1.0～3.0倍，波浪足夠陡峭，且船速與波速相當(dāng)，船速約為波速的75%[2]。

陶醉等[3-4]建立隨浪條件下船舶的單自由度縱蕩運(yùn)動(dòng)方程，通過將船舶騎浪運(yùn)動(dòng)與單擺運(yùn)動(dòng)同宿分岔現(xiàn)象進(jìn)行類比，研究船舶縱蕩非線性系統(tǒng)中所存在的同宿分岔現(xiàn)象，并預(yù)報(bào)船舶是否發(fā)生騎浪運(yùn)動(dòng)。儲(chǔ)紀(jì)龍等[5]建立船舶騎浪運(yùn)動(dòng)方程，將其簡化為哈密頓系統(tǒng)，求出其對(duì)應(yīng)的異宿軌道方程，計(jì)算船舶發(fā)生騎浪運(yùn)動(dòng)的螺旋槳臨界轉(zhuǎn)速。顧民等[6]以內(nèi)傾船型為例，開展隨浪和尾斜浪中的船舶騎浪/橫甩試驗(yàn)，分析了內(nèi)傾船型在不同波浪條件下的運(yùn)動(dòng)特性。

Kan等[7]通過模型試驗(yàn)和數(shù)值模擬結(jié)果相對(duì)比，給出一種獲得臨界航速和臨界波高的簡便方法，判斷船舶是否發(fā)生騎浪運(yùn)動(dòng)：當(dāng)船體的速度低于第一臨界速度時(shí)，船舶在縱蕩方向上作往復(fù)的周期運(yùn)動(dòng)，并不會(huì)發(fā)生騎浪運(yùn)動(dòng)；在第二臨界速度下船舶被波浪捕獲，與波浪之間的相對(duì)速度為零，發(fā)生騎浪運(yùn)動(dòng)。Spyrou等[8]研究在大波陡隨浪條件下船舶縱蕩運(yùn)動(dòng)的強(qiáng)非線性特性，研究表明在船舶騎浪條件的邊緣，船舶會(huì)發(fā)生大幅縱蕩運(yùn)動(dòng)。

本文基于非線性動(dòng)力學(xué)理論，研究船舶騎浪運(yùn)動(dòng)特性。建立縱浪條件下船舶縱蕩運(yùn)動(dòng)方程，分別采用梅林科夫方法和龍格庫塔方法對(duì)該方程進(jìn)行求解。通過解析分析方法，得到發(fā)生騎浪的臨界值，給出梅林科夫函數(shù)閾值圖。在數(shù)值計(jì)算中，分別繪制船舶騎浪時(shí)和未騎浪時(shí)的時(shí)歷曲線和運(yùn)動(dòng)相圖，并分析其非線性特性。以螺旋槳轉(zhuǎn)速為分岔參數(shù)畫出船舶騎浪運(yùn)動(dòng)分岔圖，分析了船舶騎浪的失穩(wěn)條件。

1 縱蕩運(yùn)動(dòng)方程建立

船舶在波浪上航行的坐標(biāo)系如圖1所示。O(x,y,z)為大地固定坐標(biāo)系，O′(x′,y′,z′)定義為隨船坐標(biāo)系，其中以船舶的重心G為隨船坐標(biāo)系坐標(biāo)原點(diǎn)，O″(x″,y″,z″)為波浪坐標(biāo)系，其原點(diǎn)固定于波浪的一個(gè)波谷處，并以波速c在大地坐標(biāo)系中移動(dòng)。

圖1 船舶運(yùn)動(dòng)坐標(biāo)系

船舶在規(guī)則波波浪激勵(lì)下縱蕩方向的經(jīng)典的運(yùn)動(dòng)方程為[9]：

(1)

式中：m為船舶的質(zhì)量；mx為縱蕩附連水質(zhì)量；ξG為大地坐標(biāo)系下船舶重心的縱坐標(biāo)；R(u)為船舶的縱向阻力，其表達(dá)式為：

R(u)=r1u+r2u2+r3u3

式中：r1、r2、r3為阻力系數(shù)；u為縱蕩瞬時(shí)速度。

T(u,n)為縱蕩方向螺旋槳推力，其表達(dá)式為：

式中：tp表示為螺旋槳推力減額系數(shù)；ρ為海水密度；n為螺旋槳轉(zhuǎn)速；Dp為螺旋槳槳葉直徑；KT表示為螺旋槳推力系數(shù)，其表達(dá)式為：

KT(u,n)=k1(u,n)+k2(u,n)+k3(u,n)

式中：k1、k2、k3為螺旋槳性能參數(shù)，代表螺旋槳敞水特性曲線的回歸系數(shù)；JP為螺旋槳進(jìn)速系數(shù)，表達(dá)式為：

JP=u(1-ωp)/(nDP)

式中ωp為螺旋槳性能參數(shù)中的伴流系數(shù)。

Xw(ξG,t)為船舶在縱蕩方向上受到的波浪力，按切片理論計(jì)算[10]：

(2)

式中；g表示重力加速度；A代表波浪幅值；k為波數(shù)；χ為船舶航向角；S(x)為船舶橫剖面浸水面積；d(x)為各站對(duì)應(yīng)吃水；B(x)為各站船寬；dx為站距;α為考慮繞射效應(yīng)的波浪力修正系數(shù)，按照下式進(jìn)行計(jì)算：

式中：Cb為船體方形系數(shù)；Cm為船體中橫剖面系數(shù)。

考慮規(guī)則波條件計(jì)算船體受到的波浪力，波浪力隨時(shí)間正弦變化。由于各時(shí)刻波浪力大小與船-波相對(duì)位置有關(guān)，因此將波浪力寫作波浪力幅值乘以船舶與波浪之間相對(duì)位置的正弦函數(shù)的形式：

(3)

2 梅林科夫方法

梅林科夫(Melnikov)方法始于20世紀(jì)60年代，通過度量龐加萊(Poincare)映射的雙曲不動(dòng)點(diǎn)的穩(wěn)定與不穩(wěn)定流形之間的距離來確定非線性動(dòng)力系統(tǒng)是否存在橫截同(異)宿點(diǎn)，從而判斷動(dòng)力系統(tǒng)是否發(fā)生運(yùn)動(dòng)的突變，在確定性系統(tǒng)和隨機(jī)系統(tǒng)中都有大量應(yīng)用[11-12]。

首先假設(shè)動(dòng)力系統(tǒng)的參數(shù)：

x=f(x)+εg(x,t;μ)

(4)

其中：f:R2→R2是系統(tǒng)對(duì)應(yīng)的哈密頓矢量場，εg:R2×R×R→R2代表作用在系統(tǒng)上的一個(gè)小擾動(dòng)，μ是系統(tǒng)參數(shù)。此后計(jì)算系統(tǒng)的梅林可夫積分函數(shù)：

(5)

式中xh(t)表示動(dòng)力系統(tǒng)的異宿軌道方程。

在梅林可夫函數(shù)表達(dá)式中隱含時(shí)間t，可通過求解梅林可夫函數(shù)的零點(diǎn)來找到船舶縱蕩運(yùn)動(dòng)系統(tǒng)失穩(wěn)時(shí)的參數(shù)閾值，從而得到船舶發(fā)生騎浪失穩(wěn)的條件。

(6)

T(c;n)-R(c)

(7)

其中：

(8)

其中：

引入小參數(shù)ε，其中0<ε?1，則式(8)轉(zhuǎn)化為：

(9)

當(dāng)小參數(shù)ε=0時(shí)，式(9)退化為沒有外激勵(lì)作用，沒有阻尼的哈密頓系統(tǒng)：

(10)

(11)

將哈密頓系統(tǒng)式(11)中的右端項(xiàng)取不定積分，并由第1項(xiàng)減去第2項(xiàng)，則得到哈密頓系統(tǒng)(11)對(duì)應(yīng)的Hamilton量H(y,z)：

(12)

當(dāng)Hamilton量H(±π,0)=1時(shí)，整個(gè)哈密頓系統(tǒng)存在著連接(π,0)和(-π,0)這2個(gè)鞍點(diǎn)的2條異宿軌道，在異宿軌道之間，存在著幾條圍繞著(0,0)中心的同宿軌道，從同宿軌道變換成的第1條的異宿軌道則表示此縱蕩運(yùn)動(dòng)系統(tǒng)中的船舶是否發(fā)生騎浪失穩(wěn)運(yùn)動(dòng)的分界線。其中系統(tǒng)的異宿軌道方程為：

z=±2cos(y/2)

(13)

圖2中，虛線為本系統(tǒng)的異宿軌道線，虛線內(nèi)部的實(shí)線則表示系統(tǒng)的同宿軌道方程。

圖2 異宿軌道示意

省略式(9)右端項(xiàng)中的小參數(shù)，并將其沿著異宿軌道方程進(jìn)行積分，得到縱蕩運(yùn)動(dòng)方程的梅林可夫函數(shù)為：

(14)

(15)

改變函數(shù)的積分路線，使梅林可夫函數(shù)沿著式(13)所示的異宿軌道進(jìn)行積分，并使積分后的結(jié)果M=0，則：

(16)

(17)

對(duì)式(17)進(jìn)行積分，得到I1=4，I2=π，I3=8/3，代入式(16)得：

(18)

由此，根據(jù)式(18)可求出使等式成立的臨界螺旋槳轉(zhuǎn)速n。再根據(jù)船舶在波浪中前行時(shí)的阻力和推力相等的關(guān)系式，求出船舶在波浪中行駛時(shí)的螺旋槳臨界船速，便確定了船舶騎浪失穩(wěn)的條件。

3 某內(nèi)傾船算例

以某內(nèi)傾船為例，分別采用解析方法和數(shù)值方法計(jì)算內(nèi)傾船的縱蕩運(yùn)動(dòng)響應(yīng)，并判斷其是否發(fā)生騎浪運(yùn)動(dòng)，給出內(nèi)傾船發(fā)生騎浪運(yùn)動(dòng)的臨界值條件，并通過2種方法結(jié)果的比較，驗(yàn)證解析方法的正確性。內(nèi)傾船垂線間長為112 m，吃水為3.865 m，排水體積為2 789.72 m3，水線寬為13.399 m，水線長為111.997 m，濕表面積為1 518.31 m2，浮心縱向位置距舯后-2.456 m，重心垂向位置距基線5.785 m，中橫剖面系數(shù)為0.794 6。內(nèi)傾船的縱蕩運(yùn)動(dòng)方程中的推力系數(shù)和阻力系數(shù)由模型試驗(yàn)的數(shù)據(jù)給出。

3.1 臨界轉(zhuǎn)速與臨界船速計(jì)算結(jié)果

選取波浪條件為λ/Lpp=1.5，航向角χ=0°?？紤]波幅為漸變量，采用解析方法求解船舶發(fā)生騎浪運(yùn)動(dòng)時(shí)螺旋槳臨界轉(zhuǎn)速與波幅之間的關(guān)系，采用梅林科夫方法計(jì)算得到的臨界螺旋槳轉(zhuǎn)速閾值曲線如圖3所示。由圖3可以看出當(dāng)波長保持不變、波幅增加時(shí)，梅林科夫方法得出的螺旋槳臨界轉(zhuǎn)速逐漸減低。即在波幅較大時(shí)，螺旋槳臨界轉(zhuǎn)速較小，船舶易于發(fā)生騎浪現(xiàn)象。相反，當(dāng)波幅較小時(shí)，螺旋槳臨界轉(zhuǎn)速則更高，騎浪現(xiàn)象更難發(fā)生。

圖3 螺旋槳臨界轉(zhuǎn)速隨波幅變化曲線

根據(jù)船舶的推力和阻力相等條件，可以進(jìn)一步求出發(fā)生騎浪的臨界船速。船舶臨界船速隨波幅變化的關(guān)系圖如圖4所示。由圖4可以看出，船舶臨界船速隨波幅變化趨勢與船舶螺旋槳臨界轉(zhuǎn)速隨波幅的變化趨勢相同。即當(dāng)波幅值增加時(shí)，梅林科夫方法得出船舶臨界航速逐漸減低。這一結(jié)果說明，在波幅較大時(shí)，船舶臨界航速較小，船舶就會(huì)易于發(fā)生騎浪現(xiàn)象。相反，當(dāng)波幅較小時(shí)，船舶臨界航速則更高，不易發(fā)生騎浪。

圖4 船舶臨界船速隨波幅變化曲線

3.2 數(shù)值模擬計(jì)算結(jié)果

采用龍格庫塔方法對(duì)船舶縱蕩運(yùn)動(dòng)進(jìn)行數(shù)值求解，計(jì)算不同螺旋槳轉(zhuǎn)速下船舶縱蕩運(yùn)動(dòng)，通過對(duì)比縱蕩運(yùn)動(dòng)的數(shù)值解形式，判斷螺旋槳的臨界轉(zhuǎn)速。

選取波浪條件為：λ/Lpp=1.5，航向角χ=0°。H/λ=1/10。在此波浪條件下，由梅林科夫方法得到的螺旋槳臨界轉(zhuǎn)速為3.38 r/s，對(duì)應(yīng)臨界航速為10.38 m/s，采用龍格庫塔方法得出的螺旋槳臨界轉(zhuǎn)速為3.5 r/s，對(duì)應(yīng)的臨界船速為10.69 m/s。將2種方法得出的結(jié)果進(jìn)行對(duì)比，得出2種算法預(yù)報(bào)的螺旋槳臨界轉(zhuǎn)速誤差為3.43%，臨界船速誤差為2.90%，均在5%以內(nèi)，驗(yàn)證了解析方法的準(zhǔn)確性。

當(dāng)螺旋槳轉(zhuǎn)速為3.3 r/s時(shí)，船舶與波浪間的相對(duì)速度如圖5(a)所示，與波浪的相對(duì)位移圖5(b)所示。從2個(gè)時(shí)間歷程曲線圖中可以看出，當(dāng)螺旋槳轉(zhuǎn)速為3.3 r/s時(shí)，船舶的縱蕩速度相對(duì)于波浪速度是周期變化的，同時(shí)船舶與波浪的相對(duì)位移隨時(shí)間是逐步增大的。此時(shí)船由螺旋槳進(jìn)行推進(jìn)，并未被波浪所捕獲，所以在螺旋槳轉(zhuǎn)速為3.3 r/s時(shí)，船舶未發(fā)生騎浪現(xiàn)象。

圖5 轉(zhuǎn)速3.3 r/s時(shí)船舶縱蕩運(yùn)動(dòng)響應(yīng)

當(dāng)螺旋槳轉(zhuǎn)速增加到3.5 r/s時(shí)，船與波的相對(duì)速度及相對(duì)位移如圖6(a)和圖6(b)所示。由2個(gè)時(shí)間歷程曲線圖中可以看出，當(dāng)螺旋槳轉(zhuǎn)速為3.5 r/s時(shí)，船舶在波浪中運(yùn)動(dòng)約200 s后，船舶與波浪之間的相對(duì)速度變?yōu)榱悖瑫r(shí)船舶與波浪之間的相對(duì)位移也保持不變。從這2個(gè)時(shí)歷曲線圖中可以看出，當(dāng)螺旋槳轉(zhuǎn)速為3.5 r/s時(shí)，船舶在運(yùn)動(dòng)一段時(shí)間后會(huì)以波浪速度前進(jìn)，即船舶被波浪所捕獲。此時(shí)船的螺旋槳推力系統(tǒng)失去對(duì)船舶的控制，發(fā)生了船舶騎浪運(yùn)動(dòng)現(xiàn)象。

圖6 轉(zhuǎn)速3.3 r/s時(shí)船舶縱蕩運(yùn)動(dòng)響應(yīng)

圖7 船舶運(yùn)動(dòng)相圖

3.3 船舶騎浪運(yùn)動(dòng)分岔特性

取波浪條件為λ/Lpp=1.5，航向角χ=0°，H/λ=1/10。給定不同的螺旋槳轉(zhuǎn)速n，求船舶縱蕩運(yùn)動(dòng)的時(shí)間歷程。船舶相對(duì)速度隨螺旋槳轉(zhuǎn)速變化的分岔圖如圖8所示。從圖8中看出，此波浪條件下，當(dāng)螺旋槳轉(zhuǎn)速低于3.4 r/s時(shí)，系統(tǒng)發(fā)生的是穩(wěn)定的周期運(yùn)動(dòng)，在分岔圖的表現(xiàn)上表示為一個(gè)點(diǎn)，即一個(gè)圓圈。當(dāng)螺旋槳轉(zhuǎn)速為3.5 r/s時(shí)。分岔圖上不能畫出點(diǎn)，此處用點(diǎn)集區(qū)域示意船舶騎浪的非穩(wěn)定區(qū)域的范圍，不再是穩(wěn)定的往復(fù)周期性運(yùn)動(dòng)，此時(shí)船舶與波浪之間的相對(duì)速度為0。所以此條件下船舶縱蕩運(yùn)動(dòng)喪失穩(wěn)定性，即船舶發(fā)生騎浪失穩(wěn)運(yùn)動(dòng)。

圖8 船舶相對(duì)速度隨螺旋槳轉(zhuǎn)速變化分岔圖

4 結(jié)論

1)基于非線性動(dòng)力學(xué)方法得到的騎浪條件和騎浪失穩(wěn)域，與數(shù)值方法得到的結(jié)論基本吻合，表明非線性動(dòng)力學(xué)方法分析船舶騎浪運(yùn)動(dòng)的可行性。

2)當(dāng)船舶在規(guī)則波中未發(fā)生騎浪運(yùn)動(dòng)時(shí)，船速相對(duì)于波速呈周期性變化的規(guī)律；當(dāng)船舶發(fā)生騎浪運(yùn)動(dòng)時(shí)，船舶被波浪捕獲，相對(duì)速度為零。

3)當(dāng)波幅較大，螺旋槳臨界轉(zhuǎn)速較小時(shí)，船舶容易發(fā)生騎浪運(yùn)動(dòng)，并且易發(fā)生騎浪失穩(wěn)導(dǎo)致橫甩。