周光耀，唐友剛，王麗元，李 妍

（1.天津大學(xué) 建筑工程學(xué)院，天津 300072；2.天津大學(xué) 水利工程與仿真國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津 300072）

引 言

參數(shù)橫搖是造成船舶的傾覆的主要原因之一，是近年來(lái)船舶穩(wěn)性和耐波性領(lǐng)域研究的熱點(diǎn)。

船舶參數(shù)激勵(lì)問(wèn)題早在1955年就被Kervin提出，20世紀(jì) 80年代，B.Roberts和 A.Bruce Dunwoody等首次對(duì)隨機(jī)參激橫搖進(jìn)行了概率分析[1,2]。1998年，巴拿馬型 C11級(jí)集裝箱船 APL CHINA號(hào)事故引起了學(xué)者對(duì)參數(shù)橫搖的重大關(guān)注[3]。由此國(guó)內(nèi)外學(xué)者開(kāi)始了對(duì)參數(shù)橫搖進(jìn)行了大量理論和實(shí)驗(yàn)方面的研究。早在20世紀(jì)90年代，天津大學(xué)唐友剛、鄭俊武、田凱強(qiáng)等就研究了船舶橫搖與縱搖、升沉的耦合[4,5]。李紅霞于2007年對(duì)船舶參數(shù)橫搖進(jìn)行了時(shí)歷運(yùn)動(dòng)計(jì)算和概率分析[6]。魯江對(duì)于縱向行駛在具有隨機(jī)性質(zhì)波浪中的航行體的橫向搖動(dòng)和波浪對(duì)于船舶的附加外力做出了數(shù)值上的模擬和探索[7,8]。2015年，張曉等從參數(shù)敏感性的角度對(duì)航行在具有隨機(jī)性質(zhì)波浪中的船舶的橫向搖擺進(jìn)行了探索[9]。同時(shí)國(guó)外學(xué)者也做了大量工作。Bulian等人對(duì)于船舶性質(zhì)參數(shù)和運(yùn)動(dòng)的隨機(jī)性質(zhì)進(jìn)行了探索[10,11]，他們的探索對(duì)象亦是縱向行駛在具有隨機(jī)性質(zhì)的波浪中的航行體；2009～2011年，Jovanoski Z、Atsuo Maki等均對(duì)縱向行駛在隨機(jī)環(huán)境中的航行體的運(yùn)動(dòng)概率特性和力學(xué)特性進(jìn)行了探索[12,13]。

從以上一些研究現(xiàn)狀來(lái)看，對(duì)于船舶在隨機(jī)波中的多自由度耦合的運(yùn)動(dòng)分析和船舶在隨機(jī)波中的橫搖角概率分布研究較少。且多是利用歐拉旋轉(zhuǎn)矩陣對(duì)船舶多自由度運(yùn)動(dòng)進(jìn)行研究。本文利用瞬時(shí)等效旋轉(zhuǎn)矩陣推導(dǎo)出考慮多自由度同時(shí)運(yùn)動(dòng)的運(yùn)動(dòng)方程，在此基礎(chǔ)上計(jì)算實(shí)時(shí)位置的波浪力。通過(guò)無(wú)理數(shù)分波浪譜模擬隨機(jī)波，計(jì)算實(shí)船在時(shí)歷上的運(yùn)動(dòng)響應(yīng)。并進(jìn)一步其進(jìn)行力學(xué)分析，分析運(yùn)動(dòng)不穩(wěn)定性和混沌產(chǎn)生的原因。同時(shí)對(duì)多組計(jì)算結(jié)果進(jìn)行直方圖統(tǒng)計(jì)和概率密度函數(shù)的模擬。創(chuàng)新引入風(fēng)險(xiǎn)評(píng)估與數(shù)學(xué)函數(shù)模擬相結(jié)合的思想，研究船舶在隨機(jī)波中的橫搖角概率。并深入研究，推導(dǎo)出研究其橫搖概率的廣義方法，擬合得到以波高和橫搖角為自變量的三維概率分布函數(shù)。以此估算任意波高下發(fā)生某個(gè)橫搖角度的概率。

1 坐標(biāo)系的建立及瞬時(shí)等效轉(zhuǎn)換關(guān)系

1.1 坐標(biāo)系的建立

坐標(biāo)系的選取與船舶運(yùn)動(dòng)和復(fù)原力計(jì)算有重大影響，為了方便表述，如圖1建立如下坐標(biāo)描述系統(tǒng)：

1）相對(duì)航行體平移的平移坐標(biāo)描述系統(tǒng)OUxUyUzU（equilibrium coordinate），記號(hào)U；

2）OTxTyTzT為聯(lián)動(dòng)坐標(biāo)描述系統(tǒng)（onboard coordinate），記為T，指坐標(biāo)描述系統(tǒng)隨著物體的旋轉(zhuǎn)運(yùn)動(dòng)而相應(yīng)地改變；

3）坐標(biāo)系Oxyz為恒定的坐標(biāo)描述系統(tǒng)（fixed coordinate），記為O，指坐標(biāo)描述系統(tǒng)不隨航行體的運(yùn)動(dòng)而改變的大地坐標(biāo)描述系統(tǒng)；

4）相對(duì)坐標(biāo)描述系ORxRyRzR（reference coordinate），記號(hào)R。

圖1 坐標(biāo)系示意

1.2 坐標(biāo)系的轉(zhuǎn)換關(guān)系

設(shè)A為旋轉(zhuǎn)矩陣用以表示T坐標(biāo)系中的點(diǎn)(x,y,z)與其在R坐標(biāo)系中(x1,y1,z1)的坐標(biāo)的關(guān)系為：

設(shè)P(t)為t時(shí)刻瞬時(shí)（Δt）坐標(biāo)轉(zhuǎn)換矩陣表示在Δt時(shí)刻內(nèi)的坐標(biāo)轉(zhuǎn)換關(guān)系，則有：

歐拉等效多自由度轉(zhuǎn)動(dòng)定律：航行體的運(yùn)動(dòng)可以通過(guò)矢量轉(zhuǎn)換分解為以X、Y、Z軸進(jìn)行旋轉(zhuǎn)的運(yùn)動(dòng)。加入航行體在三維空間中的旋轉(zhuǎn)角速度為ωx、ωy、ωz，則在每一個(gè)固定時(shí)刻航行體的旋轉(zhuǎn)運(yùn)動(dòng)都能夠用繞著某條直線l且直線經(jīng)過(guò)坐標(biāo)系統(tǒng)的原點(diǎn)的單軸旋轉(zhuǎn)運(yùn)動(dòng)，其單軸轉(zhuǎn)動(dòng)角速度為ωl，l與X、Y、Z軸分別成角度α、β、γ，且關(guān)系服從式(4)：

設(shè)λ=ωlΔt，則瞬時(shí)坐標(biāo)轉(zhuǎn)換矩陣為：

式(5)的矩陣確保了航行體在空間多自由度中坐標(biāo)位置的計(jì)算是同一時(shí)刻的。

2 多自由度耦合運(yùn)動(dòng)模型

聯(lián)動(dòng)坐標(biāo)系統(tǒng)的原點(diǎn)可以隨意選取，本文取航行體的質(zhì)量中心。根據(jù)航行體的結(jié)構(gòu)特性，我們能夠得到xG=yG=zG=Ixy=Iyz=0?？紤]航行體在空間中的運(yùn)動(dòng)的復(fù)雜特性和變系數(shù)微分方程的性質(zhì)，通過(guò)基本動(dòng)力學(xué)的推導(dǎo)，筆者建立了航行體的四自由度的運(yùn)動(dòng)方程為：

其中，∑Mxx表示關(guān)于T坐標(biāo)系x軸的波浪力的合力矩，∑Myy和∑Mzz類似。∑FOz表示O坐標(biāo)系中z方向上的合外力，I表示轉(zhuǎn)動(dòng)慣量。此運(yùn)動(dòng)模型不再用歐拉角進(jìn)行處理，而是通過(guò)式(3)和式(5)直接計(jì)算實(shí)時(shí)位置。

3 隨機(jī)縱浪及其波浪力的處理

設(shè)船航行方向?yàn)楣潭ㄗ鴺?biāo)系X軸正方向。對(duì)于不規(guī)則波的模擬利用譜密度函數(shù)S(ω)，采用下式[14]：

由于各頻率是等差數(shù)列，在長(zhǎng)時(shí)間后模擬的波形可能重復(fù)。一些處理方法是將頻率按波浪譜等面積劃分，在實(shí)際編程中需用到多次迭代計(jì)算各頻率。引入無(wú)理數(shù)采用較為簡(jiǎn)單的劃分方式，如下：

其中k為某無(wú)理數(shù)。取k為作為一個(gè)劃分倍數(shù)。

對(duì)(7)式可得一般的波中水壓強(qiáng)為：

其中，ρ為水密度，g為重力加速度。

4 計(jì)算實(shí)例與分析

4.1 運(yùn)動(dòng)結(jié)果及力學(xué)分析

采用某C11集裝箱船為實(shí)例研究對(duì)象，其主要參數(shù)如表1。

表1 C11集裝船主要參數(shù)

用JONSWAP譜模擬基于隨機(jī)性質(zhì)的波浪的譜密度函數(shù)，隨機(jī)參數(shù)組數(shù)為10，通過(guò)正余弦函數(shù)的線性和表示基于隨機(jī)性質(zhì)的波浪在空間中的形狀函數(shù)。改變隨機(jī)因子，對(duì)波長(zhǎng)接近船長(zhǎng)的不同波高的波浪下的運(yùn)動(dòng)進(jìn)行了多組計(jì)算。用隨船坐標(biāo)系的z軸與參考系的Oxz平面夾角 Angle表示橫搖運(yùn)動(dòng)。只給出各個(gè)波高下的一組時(shí)歷運(yùn)動(dòng)結(jié)果，如圖2～圖5。

圖2 波高2.5 m隨機(jī)波橫搖運(yùn)動(dòng)響應(yīng)

圖3 波高3.5 m隨機(jī)波橫搖運(yùn)動(dòng)響應(yīng)

圖4 波高4 m隨機(jī)波橫搖運(yùn)動(dòng)響應(yīng)

圖5 波高5 m隨機(jī)波橫搖運(yùn)動(dòng)響應(yīng)

在運(yùn)動(dòng)的計(jì)算中，我們選取了262 m為具有隨機(jī)性質(zhì)的波浪的波長(zhǎng)，縱向行駛的航行體在具有隨機(jī)性質(zhì)的波浪中的航行速度是2 kn，航行體在零時(shí)刻的橫向轉(zhuǎn)動(dòng)的角速度是0.01 rad/s。隨機(jī)波浪在零時(shí)刻的相位是通過(guò)計(jì)算機(jī)隨機(jī)產(chǎn)生的。

從圖2～圖5的計(jì)算結(jié)果可以看出，船舶在隨機(jī)波下的橫搖運(yùn)動(dòng)是不穩(wěn)定的。圖5和圖4比較可以看出，縱向行駛在具有隨機(jī)性質(zhì)的波浪的航行體的橫向搖擺的幅度最大值和波浪幅度并非是正相關(guān)的，這是因?yàn)椴ɡ说碾S機(jī)性質(zhì)存在一定的不確定性和奇異性。但也能夠看出航行體在整個(gè)橫向搖擺系統(tǒng)中的能量與波浪幅度是正相關(guān)的。

為了研究其非線性運(yùn)動(dòng)的不規(guī)則性和運(yùn)動(dòng)能量的分布情況，給出圖2～圖5對(duì)應(yīng)的相圖如圖6～圖9。

圖6 波高2.5 m隨機(jī)波橫搖運(yùn)動(dòng)相圖

圖7 波高3.5 m隨機(jī)波橫搖運(yùn)動(dòng)相圖

圖8 波高4 m隨機(jī)波橫搖運(yùn)動(dòng)相圖

圖9 波高5 m隨機(jī)波橫搖運(yùn)動(dòng)相圖

從圖6～圖9可以看出，其運(yùn)動(dòng)速度極值點(diǎn)在波高2 m是0.1 rad/s左右，到5 m波高時(shí)達(dá)到0.3 m/s左右。這表明其運(yùn)動(dòng)的能量是隨波高增大而增大的，這個(gè)現(xiàn)象從橫搖幅值不能完全斷言。船舶橫搖運(yùn)動(dòng)的能量，主要還是來(lái)源于波浪的能量，其中可能有一些轉(zhuǎn)換過(guò)程，但最終表現(xiàn)出來(lái)的是波浪的能量帶動(dòng)了船舶的橫搖。從圖6～圖9各相圖可以看出，在隨機(jī)波下的參數(shù)橫搖運(yùn)動(dòng)很容易就產(chǎn)生了混沌。觀察圖8和圖9可以發(fā)現(xiàn)，運(yùn)動(dòng)變得劇烈，相圖更是不規(guī)則，其運(yùn)動(dòng)越是不穩(wěn)定。其產(chǎn)生的主要原因來(lái)自于，大波高情況下船體表面的不規(guī)則性表現(xiàn)得更加明顯，加上隨機(jī)波的難以預(yù)料性，使其運(yùn)動(dòng)更為不穩(wěn)定。對(duì)于各相圖的集中區(qū)域可以看出，其集中區(qū)域是在不斷擴(kuò)大的，這是其分叉運(yùn)動(dòng)的一個(gè)縮影，也是其運(yùn)動(dòng)能量轉(zhuǎn)換的一個(gè)標(biāo)志。

為了研究其混沌運(yùn)動(dòng)的內(nèi)在頻率關(guān)系，給出各運(yùn)動(dòng)的頻譜圖，如圖10～圖13。

圖10 波高2.5 m隨機(jī)波橫搖運(yùn)動(dòng)頻譜

圖11 波高3.5 m隨機(jī)波橫搖運(yùn)動(dòng)頻譜

圖12 波高4 m隨機(jī)波橫搖運(yùn)動(dòng)頻譜

圖13 波高5 m隨機(jī)波橫搖運(yùn)動(dòng)頻譜

其中ω表示圓頻率。船舶橫搖固有頻率為0.26 rad/s。在各波高下的運(yùn)動(dòng)的主頻中都出現(xiàn)了接近固有頻率的主頻。此頻率在能量方面對(duì)于水面航行體的橫向搖擺運(yùn)動(dòng)的影響隨著具有隨機(jī)性質(zhì)的波浪的波幅增大而減小。同時(shí)，從功率譜圖可以看出水面航行體的橫向搖擺運(yùn)動(dòng)具有多個(gè)復(fù)雜的頻率，這些復(fù)雜的頻率來(lái)自于隨機(jī)性質(zhì)的波浪中占能量的比例相對(duì)較大的正余弦函數(shù)的頻率。然而還有一些偏離了隨機(jī)波浪中的正余弦函數(shù)的頻率，這些值來(lái)自于船舶運(yùn)動(dòng)中的多種不確定性，如航行體和波的不規(guī)則性、耦合影響和仿真中數(shù)值不穩(wěn)定性。

4.2 隨機(jī)運(yùn)動(dòng)結(jié)果的統(tǒng)計(jì)分布擬合和三維概率分析

由于船舶復(fù)原力不能用一個(gè)解析的表達(dá)式去相對(duì)準(zhǔn)確的模擬出來(lái)，所以這些不規(guī)則的頻率是無(wú)法通過(guò)一個(gè)完整的解析表達(dá)出來(lái)的，因此對(duì)于其隨機(jī)波橫搖運(yùn)動(dòng)的判定就需要通過(guò)數(shù)值的模擬。以下將通過(guò)隨機(jī)統(tǒng)計(jì)的方法對(duì)航行體的橫向搖擺運(yùn)動(dòng)性質(zhì)進(jìn)行概率分析。每次計(jì)算選取諧波的零時(shí)刻的相位的隨機(jī)種子個(gè)數(shù)為10，相位是由計(jì)算機(jī)隨機(jī)產(chǎn)生的。對(duì)各個(gè)波高下的船舶橫搖運(yùn)動(dòng)都進(jìn)行多組計(jì)算，得到不同初始相位下的運(yùn)動(dòng)計(jì)算結(jié)果。提取各橫搖時(shí)歷運(yùn)動(dòng)的幅值，統(tǒng)計(jì)其絕對(duì)值的直方圖。下面使用左端點(diǎn)接近于0的log-logistics分布對(duì)在航行在具有隨機(jī)性的波浪的航行體的橫向搖擺運(yùn)動(dòng)的角度進(jìn)行數(shù)值擬合。如式（10）：

其中：μ表示位置參數(shù)（location），相當(dāng)于分布的均值；σ表示尺度參數(shù)（scale），相當(dāng)于分布的標(biāo)準(zhǔn)差；λ表示閾值參數(shù)(threshold)，表示分布的左極限。在研究中均取接近于0為其左極限。

每個(gè)波高下的統(tǒng)計(jì)結(jié)果和分布擬合，作出其百分比（percent）和角度（angel）關(guān)系如圖14～圖17。

圖14 波高2.5 m統(tǒng)計(jì)擬合分布

圖15 波高3.5 m統(tǒng)計(jì)擬合分布

圖16 波高4 m統(tǒng)計(jì)擬合分布

圖17 波高5 m統(tǒng)計(jì)擬合分布

其中，loc表示參數(shù)μ，scale表示參數(shù)σ，thresh表示參數(shù)λ。

這樣一組統(tǒng)計(jì)擬合分布圖可以用來(lái)預(yù)測(cè)相應(yīng)波高下的船舶橫搖某個(gè)角度的發(fā)生概率。從統(tǒng)計(jì)擬合分布圖中可以看出，其概率密度函數(shù)的集中區(qū)域隨著隨機(jī)波的波高增大而右移。這表明隨波高增大，出現(xiàn)大角度橫搖角的概率越來(lái)越大，并存在一定概率達(dá)到極限橫搖角而發(fā)生傾覆。從圖中亦可發(fā)現(xiàn)，在小波高情況下，船舶很難發(fā)生大幅橫搖，在40°左右概率密度就接近于0了。在實(shí)際工程中船舶存在一個(gè)入水角，故概率分析考慮了此入水角為概率的吸收壁，即到達(dá)此角度時(shí)橫搖角幅值將停留于危險(xiǎn)狀態(tài)不再回去。所以這就需要求得某個(gè)角度下的發(fā)生概率即求其分布函數(shù)。利用(10)式求圖14～圖17中的擬合結(jié)果的分布函數(shù)，并作出其概率（probability）和角度（angle）的分布函數(shù)，如圖18。

圖18 各波高分布函數(shù)示意

從圖18可以看到在等角度的情況下，水面航行體的橫向搖擺概率與具有隨機(jī)性質(zhì)的波浪的能量是正相關(guān)的。對(duì)于概率密度函數(shù)的參數(shù)繼續(xù)進(jìn)行擬合。根據(jù)仿真計(jì)算結(jié)果和理論推導(dǎo)，可用下面的式子對(duì)概率密度函數(shù)的均值進(jìn)行近似表達(dá)：

其中bi表示擬合的參數(shù)。這組參數(shù)與航行體的結(jié)構(gòu)特性相關(guān)，若這組參數(shù)被確定，則航行體在隨機(jī)復(fù)雜環(huán)境中的運(yùn)動(dòng)的概率性規(guī)律就可以準(zhǔn)確描述，甚至（11）式的形式并不唯一，亦可通過(guò)大量的仿真數(shù)據(jù)對(duì)其進(jìn)行訓(xùn)練從而進(jìn)行調(diào)整。而這組參數(shù)是可以通過(guò)機(jī)器學(xué)習(xí)進(jìn)行不斷的訓(xùn)練而得的，從而能夠達(dá)到準(zhǔn)確預(yù)測(cè)水面航行體在不確定環(huán)境下的運(yùn)動(dòng)規(guī)律；H是隨機(jī)波浪的波高有義值。引入機(jī)器訓(xùn)練的概念是另一個(gè)計(jì)算機(jī)科學(xué)和神經(jīng)網(wǎng)絡(luò)科學(xué)與水面航行體運(yùn)動(dòng)預(yù)測(cè)相結(jié)合的一個(gè)重要的研究方向，本文還是側(cè)重于航行體在具有隨機(jī)性質(zhì)的波浪中的運(yùn)動(dòng)的分析方法框架和理論模型的建立，不在于參數(shù)訓(xùn)練原則和方法，因此不作贅述。對(duì)參數(shù)做出簡(jiǎn)易的訓(xùn)練處理，可得：b1=0.013928;b2= 0 .660010。結(jié)合參數(shù)訓(xùn)練結(jié)果做出概率分布函數(shù)的三維視圖如圖19。

圖19 波高、橫搖角三維分布函數(shù)

圖19是基于統(tǒng)計(jì)概率學(xué)的模型對(duì)描述水面航行體的運(yùn)動(dòng)的概率公式進(jìn)行參數(shù)分析而得到的解果，這個(gè)結(jié)果的價(jià)值在于可以對(duì)水面航行體在具有隨機(jī)性質(zhì)的環(huán)境中的運(yùn)動(dòng)進(jìn)行概率預(yù)測(cè)。同時(shí)，引入概率學(xué)模型和計(jì)算機(jī)科學(xué)機(jī)器學(xué)習(xí)的思路對(duì)于宏觀公式中的參數(shù)進(jìn)行訓(xùn)練從而達(dá)到描述和預(yù)測(cè)海洋港口工程項(xiàng)目中的對(duì)象的運(yùn)動(dòng)規(guī)律的思路，是一個(gè)跨學(xué)科理論體系的建立，對(duì)于后續(xù)的海洋港口中的科學(xué)系統(tǒng)研究可以起到很好的指導(dǎo)作用。

5 結(jié) 語(yǔ)

1）基于瞬時(shí)轉(zhuǎn)換矩陣考慮坐標(biāo)軸的瞬變性建立了多自由度耦合非線性運(yùn)動(dòng)方程?？紤]控制航向，對(duì)運(yùn)動(dòng)方程進(jìn)行退化。基于勢(shì)流理論推導(dǎo)波浪力的計(jì)算方法，并通過(guò)編程實(shí)現(xiàn)考慮瞬時(shí)位置的波浪力的實(shí)時(shí)計(jì)算。對(duì)于隨機(jī)波利用無(wú)理數(shù)倍數(shù)關(guān)系劃分波浪譜。最終實(shí)現(xiàn)對(duì)船舶在隨機(jī)縱浪中的多自由度耦合下的參數(shù)橫搖的研究。

2）研究船舶在隨機(jī)波中的參數(shù)橫搖運(yùn)動(dòng)發(fā)現(xiàn)，船舶在隨機(jī)波中的參數(shù)橫搖總是不穩(wěn)定的混沌運(yùn)動(dòng)。其分叉區(qū)域在不斷擴(kuò)大。在某些情況下船舶橫搖最大幅值不一定隨隨機(jī)波的波高增大而增大，但其能量大致是隨波高增大的。

3）影響船舶在隨機(jī)波中的參數(shù)橫搖運(yùn)動(dòng)的主要頻率來(lái)自于船舶的固有橫搖頻率和由二階耦合項(xiàng)產(chǎn)生的頻率。同時(shí)伴隨著由船表面不規(guī)則性和波浪的隨機(jī)性引起的其它一些無(wú)規(guī)律的頻率，這些頻率最終導(dǎo)致船舶參數(shù)橫搖運(yùn)動(dòng)產(chǎn)生了不穩(wěn)定性和混沌性。

4）通過(guò)統(tǒng)計(jì)和分布模擬對(duì)船舶在某個(gè)特定波長(zhǎng)下的隨機(jī)波中的橫搖角發(fā)生的概率分布進(jìn)行擬合和分析。并且以波高為自變量對(duì)其進(jìn)行繼續(xù)擬合，最終得到波高、橫搖角下三維概率分布。發(fā)現(xiàn)這個(gè)三維分布曲面是隨波高增大而向下滑移的，這符合了在大波高情況下大橫搖角發(fā)生的概率更大。首次引入三參數(shù) log-logistics分布進(jìn)行擬合，其特性符合船舶在隨機(jī)波中的參數(shù)橫搖時(shí)角度發(fā)生的概率分布。通過(guò)這種分布模擬出的三維分布可以估算此波長(zhǎng)下任意波高下的某個(gè)橫搖角發(fā)生的概率。