隨機縱浪中多自由度耦合的船舶橫搖運動研究

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

（1.天津大學(xué) 建筑工程學(xué)院，天津 300072；2.天津大學(xué) 水利工程與仿真國家重點實驗室，天津 300072）

引 言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 多自由度耦合運動模型

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

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

3 隨機縱浪及其波浪力的處理

設(shè)船航行方向為固定坐標(biāo)系X軸正方向。對于不規(guī)則波的模擬利用譜密度函數(shù)S(ω)，采用下式[14]：

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

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

對(7)式可得一般的波中水壓強為：

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

4 計算實例與分析

4.1 運動結(jié)果及力學(xué)分析

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

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

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

圖2 波高2.5 m隨機波橫搖運動響應(yīng)

圖3 波高3.5 m隨機波橫搖運動響應(yīng)

圖4 波高4 m隨機波橫搖運動響應(yīng)

圖5 波高5 m隨機波橫搖運動響應(yīng)

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

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

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

圖6 波高2.5 m隨機波橫搖運動相圖

圖7 波高3.5 m隨機波橫搖運動相圖

圖8 波高4 m隨機波橫搖運動相圖

圖9 波高5 m隨機波橫搖運動相圖

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

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

圖10 波高2.5 m隨機波橫搖運動頻譜

圖11 波高3.5 m隨機波橫搖運動頻譜

圖12 波高4 m隨機波橫搖運動頻譜

圖13 波高5 m隨機波橫搖運動頻譜

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

4.2 隨機運動結(jié)果的統(tǒng)計分布擬合和三維概率分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 語

1）基于瞬時轉(zhuǎn)換矩陣考慮坐標(biāo)軸的瞬變性建立了多自由度耦合非線性運動方程。考慮控制航向，對運動方程進行退化?；趧萘骼碚撏茖?dǎo)波浪力的計算方法，并通過編程實現(xiàn)考慮瞬時位置的波浪力的實時計算。對于隨機波利用無理數(shù)倍數(shù)關(guān)系劃分波浪譜。最終實現(xiàn)對船舶在隨機縱浪中的多自由度耦合下的參數(shù)橫搖的研究。

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

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

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