劉利琴，唐友剛

（天津大學(xué) 水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津300072）

1 引 言

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

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

2 運(yùn)動(dòng)方程的建立

隨機(jī)橫浪中甲板上浪時(shí)船舶橫搖運(yùn)動(dòng)的微分方程可以寫為：

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

（1）式兩邊除以 （I44+ A44），并將隨機(jī)橫浪激勵(lì)簡(jiǎn)化為周期激勵(lì)加白噪聲擾動(dòng)，有：

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

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

3 隨機(jī)Melnikov方法

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

與方程（5）對(duì)應(yīng)的無阻尼、未擾動(dòng)方程為：

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

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

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

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

上式可進(jìn)一步寫為：

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

4 算 例

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

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

將恢復(fù)力臂曲線GZm擬合成十三次多項(xiàng)式，代入方程（4），有：

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

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

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

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

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

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

根據(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），分別對(duì)一百個(gè)初始點(diǎn)和一個(gè)初始點(diǎn)計(jì)算船舶的橫搖響應(yīng)，構(gòu)造龐加萊截面，結(jié)果如圖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可以看出，對(duì)于白噪聲強(qiáng)度D等于零的情況，當(dāng)阻尼系數(shù)d1=0.056時(shí)船舶發(fā)生混沌橫搖運(yùn)動(dòng)，如圖2所示；當(dāng)阻尼系數(shù)d1增大到0.1時(shí)，船舶橫搖為規(guī)則的周期運(yùn)動(dòng)，如圖3所示。由圖4和圖5可以看出，對(duì)于白噪聲強(qiáng)度D大于零的情況，當(dāng)阻尼系數(shù)d1=0.056時(shí)船舶發(fā)生隨機(jī)混沌運(yùn)動(dòng)，如圖4所示；當(dāng)阻尼系數(shù)d1增大到0.12時(shí)，船舶運(yùn)動(dòng)為隨機(jī)橫搖運(yùn)動(dòng)，如圖5所示。因此，增大橫搖阻尼抑制了船舶混沌運(yùn)動(dòng)的發(fā)生。

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

用響應(yīng)的時(shí)間歷程進(jìn)一步驗(yàn)證以上結(jié)論。對(duì)上面的四組參數(shù)，計(jì)算船舶橫搖響應(yīng)的時(shí)間歷程，結(jié)果如圖6和圖7所示。

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

圖7 船舶橫搖響應(yīng)歷程（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中的橫搖響應(yīng)歷程表明，在非混沌參數(shù)域，船舶圍繞某一個(gè)橫搖中心作小幅度橫搖，如圖6（a）和圖7（a）所示；在混沌參數(shù)域，響應(yīng)過程中船舶將隨機(jī)地由一個(gè)橫搖中心跳躍到另一個(gè)橫搖中心，如圖 6（b）和圖 7（b）所示。

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

5 結(jié) 語(yǔ)

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

[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]黃祥鹿,顧謝忡.船舶甲板上浪橫搖的時(shí)域模擬[J].中國(guó)造船,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]朱位秋.非線性隨機(jī)動(dòng)力學(xué)與控制—Hamilton理論體系框架[M].北京:科學(xué)出版社,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.