船舶非線性動力學分析方法及工程應用

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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

船舶在隨機斜浪中航行，假設船舶的縱向運動為小量，得到船舶參-強激勵單自由度橫搖運動方程如下：

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

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

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

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

作尺度變換如下：

用替代t ，則運動微分方程為：

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

式中：每個能量值（，）對應于特定的橫搖運動。例如：當=23.73、=19.74 時，得到對應的能量等值線，如圖3所示。

圖3 能量等值線

運動方程進一步寫為：

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

則（11）可進一步寫為如下隨機微分方程式：

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

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

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

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

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

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

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

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

式（25）可進一步寫為：

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

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

對于低強度噪聲激勵有=0，則式（29）可進一步寫為：

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

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

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

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

1.3 波群分析方法

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

同時，當≤0.7，f / f=5~15，當＞0.7，f / f=10~28。

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

圖5 波群模擬流程圖

JONSWAP公式如下所示：

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

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

2 船舶騎浪運動非線性分析

參照圖7坐標系，船舶在波浪中縱蕩運動方程為：

圖7 運動坐標系

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

船舶靜水阻力()為：

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

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

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

基于切片理論，船舶在波浪中運動時受到的波浪力可表達為：

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

式中：c可用多項式表達如下：

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

2.1 規(guī)則波中騎浪運動分岔研究

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

將式（52）中的右端項取不定積分，并由第1項減去第2項，則得到方程（52）對應哈密頓的量（，）：

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

圖8 異宿軌道示意圖

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

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

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

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

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

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

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

圖9 螺旋槳臨界轉速隨波幅變化曲線圖

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

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

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

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

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

圖11 轉速3.4 r/s時相對速度時間歷程曲線

圖12 轉速3.4 r/s時相對位移時間歷程曲線

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

圖13 轉速3.5 r/s時相對速度時間歷程曲線

圖14 轉速3.5 r/s時相對位移時間歷程曲線

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

圖15 3.4 r/s船舶運動相圖（未騎浪）

圖16 3.5 r/s船舶運動相圖（騎浪發(fā)生）

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

圖17 船舶相對速度隨螺旋槳轉速變化分岔圖

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

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

通過線性疊加求解隨機波浪力：

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

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

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

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

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

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

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

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

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

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

圖19 不同轉速下的騎浪概率

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

基于累積艏搖運動原理，建立船舶艏搖運動方程如下：

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

引入時間尺度=，將其代入式（71）得到如下表達式：

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

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

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

假設方程（75）的解為：

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

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

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

將公式（80）的近似解表達為：

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

設方程（82）的解分別為：

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

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

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

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

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

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

圖24 波高1.0 m艏搖角時間歷程曲線圖

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

圖22 波高0.5 m艏搖角時間歷程曲線圖

圖23 波高0.5 m船舶艏搖運動相圖

圖25 波高1.0 m船舶艏搖運動相圖

4 結 論

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

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

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

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

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

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