劉利琴，唐友剛

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

1 引 言

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

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

2 運動方程的建立

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

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

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

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

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

3 隨機Melnikov方法

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

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

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

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

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

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

上式可進一步寫為：

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

4 算 例

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

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

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

表1 船舶主要參數(shù)Tab.1 Principal particulars of the vessel

其中，ki（i=1,3,…,13）分別為線性和非線性恢復力矩系數(shù)。根據(jù)實驗[9]，對于一定量的甲板上浪，與方程（15）對應的船舶參數(shù)為：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，參數(shù) E、ω和D依據(jù)不同的海況而定。

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

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

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

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

根據(jù)圖1劃分的參數(shù)區(qū)域，在混沌參數(shù)域中取兩組參數(shù)（d1=0.056,D=0,E=0.06,ω=0.5rad/s）和（d1=0.056,D=0.000 5,E=0.06,ω=0.5rad/s），在非混沌參數(shù)域中取兩組參數(shù)（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等于零的情況，當阻尼系數(shù)d1=0.056時船舶發(fā)生混沌橫搖運動，如圖2所示；當阻尼系數(shù)d1增大到0.1時，船舶橫搖為規(guī)則的周期運動，如圖3所示。由圖4和圖5可以看出，對于白噪聲強度D大于零的情況，當阻尼系數(shù)d1=0.056時船舶發(fā)生隨機混沌運動，如圖4所示；當阻尼系數(shù)d1增大到0.12時，船舶運動為隨機橫搖運動，如圖5所示。因此，增大橫搖阻尼抑制了船舶混沌運動的發(fā)生。

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

用響應的時間歷程進一步驗證以上結論。對上面的四組參數(shù)，計算船舶橫搖響應的時間歷程，結果如圖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中的橫搖響應歷程表明，在非混沌參數(shù)域，船舶圍繞某一個橫搖中心作小幅度橫搖，如圖6（a）和圖7（a）所示；在混沌參數(shù)域，響應過程中船舶將隨機地由一個橫搖中心跳躍到另一個橫搖中心，如圖 6（b）和圖 7（b）所示。

在混沌參數(shù)域中另取兩組參數(shù)，計算船舶橫搖響應的龐加萊截面，結果如圖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表明，在混沌參數(shù)域中隨機噪聲使混沌吸引域的面積有所擴散。

5 結 語

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

[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.