基于隨機Melnikov方法的甲板上浪船舶混沌運動研究

劉利琴，唐友剛

（天津大學 水利工程仿真與安全國家重點實驗室，天津300072）

1 引 言

甲板上浪后船舶的運動非常復雜，涉及船舶的大幅非線性運動、船舶與波浪之間的非線性耦合運動以及海洋環境的隨機因素。國內外學者基于數值模擬的方法研究了規則波浪中甲板上浪船舶的運動特性。Dillingham[1]最早應用淺水波理論，將甲板劃分成網格，數值模擬了甲板上浪船舶的橫搖、橫蕩耦合運動。黃祥鹿[2]將淺水波理論與切片理論相結合，在時域中分析比較了甲板上浪和無上浪兩種情況船舶的橫搖運動，結果表明，甲板上浪會導致具有小初穩性高的船舶傾覆。Belenky[3]基于淺水波假設，綜合計算了不同甲板上浪行為時船舶的橫搖運動，研究發現，甲板上浪水將增加船舶的阻尼，并引起船舶的亞諧運動。Laranjinha[4]應用隨機選取法求解甲板上浪水的三維運動，數值模擬了船舶六個自由度的響應，結果表明，小量的甲板上浪水將增加船舶的阻尼，并引起船舶的亞諧運動。然而，實際海洋環境是非規則的，研究隨機海浪中甲板上浪船舶的響應更具有實際意義。Yim[5]用概率密度函數研究了甲板上浪船舶的隨機混沌運動，發現異宿軌道與傾覆相聯系。Nakhata Tongchate[6]基于Yim的工作，在時域中數值模擬了甲板上浪船舶橫搖、垂蕩、縱搖耦合的三自由度響應。Troesch[7]將甲板上浪近似為一階靜水力，用概率統計的方法計算了隨機波浪中船舶的非線性橫搖運動，以及甲板上浪導致的傾覆問題。

本文考慮甲板上浪對船舶靜穩性的影響，建立隨機橫浪中船舶橫搖運動方程。將隨機橫浪激勵簡化為周期激勵加白噪聲擾動，應用隨機Melnikov方法和龐加萊截面研究甲板上浪時船舶的混沌運動。

2 運動方程的建立

隨機橫浪中甲板上浪時船舶橫搖運動的微分方程可以寫為：

其中，I44為船舶橫搖慣性矩；A44為附連水引起的附加慣性矩；B44為船舶橫搖線性阻尼系數；B44q為船舶橫搖非線性阻尼系數；Δ為船舶排水量；GZ（φ）為船舶的靜穩性臂；Fsea（t）為隨機橫浪引起的干擾力矩；Fwod（t）為甲板上浪引起的干擾力矩。

（1）式兩邊除以 （I44+ A44），并將隨機橫浪激勵簡化為周期激勵加白噪聲擾動，有：

考慮甲板上浪對船舶靜穩性的影響，（2）式進一步寫為：

其中，GZm（φ）為計入甲板上浪后船舶的恢復力臂曲線。分別用x和y表示橫搖角φ和橫搖角速度φ˙，將（3）式寫成以x和y表示的隨機微分方程的典型形式為：

3 隨機Melnikov方法

Melnikov函數具有簡單零點，意味著穩定流行與不穩定流行在龐加萊截面上橫截相交，系統出現Smale變化意義下的混沌。對于隨機動力學系統，需要從統計意義上討論隨機Melnikov過程是否有簡單零點，將隨機Melnikov過程和均方準則相結合確定系統的混沌運動的參數域[8]，獲得發生混沌時船舶的形狀參數與波浪參數間的臨界關系。

與方程（5）對應的無阻尼、未擾動方程為：

用-t來代替 t，（7）式可以寫成：

可以看出，Md和Mp是確定量，與隨機激勵無關，由方程（5）的阻尼和簡諧激勵產生，而Z是隨機過程，由方程（5）的隨機噪聲引起。隨機Melnikov過程（8）的均值為E[ M（ t1,t2）]=-Md+Mp（t1），方差為 σ2Z，計算公式為：

對于一般的船舶，GZm（x）是關于x的高次多項式，由（6）式很難求出橫搖角速度y0的解析表達式，本文采用數值法求解。由于y0是時間t（ t∈（-∞,+∞））的函數，（x0, y0）=（x0（ t）,y0（t））為同宿軌道，初始積分點q2選為同宿軌道與x軸的交點，從而y0（t）為t的奇函數。（8）式可進一步寫為：

其中，Mhom為同宿軌道的隨機Melnikov過程，A1，A3以及B（ω）的表達式如下：所以當y0給定時，A1和A3為常數，B為波浪頻率ω的函數。隨機Melnikov過程（8）在均方意義下出現零點的條件是：

上式可進一步寫為：

隨機Melnikov過程具有簡單零點，是隨機混沌運動的必要條件。本文應用（14）式確定船舶產生混沌運動的參數域，并結合龐加萊截面和響應歷程研究甲板上浪船舶的混沌運動。

4 算 例

算例為一艘拖網漁船，該船于1972年制造，1974年在北挪威諾爾辰角西部島海面附近失事，船舶的主要參數如表1所示[9]。

該拖網漁船具有雙層甲板，即捕魚甲板和主甲板，上層建筑設在船首，在船尾進行捕魚作業，漁獲物直接卸到下層主甲板的加工間進行處理。在船兩側的舷墻上，有24個排水孔，其中右舷13個，左舷11個；在右舷距設計水線1.6m處，有兩個0.53m長、0.46m寬的斜道，用來排除加工間產生的雜物，船舶具體結構參見文獻[10]。在一定條件下，海水會從排水孔、舷墻頂部以及右舷側的兩個斜道進入甲板并累積[10]。事故發生后，很多機構和學者分析了船舶失事的原因，通過對沉船的實際探測，推測該船的傾覆可能與甲板上浪有關[10]，目前對這一推測還沒有從理論上給予足夠的證實。下面以該拖網漁船為例，研究甲板上浪時船舶的混沌運動。

將恢復力臂曲線GZm擬合成十三次多項式，代入方程（4），有：

表1 船舶主要參數Tab.1 Principal particulars of the vessel

其中，ki（i=1,3,…,13）分別為線性和非線性恢復力矩系數。根據實驗[9]，對于一定量的甲板上浪，與方程（15）對應的船舶參數為：k1=-0.331 1,k3=2.644 6,k5=-6.5698,k7=8.242 1,k9=-5.339 4,k11=1.708 1,k13=-0.214 3,d1=0.056,d3=0.165 9,β0=0.8，參數 E、ω和D依據不同的海況而定。

當周期波浪頻率ω=0.5rad/s時，用（14）式計算不同白噪聲強度D時船舶發生混沌運動的阻尼系數d1與波浪幅值E間的臨界關系，結果如圖1所示。

圖1 船舶混沌運動的參數區域（ω=0.5rad/s）Fig.1 The parameter domain for onset of chaotic motion of ship(ω=0.5rad/s)

圖1中橫坐標為波浪幅值E，縱坐標為線性阻尼系數d1，圖中曲線為混沌域和非混沌域的分界線，不同的線型表示不同的白噪聲強度。由圖1可以看出，對于確定的波浪幅值E，增加船舶阻尼將抑制混沌運動的發生；增加白噪聲強度，將使非混沌參數域減小、混沌參數域增大；增加波浪幅值，船舶不發生混沌運動的臨界阻尼值將隨之增大。

以下計算船舶橫搖響應的龐加萊截面和時間歷程。龐加萊截面的構造方法如下：選定n個不同的初始條件，對每一個初始條件用四階Runge-Kutta法數值求解方程（15），每經過一個波浪周期取一個相點，計算200個周期，選后50個周期的相點構造龐加萊截面。

根據圖1劃分的參數區域，在混沌參數域中取兩組參數（d1=0.056,D=0,E=0.06,ω=0.5rad/s）和（d1=0.056,D=0.000 5,E=0.06,ω=0.5rad/s），在非混沌參數域中取兩組參數（d1=0.1,D=0,E=0.06,ω=0.5rad/s）和（d1=0.12,D=0.000 5,E=0.06,ω=0.5rad/s），分別對一百個初始點和一個初始點計算船舶的橫搖響應，構造龐加萊截面，結果如圖2～5所示。

圖2 龐加萊截面（d1=0.056,D=0,E=0.06,ω=0.5rad/s）Fig.2 Poincare’ maps(d1=0.056,D=0,E=0.06,ω=0.5rad/s)

圖3 龐加萊截面（d1=0.1,D=0,E=0.06,ω=0.5rad/s）Fig.3 Poincare’ maps(d1=0.1,D=0,E=0.06,ω=0.5rad/s)

圖4 龐加萊截面（d1=0.056,D=0.000 5,E=0.06,ω=0.5rad/s）Fig.4 Poincare’ maps(d1=0.056,D=0.000 5,E=0.06,ω=0.5rad/s)

圖5 龐加萊截面（d1=0.12,D=0.000 5,E=0.06,ω=0.5rad/s）Fig.5 Poincare’ maps(d1=0.12,D=0.000 5,E=0.06,ω=0.5rad/s)

由圖2和圖3可以看出，對于白噪聲強度D等于零的情況，當阻尼系數d1=0.056時船舶發生混沌橫搖運動，如圖2所示；當阻尼系數d1增大到0.1時，船舶橫搖為規則的周期運動，如圖3所示。由圖4和圖5可以看出，對于白噪聲強度D大于零的情況，當阻尼系數d1=0.056時船舶發生隨機混沌運動，如圖4所示；當阻尼系數d1增大到0.12時，船舶運動為隨機橫搖運動，如圖5所示。因此，增大橫搖阻尼抑制了船舶混沌運動的發生。

圖2～5還表明，甲板上浪時船舶橫搖響應的龐加萊截面上有兩個吸引域，船舶橫搖過程有兩個橫搖中心，分別位于左舷側和右舷側。在非混沌參數域，對于任意一組確定的初始條件，船舶圍繞其中的一個橫搖中心運動，如圖3（b）和圖5（b）所示；在混沌參數域，響應過程中發生由一個橫搖中心到另一個橫搖中心的隨機跳躍，如圖2（b）和圖4（b）所示。

用響應的時間歷程進一步驗證以上結論。對上面的四組參數，計算船舶橫搖響應的時間歷程，結果如圖6和圖7所示。

圖6 船舶橫搖響應歷程（D=0,E=0.06,ω=0.5rad/s）Fig.6 Roll response history of ship(D=0,E=0.06,ω=0.5rad/s)

圖7 船舶橫搖響應歷程（D=0.000 5,E=0.06,ω=0.5rad/s）Fig.7 Roll response history of ship(D=0.000 5,E=0.06,ω=0.5rad/s)

圖6和圖7中的橫搖響應歷程表明，在非混沌參數域，船舶圍繞某一個橫搖中心作小幅度橫搖，如圖6（a）和圖7（a）所示；在混沌參數域，響應過程中船舶將隨機地由一個橫搖中心跳躍到另一個橫搖中心，如圖 6（b）和圖 7（b）所示。

在混沌參數域中另取兩組參數，計算船舶橫搖響應的龐加萊截面，結果如圖8和圖9所示。

圖8 龐加萊截面（d1=0.056,E=0.1,ω=0.5rad/s）Fig.8 Poincare’ maps(d1=0.056,E=0.1,ω=0.5rad/s)

圖9 龐加萊截面（d1=0.056,E=0.12,ω=0.5rad/s）Fig.9 Poincare’ maps(d1=0.056,E=0.12,ω=0.5rad/s)

圖8和圖9表明，在混沌參數域中隨機噪聲使混沌吸引域的面積有所擴散。

5 結 語

考慮甲板上浪對船舶靜穩性的影響，應用隨機Melnikov方法和龐加萊截面研究隨機橫浪中甲板上浪時船舶的混沌運動。由隨機Melnikov過程結合均方準則確定產生混沌運動時船舶的阻尼系數與波浪參數間的臨界關系，計算不同參數域中船舶橫搖響應的龐加萊截面和時間歷程。結果表明，增加船舶阻尼將抑制混沌運動的發生；甲板上浪時船舶橫搖響應的龐加萊截面上有兩個吸引域，船舶運動過程有兩個橫搖中心；在非混沌參數域，船舶只圍繞其中的一個橫搖中心運動；在混沌參數域，響應過程中發生由一個橫搖中心到另一個橫搖中心的隨機跳躍。

[1]Dillingham J T,Falzarano J M.A numerical method for three-dimensional sloshing[C]//SNAME Spring Meeting Technical and Research Symposium(STAR).Portland,Oregon:1986:75-85.

[2]黃祥鹿,顧謝忡.船舶甲板上浪橫搖的時域模擬[J].中國造船,1993,26(3):27-36.

[3]Belenky V L,Liut D,et al.Nonlinear ship roll simulation with water-on-deck[C]//Proceedings of the Stability Workshop.New York,2002.

[4]Laranjinha M,Falzarano J M,et al.Analysis of the dynamical behavior of an offshore supply vessel with water on deck[C]//Proceedings of the 21st International Conference on Offshore Mechanics and Arctic Engineering(OMAE).Oslo,Norway,2002.

[5]Yim S C S,Lin H.Unified analysis of complex nonlinear motions via densities[J].Nonlinear Dynamics,2001,24:103-127.

[6]Nakhata Tongchate.Stability analysis of nonlinear coupled barge motions[D].Salem:Oregon State University,2003.

[7]McCue L,Troesch A.Probabilistic determination of critical wave height for a multi-degree of freedom capsize model[J].Ocean Engineering,2005,32(13):1608-1622.

[8]朱位秋.非線性隨機動力學與控制—Hamilton理論體系框架[M].北京:科學出版社,2003.

[9]Morrall A.The Gaul disaster:A investigation into the loss of a large stern trawler[J].Naval Architect,1981,5:391-440.

[10]Marine accident investigation branch.Report of the re-opened formal investigation into the loss of the FV Gaul,part 2:The subsea surveys from 1997 and the new evidence[R].Published by the Stationery Office,2004.