唐友剛 王麗元 劉 崢 劉利琴 李 焱

（1. 天津大學(xué) 建筑工程學(xué)院 天津 300350； 2. 江蘇科技大學(xué) 海洋裝備研究院 鎮(zhèn)江 212003）

0 引 言

針對(duì)船舶二代穩(wěn)性中的失效模式，目前主要的研究方法和手段有非線性動(dòng)力學(xué)方法、模型試驗(yàn)方法和數(shù)值計(jì)算分析方法等。研究內(nèi)容主要包括各種失效模式的發(fā)生機(jī)理以及如何避免船舶發(fā)生失穩(wěn)。非線性動(dòng)力學(xué)在研究船舶二代穩(wěn)性中，有著不可替代的作用，是研究船舶失穩(wěn)機(jī)理的主要手段。本文基于作者所在團(tuán)隊(duì)多年科研工作，總結(jié)船舶非線性動(dòng)力學(xué)分析方法在船舶參數(shù)橫搖、騎浪橫甩方面的工程應(yīng)用，探討船舶非線性運(yùn)動(dòng)響應(yīng)和穩(wěn)性的研究方法，揭示船舶大幅運(yùn)動(dòng)響應(yīng)、運(yùn)動(dòng)失穩(wěn)直至傾覆全過程的非線性動(dòng)力特性。

1 船舶參數(shù)橫搖非線性分析

1.1 規(guī)則波中非線性響應(yīng)

規(guī)則縱浪中船舶參數(shù)激勵(lì)非線性橫搖方程為：

采用多尺度法求解上述方程，得到近似解析解。當(dāng)φ = 0時(shí)，得到方程的平凡解；但當(dāng)參數(shù)激勵(lì)的幅值、頻率滿足一定的條件時(shí)，將得到方程的非平凡解，橫搖就會(huì)被激起。在參數(shù)空間內(nèi)確定平凡解的穩(wěn)定域與非穩(wěn)定域，也就確定了參數(shù)激勵(lì)橫搖發(fā)生的條件。

由圖1可見，當(dāng)=0.4時(shí)，若遭遇頻率與橫搖固有頻率的比值Ω從0開始增加，對(duì)于主參數(shù)共振來講，漁政船的穩(wěn)定橫搖幅值將在T處發(fā)生跳躍，出現(xiàn)大幅橫搖。繼續(xù)增加Ω，橫搖幅值逐漸變小，在T之后周期解失穩(wěn)，平凡解穩(wěn)定。

圖1 漁政船幅頻曲線h0=0.4

由圖2可見，當(dāng)=0.9時(shí)，若Ω從1.5開始增加，則船舶將在T處傾覆，平凡解失穩(wěn)，而沒有穩(wěn)定的周期解。

圖2 漁政船幅頻曲線h0=0.9

1.2 非規(guī)則波中非線性響應(yīng)

船舶在隨機(jī)斜浪中航行，假設(shè)船舶的縱向運(yùn)動(dòng)為小量，得到船舶參-強(qiáng)激勵(lì)單自由度橫搖運(yùn)動(dòng)方程如下：

式中：為橫搖角，rad；I為橫搖慣性矩，kg·m；A()為橫搖附加慣性矩，kg · m；為線性阻尼系數(shù)，N· m· s；為立方阻尼系數(shù)，N· m· s； Δ為排水量，kg；為重力加速度，m/s；為船舶復(fù)原力臂，m，由、和決定；()為波浪力矩，N·m。

采用能量包線隨機(jī)平均法研究船舶隨機(jī)參-強(qiáng)激勵(lì)橫搖運(yùn)動(dòng)概率密度函數(shù)。將船舶運(yùn)動(dòng)過程處理為馬爾可夫過程，將運(yùn)動(dòng)響應(yīng)分成快變量與慢變量，通過對(duì)快變量的平均得到慢變量的近似方程。

將隨機(jī)波面升高處理為窄帶隨機(jī)過程，即用有效波面Z (，)來代替隨機(jī)波面(，)。將Z (，)展開為2個(gè)隨機(jī)過程η()和η()的組合，具體如下：

式中：η()和η()均為高斯平穩(wěn)隨機(jī)過程，經(jīng)證明兩者相互獨(dú)立。

作尺度變換如下：

用替代t ，則運(yùn)動(dòng)微分方程為：

將響應(yīng)分量分為快變量和慢變量，即系統(tǒng)的狀態(tài)空間由（，）變換為（，），其中為系統(tǒng)總能量：

式中：每個(gè)能量值（，）對(duì)應(yīng)于特定的橫搖運(yùn)動(dòng)。例如：當(dāng)=23.73、=19.74 時(shí)，得到對(duì)應(yīng)的能量等值線，如圖3所示。

圖3 能量等值線

運(yùn)動(dòng)方程進(jìn)一步寫為：

引入中間變量（，），令：

則（11）可進(jìn)一步寫為如下隨機(jī)微分方程式：

根據(jù)隨機(jī)平均法，其漂移與擴(kuò)散系數(shù)分別為：

式（13）與式（14）中，（）為時(shí)間區(qū)間：

針對(duì)0≤≤H，采用雅可比橢圓函數(shù)求解運(yùn)動(dòng)微分方程，得到：

式中：、、為雅克比橢圓函數(shù)。符號(hào)定義為：

將式（15）至（18）代入式（13）與式（14），并考慮上述的符號(hào)定義，得到漂移與擴(kuò)散系數(shù)為：

系統(tǒng)能量滿足如下伊藤隨機(jī)微分方程：

式中：（，）表示初始狀態(tài)，該式可表達(dá)為式（23）所示形式：

引入尺度函數(shù)()與速度函數(shù)()：

式（25）可進(jìn)一步寫為：

如果系統(tǒng)能量＞H，方程（28）不存在穩(wěn)定解；如果≤H，則式（28）的平穩(wěn)概率密度函數(shù)為：

式中：參數(shù)、由邊界條件確定。

對(duì)于低強(qiáng)度噪聲激勵(lì)有=0，則式（29）可進(jìn)一步寫為：

應(yīng)用傳遞函數(shù)關(guān)系p()=p()(d/d)，得到基于變量的平穩(wěn)概率密度函數(shù)為：

對(duì)式（30）在某一角度范圍[，]內(nèi)積分，可進(jìn)一步得到()在該范圍內(nèi)的概率。

應(yīng)用蒙特卡羅方法求解C11型集裝箱船在實(shí)尺度下的橫搖響應(yīng)概率密度函數(shù)。隨機(jī)種子數(shù)=25、迭代步數(shù)=5 000 000、時(shí)間步長Δ=0.01 s，得到的概率密度函數(shù)的數(shù)值解，其與隨機(jī)平均法得到的結(jié)果對(duì)比如下頁圖4所示。對(duì)比結(jié)果表明，兩者吻合較好，同時(shí)也證實(shí)本文采用的解析解法和數(shù)值解法的正確性。

圖4 概率密度函數(shù)對(duì)比

1.3 波群分析方法

波浪頻譜不能反映波群的特性，本文用波包來描述波群，用波包譜作為靶譜來模擬相位譜，波包譜S()的靶譜：

同時(shí)，當(dāng)≤0.7，f / f=5~15，當(dāng)＞0.7，f / f=10~28。

采用波包譜數(shù)值模擬不規(guī)則波的波群過程如圖5所示。

圖5 波群模擬流程圖

JONSWAP公式如下所示：

GFH和GLF分別表征波群的高度和長度特征， 參見圖6。對(duì)于給定的頻譜，當(dāng)GFH不變，隨著 GLF增大，波包靶譜面積基本不變，峰頻減小，波包變平坦，群長增大；當(dāng)GLF不變，隨著GFH增大，峰頻不變，波包譜面積增大，波包線的波動(dòng)幅度增大，群高有明顯變化。群高和群長對(duì)參數(shù)橫搖有明顯影響，群長的影響更顯著。

圖6 波群波面升高及船舶橫搖角（GFH=0.8，GLF=9.82）

2 船舶騎浪運(yùn)動(dòng)非線性分析

參照?qǐng)D7坐標(biāo)系，船舶在波浪中縱蕩運(yùn)動(dòng)方程為：

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

式中：為船舶的質(zhì)量，kg；m為縱蕩附加質(zhì)量，kg；ξ為船舶重心的縱坐標(biāo)，m；()為船舶阻力，N；為縱向瞬時(shí)速度，m/s；(，)為船舶螺旋 槳推力，N；為螺旋槳轉(zhuǎn)速，r/s；F(ξ，)為船 舶在縱向受到的波浪力，N。

船舶靜水阻力()為：

式中：、、為多項(xiàng)式的系數(shù)。

船舶靜水中螺旋槳推力（，）為：

式中：t表示為推力減額分?jǐn)?shù)；為海水密度，kg/m；為轉(zhuǎn)速，r/s；D為直徑，m；K為推力系數(shù)。

基于切片理論，船舶在波浪中運(yùn)動(dòng)時(shí)受到的波浪力可表達(dá)為：

式中：為繞射效應(yīng)修正系數(shù)；C為方形系數(shù)；C為中橫剖面系數(shù)；為海水密度，kg/m；為重力加速度，m/s；為波浪幅值，m；為波浪的波數(shù)，m；為船舶航向角，rad；為船舶每站橫剖面的濕面積，m；為每站的吃水，m；d為每站間隔距離，m；為每站的船體寬度，m。

式中：c可用多項(xiàng)式表達(dá)如下：

在方程（48）中引入小參數(shù)，其中0＜1，將縱蕩方程（48）轉(zhuǎn)為：

2.1 規(guī)則波中騎浪運(yùn)動(dòng)分岔研究

當(dāng)小參數(shù)0時(shí)，式（50）退化為無阻尼和外激勵(lì)作用的哈密頓系統(tǒng)：

將式（52）中的右端項(xiàng)取不定積分，并由第1項(xiàng)減去第2項(xiàng)，則得到方程（52）對(duì)應(yīng)哈密頓的量（，）：

當(dāng)哈密頓的量（+π，0）=1時(shí)，系統(tǒng)軌道如圖8所示。

圖8 異宿軌道示意圖

其中，紅色線為系統(tǒng)的異宿軌道線，紅色線內(nèi)部的藍(lán)線則是表示系統(tǒng)的同宿軌道線。整個(gè)哈密頓系統(tǒng)存在著連接（π，0）和（-π，0）2個(gè)鞍點(diǎn)的 2條異宿軌道，在異宿軌道之間，存在著幾條圍繞著 （0，0）中心的同宿軌道。同宿軌道與異宿軌道的分界線，即為船舶縱蕩運(yùn)動(dòng)中發(fā)生騎浪現(xiàn)象的分界線。系統(tǒng)的異宿軌道方程為：

忽略式（50）右端項(xiàng)中的小參數(shù)，將其沿異宿軌道積分，得到船舶縱蕩運(yùn)動(dòng)方程的梅林可夫函數(shù)：

沿著（54）所示的異宿軌道進(jìn)行積分，并使=0，則有：

式（59）當(dāng)中除了參數(shù)未知，所有其他參數(shù)都是已知的，所以根據(jù)式（59）可求出臨界螺旋槳轉(zhuǎn)速；再根據(jù)船舶在波浪中前行時(shí)的阻力和推力相等，求出船舶的臨界船速，確定船舶騎浪失穩(wěn)的條件。

以某內(nèi)傾船為例開展研究，其主尺度如表1所示，推力系數(shù)和阻力系數(shù)由模型試驗(yàn)得到。分別采用解析方法和數(shù)值方法判斷該船是否發(fā)生騎浪運(yùn)動(dòng)，給出對(duì)應(yīng)的臨界值條件。通過2種方法的比較，驗(yàn)證解析方法的正確性。

表1 內(nèi)傾船主尺度

當(dāng)波浪船長比/L=1.5，航向角χ=0°時(shí)，船舶發(fā)生騎浪運(yùn)動(dòng)的螺旋槳臨界轉(zhuǎn)速與波幅之間的關(guān)系如圖9所示。由圖可以看出，隨著波幅值增加，螺旋槳臨界轉(zhuǎn)速減低。

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

船舶臨界船速隨波幅的變化如圖10所示。

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

當(dāng)波幅值增加時(shí)，船舶臨界航速減低。波幅較大時(shí)，臨界航速較小，船舶更容易發(fā)生騎浪現(xiàn)象。

當(dāng)波長船長比λ/L=1.5，航向角=0°時(shí)，取1/10，基于解析方法和數(shù)值方法分析規(guī)則波中的船舶縱蕩運(yùn)動(dòng)。由解析方法得到的螺旋槳臨界轉(zhuǎn)速為3.38 r/s，對(duì)應(yīng)臨界航速為10.38 m/s。由數(shù)值方法得到的螺旋槳臨界轉(zhuǎn)速為3.5 r/s，對(duì)應(yīng)的臨界船速為10.69 m/s。2種方法得出的結(jié)果相差很小，誤差在5%左右，從而驗(yàn)證了2種方法的正確性。

當(dāng)螺旋槳轉(zhuǎn)速為3.4 r/s時(shí)，船舶相對(duì)于波浪的相對(duì)速度和位移分別如圖11和圖12所示。由圖中可以看出，船舶與波浪的相對(duì)速度是周期變化的，相對(duì)位移隨時(shí)間增大，即船舶尚未被波浪所捕獲，未發(fā)生騎浪現(xiàn)象。

圖11 轉(zhuǎn)速3.4 r/s時(shí)相對(duì)速度時(shí)間歷程曲線

圖12 轉(zhuǎn)速3.4 r/s時(shí)相對(duì)位移時(shí)間歷程曲線

當(dāng)螺旋槳轉(zhuǎn)速增加到3.5 r/s時(shí)，船與波的相對(duì)速度和相對(duì)位移分別如圖13和圖14所示。由圖中可以看出，約200 s后，船舶與波浪之間的相對(duì)速度為0，船舶與波浪之間的相對(duì)位移也保持為一個(gè)常數(shù)。表明船舶在運(yùn)動(dòng)一段時(shí)間后會(huì)被波浪捕獲，以波浪速度前進(jìn)。此時(shí)，螺舵系統(tǒng)失去對(duì)船舶的控制，發(fā)生了船舶騎浪運(yùn)動(dòng)現(xiàn)象。

圖13 轉(zhuǎn)速3.5 r/s時(shí)相對(duì)速度時(shí)間歷程曲線

圖14 轉(zhuǎn)速3.5 r/s時(shí)相對(duì)位移時(shí)間歷程曲線

為了更加準(zhǔn)確地判斷船舶騎浪運(yùn)動(dòng)的發(fā)生，以cos(2π·ξ″/)為相圖橫坐標(biāo)，船的絕對(duì)速度為相圖縱坐標(biāo)，分別給出螺旋槳轉(zhuǎn)速為3.4 r/s和3.5 r/s時(shí)船舶運(yùn)動(dòng)的相圖。由下頁圖15和圖16可以看出：當(dāng)螺旋槳轉(zhuǎn)速為3.4 r/s時(shí)，船舶縱蕩運(yùn)動(dòng)系統(tǒng)的相圖是封閉的，此時(shí)船舶運(yùn)動(dòng)是穩(wěn)定的，未發(fā)生騎浪現(xiàn)象；當(dāng)螺旋槳轉(zhuǎn)速為3.5 r/s時(shí)，船舶縱蕩運(yùn)動(dòng)系統(tǒng)的相圖不是封閉的，此時(shí)船舶運(yùn)動(dòng)喪失穩(wěn)定性，發(fā)生了船舶騎浪現(xiàn)象。

圖15 3.4 r/s船舶運(yùn)動(dòng)相圖（未騎浪）

圖16 3.5 r/s船舶運(yùn)動(dòng)相圖（騎浪發(fā)生）

取波長船長比/L=1.5，航向角=0°，/= 1/10。針對(duì)不同的螺旋槳轉(zhuǎn)速，數(shù)值求解船舶縱蕩運(yùn)動(dòng)，得到運(yùn)動(dòng)時(shí)歷曲線。以螺旋槳轉(zhuǎn)速為橫坐標(biāo)，船舶在波浪中的相對(duì)速度運(yùn)動(dòng)的穩(wěn)定響應(yīng)幅值為縱坐標(biāo)，給出船舶相對(duì)速度隨螺旋槳轉(zhuǎn)速變化的分岔圖，如圖17所示。由圖可見：當(dāng)螺旋槳轉(zhuǎn)速低于3.4 r/s時(shí)，船舶縱蕩運(yùn)動(dòng)是穩(wěn)定的周期運(yùn)動(dòng)，在分岔圖的表現(xiàn)上表示為一個(gè)點(diǎn)；當(dāng)螺旋槳轉(zhuǎn)速為3.5 r/s時(shí)，分岔圖上不能畫出點(diǎn)，此處用藍(lán)色的圓點(diǎn)示意非穩(wěn)定區(qū)域，此時(shí)船舶縱蕩運(yùn)動(dòng)開始發(fā)散，即船舶發(fā)生騎浪失穩(wěn)運(yùn)動(dòng)。

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

2.2 非規(guī)則波中騎浪概率研究

建立船舶縱蕩運(yùn)動(dòng)微分方程，采用隨機(jī)Melnikov非線性動(dòng)力學(xué)方法，求解船舶在隨機(jī)波浪中發(fā)生騎浪的概率。

通過線性疊加求解隨機(jī)波浪力：

是一個(gè)隨機(jī)過程，其譜密度函數(shù)為S。假定不規(guī)則波浪力為在規(guī)則波的基礎(chǔ)上疊加一隨機(jī)的波浪力，規(guī)則波周期和波高為對(duì)應(yīng)波浪譜的譜峰周期和有義波高?；谝陨霞僭O(shè)，可以表示為如下 形式：

式中：ω為波浪譜的譜峰頻率，rad/s；()為隨機(jī)相位角，rad。

采用Melnikov函數(shù)來表征穩(wěn)定流形與不穩(wěn)定流形之間的距離，如下表示：

船舶在隨機(jī)波浪中發(fā)生騎浪的必要條件為()<0。隨機(jī)波浪產(chǎn)生的激勵(lì)可近似為高斯隨機(jī)過程，因此通過線性變化得到的Melnikov過程()也是高斯過程，均值μ和方差σ的表達(dá)式為：

式中：()為頻率響應(yīng)函數(shù)：

正態(tài)分布的概率密度函數(shù)為：

通過積分可得到該正態(tài)分布的概率分布函數(shù)：

根據(jù)上述概率分布函數(shù)，隨機(jī)波中騎浪發(fā)生的概率為：

分別研究波高和轉(zhuǎn)速對(duì)船舶騎浪發(fā)生概率的影響，計(jì)算結(jié)果如圖18和圖19所示。隨著有義波高的增加，船舶所受到的波浪載荷增加，騎浪概率也隨之增加。隨著螺旋槳轉(zhuǎn)速的增加，船舶的航行速度增加，當(dāng)航速越接近波速時(shí)，騎浪發(fā)生概率 變大。

圖18 不同波高下的騎浪概率

圖19 不同轉(zhuǎn)速下的騎浪概率

3 船舶橫甩非線性分析方法

基于累積艏搖運(yùn)動(dòng)原理，建立船舶艏搖運(yùn)動(dòng)方程如下：

將船舶的艏向角以艏搖偏差角和目標(biāo)艏向角的方式來表示，引入艏搖偏差角，=-Ψ，將其代入公式（70）可得：

引入時(shí)間尺度=，將其代入式（71）得到如下表達(dá)式：

引入小參數(shù)，并在計(jì)算完成后令小參數(shù)=1，式（72）變換為：

為研究艏搖運(yùn)動(dòng)的亞諧共振，令船舶遭遇頻率與固有頻率的比值滿足Ω≈2，并引入頻率比調(diào)諧因子來描述頻率比與2之間的差異，令：

將式（74）代入式（73），整理可得：

假設(shè)方程（75）的解為：

根據(jù)多尺度法，可得到方程的一階近似解為：

式中：變量和變量由如下兩式確定：

根據(jù)消除永年項(xiàng)條件，可得：

將公式（80）的近似解表達(dá)為：

式中：B、B都表示為實(shí)數(shù)，將表達(dá)式（81）代入到表達(dá)式（80）中，將方程的實(shí)部和虛部分離開，可得：

設(shè)方程（82）的解分別為：

式中：b和b都為常數(shù)，將公式（83）代入到公式（82）中可得到：

由此得到艏搖運(yùn)動(dòng)的穩(wěn)定域和不穩(wěn)定域如圖20所示，在頻率比為2時(shí)，船舶艏搖運(yùn)動(dòng)方程中的參數(shù)激勵(lì)項(xiàng)在達(dá)到0.62左右時(shí)就會(huì)發(fā)生失穩(wěn)。

圖20 解的穩(wěn)定性區(qū)域圖

以波高為分岔參數(shù)，得到船舶艏搖運(yùn)動(dòng)響應(yīng)幅值隨波高變化的分岔圖，見下頁圖21。在 2倍頻波浪條件下，當(dāng)波高<0.7 m時(shí)，系統(tǒng)為穩(wěn)定的周期運(yùn)動(dòng)，在分岔圖的表現(xiàn)上表示為1個(gè)紅點(diǎn)。當(dāng)波高＞0.7 m時(shí)，發(fā)生船舶發(fā)生艏搖失穩(wěn)運(yùn)動(dòng)。

圖21 2倍頻下船舶艏搖角隨波高變化分岔圖

由下頁圖23和圖24可見，當(dāng)波高取為1 m時(shí)，船舶在經(jīng)過幾個(gè)周期后，艏搖運(yùn)動(dòng)幅值越來越大。此時(shí)參數(shù)激勵(lì)的作用明顯增加，引起艏搖運(yùn)動(dòng)發(fā)生亞諧共振，此時(shí)相圖不再是1個(gè)閉合的圓圈，而是螺旋線，即船舶艏搖運(yùn)動(dòng)不斷增大導(dǎo)致喪失穩(wěn)定性，發(fā)生船舶橫甩現(xiàn)象。

圖24 波高1.0 m艏搖角時(shí)間歷程曲線圖

由圖22和圖23可見，當(dāng)波高為0.5 m時(shí)，參數(shù)激勵(lì)較小，不能激起系統(tǒng)的亞諧共振，船舶艏搖運(yùn)動(dòng)的相圖構(gòu)成1個(gè)閉合的圓圈，船舶艏搖運(yùn)動(dòng)是穩(wěn)定的，不會(huì)發(fā)生艏搖失穩(wěn)現(xiàn)象。

圖22 波高0.5 m艏搖角時(shí)間歷程曲線圖

圖23 波高0.5 m船舶艏搖運(yùn)動(dòng)相圖

圖25 波高1.0 m船舶艏搖運(yùn)動(dòng)相圖

4 結(jié) 論

本文回顧總結(jié)了基于非線性動(dòng)力學(xué)理論和方法開展船舶在波浪中的運(yùn)動(dòng)、傾覆過程和運(yùn)動(dòng)失穩(wěn)機(jī)理研究的主要進(jìn)展。主要結(jié)論如下：

（1）船舶運(yùn)動(dòng)（尤其大幅運(yùn)動(dòng)）是非線性動(dòng)力系統(tǒng)、非線性動(dòng)力學(xué)方法與耐波性理論的結(jié)合，是船舶穩(wěn)性研究的重要發(fā)展方向；

（2）非線性動(dòng)力學(xué)方法可以更為精準(zhǔn)和方便地描述船舶的瞬時(shí)運(yùn)動(dòng)狀態(tài)，比如采用相圖可以給出船舶瞬時(shí)的運(yùn)動(dòng)速度、幅值和趨勢(shì)等；

（3）采用哈爾米頓系統(tǒng)分析軌道的同宿與異宿，可以確定船舶不同運(yùn)動(dòng)形式（例如騎浪、橫甩）的穩(wěn)定與失穩(wěn)域的范圍；

（4）船舶失穩(wěn)導(dǎo)致運(yùn)動(dòng)過程出現(xiàn)不同狀態(tài)，采用分叉方法研究船舶運(yùn)動(dòng)失穩(wěn)的不同狀態(tài)（包括傾覆）是可行的方法；

（5）今后研究的方向是將耐波性理論、波浪理論及非線性動(dòng)力學(xué)理論相結(jié)合，探索船舶運(yùn)動(dòng)和穩(wěn)性研究的新概念和新方法，開辟船舶運(yùn)動(dòng)和穩(wěn)性研究的新途徑和技術(shù)手段，使船舶運(yùn)動(dòng)安全的研究和表達(dá)更為精準(zhǔn)。