馬騁遠(yuǎn)，馬 寧，王廷昊，顧解忡

(1. 上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海 200240；2. 上海交通大學(xué) 海洋工程國家重點實驗室，上海 200240；3. 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

0 引 言

IMO正在制定包括參數(shù)橫搖、純穩(wěn)性喪失、騎浪/橫甩、癱船穩(wěn)性和過度加速度這5種穩(wěn)性失效模式的第2代完整穩(wěn)性規(guī)則[1]。具有傳統(tǒng)意義的癱船穩(wěn)性是指船舶處于失去動力無法推進(jìn)或操舵狀態(tài)下，當(dāng)船自由漂移時在風(fēng)浪聯(lián)合作用下發(fā)生共振橫搖甚至傾覆[2]。IMO目前正在制定的癱船穩(wěn)性薄弱性衡準(zhǔn)草案主要以意大利和日本的提案為基礎(chǔ)，分為第1層衡準(zhǔn)和第2層衡準(zhǔn)[3–4]。針對第1層衡準(zhǔn)方法基本已達(dá)成一致，對于第2層衡準(zhǔn)方法的確定預(yù)計2018年左右才能完成[5]。

癱船穩(wěn)性失效中基于危險性評估的傾覆概率方法成為近年來研究的熱點，針對隨機(jī)激勵下復(fù)原力矩的非線性一直是處理該問題的主要難點[6]。國內(nèi)外學(xué)者針對橫浪狀態(tài)下的船舶大幅橫搖運動的穩(wěn)性進(jìn)行過許多研究，包括線性化理論和線性化處理方法、非線性動力學(xué)的應(yīng)用以及直接數(shù)值模擬方法[7]。自從1953年Denis和Pierson首次將隨機(jī)過程理論引入到線性耐波性領(lǐng)域的研究中，船舶在隨機(jī)風(fēng)浪中的傾覆概率研究也引起了各國學(xué)者的重視。

IMO草案中日本和意大利分別基于分段線性和局部線性方法，各自提出較為完善的概率計算方法。52屆SLF會議上，意大利代表建議將2008 IS Code中氣候衡準(zhǔn)作為癱船穩(wěn)性的第1層衡準(zhǔn)[8]，并提出了基于局部線性化方法的橫搖傾覆概率模型作為癱船穩(wěn)性的第2層薄弱性衡準(zhǔn)或直接評估方法[3]；日本代表提議利用基于Froude-Krylov假設(shè)的簡化公式求解有效波傾系數(shù)[9]，利用簡化的Ikeda方法求解橫搖阻尼系數(shù)[10]，同時將基于GZ曲線分段線性化處理的橫搖傾覆概率模型作為癱船穩(wěn)性的第2層衡準(zhǔn)方法[4]。Belenky[11]在1993年提出了船舶橫搖運動的分段線性化理論，通過對回復(fù)力矩的分段線性處理，求解得到單自由度橫搖運動方程的解析解，并以此評估傾覆概率。之后基于該方法，Iskandar[12]，Umeda[13]，Paroka[14–15]等日本學(xué)者進(jìn)行了大量實船評估計算和方法上的改進(jìn)。

本文基于癱船穩(wěn)性第2層衡準(zhǔn)，針對C11和CEHIPAR 2792兩艘集裝箱船，利用分段線性化方法得到橫搖運動方程的解析解，并據(jù)此計算傾覆概率。在樣船計算的基礎(chǔ)上同局部線性化結(jié)果進(jìn)行對比，并重點探討暴露時間、波浪譜以及分段分界點的選取等敏感性因素對傾覆概率的影響。

1 癱船穩(wěn)性分段線性化數(shù)值計算模型

1.1 單自由度橫搖運動方程

船舶在橫浪下運動，一般采用橫搖和橫蕩耦合的運動模型，并且考慮波浪輻射力和繞射力的影響。田才福造等[16]已經(jīng)證明單自由度橫搖方程也可以近似表達(dá)橫搖運動。因此，日本方法采用1-DOF非耦合運動方程，波浪激勵力只考慮Froude-Krylov力，方程無因次化后表示如下：

式中：ω0為橫搖固有頻率；G M為初穩(wěn)性高；為風(fēng)作用的無因次瞬時力矩；為波浪作用的無因次瞬時力矩。

無因次的Froude-Krylov波浪激勵力矩可由有效波傾系數(shù)計算得到。在線性假設(shè)下，將不同振幅和相位的規(guī)則波進(jìn)行疊加，研究船舶在不規(guī)則波中的搖蕩運動特性。波浪譜采用第15屆ITTC會議推薦的雙參數(shù)譜：

式中：

其中：H1/3為有義波高，m；Tm為不規(guī)則波的平均周期，s。

其中：

陣風(fēng)譜函數(shù)采用Davenport譜，表示如下：

非線性阻尼一般用隨機(jī)等效線性方法近似，橫搖阻尼系數(shù)α的取值按下式進(jìn)行：

圖1 風(fēng)速 m/s時波浪譜與陣風(fēng)譜密度函數(shù)Fig.1 Functions of wave and wind spectral density in wind speed m/s

1.2 復(fù)原力矩項分段線性處理

分段線性化方法最早由Belenky[11]在1993年提出，使得船舶在波浪中的傾覆概率可以通過解析解進(jìn)行研究。

GZ曲線的分段線性按照如下原則[6]進(jìn)行計算：分段線性三角形的高與原始GZ曲線的最大復(fù)原力臂保持一致，三角形的底邊長等于靜穩(wěn)性曲線的穩(wěn)距，三角形與坐標(biāo)軸所圍成的面積與靜穩(wěn)性曲線下的面積（傾斜后船舶具有的位能）近似相等，如圖2所示。

等效的穩(wěn)性高計算式如下：

圖2 分段線性化方法示意圖[6]Fig.2 Diagram of piecewise linearization method [6]

圖3 無因次化GZ曲線分段線性化示意圖Fig.3 Diagram of piecewise linearization method for dimensionless GZ curve

1.3 橫搖運動方程解析解

將1-DOF無因次化的橫搖運動方程（1）的各非線性系數(shù)線性近似后，將簡化為關(guān)于橫搖角的2階常系數(shù)線性非齊次微分方程。經(jīng)過計算，得到方程（1）在每個階段的解析解如下（式中各參數(shù)表達(dá)式詳見文獻(xiàn)[17]）：

1）Range 0階段

2）Range 1階段

2 傾覆概率計算模型

2.1 傾覆條件

分段線性化方法將船舶的橫搖運動分為2個區(qū)間進(jìn)行考慮：

對于Range 0階段的解析式（12），可以看出自由橫搖運動幅時間推移逐漸收斂，C船橫搖11角時歷如圖4所示，即認(rèn)為這種情況不屬于引起船舶傾覆的危險事件。

圖4 風(fēng)速 m/s下C11船設(shè)計吃水時橫搖角時歷變化（解析式（12）的圖像）Fig.4 Results of rolling angle’s variation in time-domain of C11 ship in wind speed m/s(graph of formula (12))

順風(fēng)條件下，A＞0為當(dāng)橫搖角從Range 0增大進(jìn)入Range 1后，橫搖運動幅值會隨著時間推移繼續(xù)增大，直至船舶傾覆。A＜0和A=0為當(dāng)橫搖角從Range 0增大進(jìn)入Range 1后，橫搖運動幅值會逐漸減小直至返回Range 0，這種情況安全。

逆風(fēng)條件下，A＞0和A=0趨于安全，A＜0代表船舶發(fā)生傾覆。

這樣船舶傾覆概率問題就轉(zhuǎn)換為求解橫搖運動解析解中的參數(shù)A的取值范圍。如果定義單位時間內(nèi)船舶傾覆為事件F，船舶橫搖進(jìn)入順風(fēng)Range 1階段為事件L，船舶橫搖進(jìn)入逆風(fēng)Range 1階段為事件W，A＞0為事件A+，A＜0為事件A-，則事件F發(fā)生概率可以寫成：

定義順風(fēng)和逆風(fēng)狀態(tài)下單位時間內(nèi)橫搖角從Range 0增大到Range 1的概率分別為和，由于在Range 1階段的橫搖運動符合高斯過程(Gaussian Process)，所以與分別由以下2個公式確定(Price & Bishop, 1974)[18]：

對于參數(shù)A取值范圍的概率問題，本文采用IMO推薦的簡化計算方法[17]，認(rèn)為橫搖運動進(jìn)入Range 1時的初始橫搖運動是一個小量，即等于0。此時參數(shù)A的概率密度函數(shù)可表達(dá)為：

最終通過計算：

即可得到參數(shù)A取值范圍的概率。由以上各式即可求解單位時間內(nèi)船舶的傾覆概率：

2.2 傾覆概率計算

經(jīng)過2.1節(jié)的公式推導(dǎo)，得到了單位時間傾覆概率計算表達(dá)式，假設(shè)船舶的傾覆事件符合泊松分布，可得到給定暴露時間T內(nèi)的傾覆概率：

3 樣船計算與敏感性因素分析

3.1 樣船模型

以C11集裝箱船（APL China號）和IMO穩(wěn)性工作組提供的研究船型CEHIPAR 2792集裝箱船為例，進(jìn)行分析與計算。2艘樣船的基本參數(shù)見表1與表2。取空氣阻尼系數(shù)Cm=1.0，線性阻尼系數(shù)為0.003 81 v/s–1。

3.2 GZ曲線分段線性化結(jié)果

分別對C11船與CEHIPAR 2792船在各自設(shè)計吃水下的GZ曲線進(jìn)行分段線性近似，如圖5所示。

3.3 傾覆概率敏感性因素分析

3.3.1 暴露時間

本節(jié)基于分段線性化方法，利用Matlab編制了相應(yīng)的計算程序，計算了CEHIPAR 2792船在設(shè)計吃水載況下，暴露時間h和年，定常風(fēng)速取范圍內(nèi)的傾覆概率。同時，將暴露時間為1 h的計算結(jié)果同局部線性化（具體計算過程見文獻(xiàn)[3]）結(jié)果進(jìn)行對比，如圖6所示。

表1 C11船基本參數(shù)Tab.1 Principle dimensions of C11 ship

表2 CEHIPAR 2792船基本參數(shù)Tab.2 Principle dimensions of CEHIPAR 2792 ship

圖5 兩艘樣船在各自設(shè)計吃水狀態(tài)下GZ曲線分段線性化結(jié)果Fig.5 Piecewise linearization method of two ships' GZ curve

通過分段線性化和局部線性化方法計算CEHIPAR 2792船的傾覆概率，都能夠得到一條具有當(dāng)定常風(fēng)速超過一定閾值后，船舶傾覆概率出現(xiàn)急劇上升特征的曲線。h的傾覆概率大幅增長對應(yīng)的風(fēng)速區(qū)間要大于年，局部線性化方法計算所得的傾覆概率突增對應(yīng)的風(fēng)速區(qū)間要遠(yuǎn)大于分段線性化結(jié)果，這可能與分段線性化處理后船舶傾斜位能（GZ曲線下的面積）的減小有關(guān)。

圖6 CEHIPAR 2792船傾覆概率曲線Fig.6 Capsizing probability’s curve of CEHIPAR 2792 ship

3.3.2 波浪譜

圖7（a）～圖7（d）所示為在定常風(fēng)速均為Um=25 m/s時4種常用的波浪譜函數(shù)圖像。

圖8所示為分別采用這4種波浪譜計算C11船在正常吃水下的傾覆概率。我國沿海波能譜的計算結(jié)果和其他3種譜的差別比較大，而且最為危險。應(yīng)用JONSWAP譜的計算結(jié)果對船舶抵抗傾覆能力的評估最為樂觀，ITTC雙參數(shù)譜和PM譜則介于中間。對于我國沿海波能譜、PM譜、ITTC雙參數(shù)譜，從三者波浪譜函數(shù)圖像來看，其形態(tài)相近，我國沿海波能譜波峰幅值最大，因此波浪激勵力矩也較大，相同定常風(fēng)速下傾覆概率計算結(jié)果最大；JONSWAP譜雖然波能譜峰值很高，但波浪頻率跨度很窄，波能很低，相同定常風(fēng)速下傾覆概率計算結(jié)果最小。

3.3.3 GZ曲線分段點位置

在1.2節(jié)中已經(jīng)簡要介紹了GZ曲線分段線性擬合的方式，但是對于分段三角形中回復(fù)力臂最大處橫傾角（即Range 0和Range 1的分界點）取法并沒有做出規(guī)定。計算表明，取值不同會影響傾覆概率的數(shù)值。圖9（a）和圖9（b）所示，分別為C11船和CEHIPAR 2792船在正常吃水下通過適當(dāng)平移分界點橫坐標(biāo)位置后得到的若干條傾覆概率曲線。

4 結(jié) 語

本文針對船舶在橫風(fēng)橫浪中癱船狀態(tài)的大幅橫搖運動，主要是基于IMO最新制定的癱船穩(wěn)性第2層衡準(zhǔn)開展分析，利用分段線性化方法得到橫搖運動方程的解析解，并據(jù)此進(jìn)行傾覆概率的計算。可以得到以下結(jié)論：

圖7 風(fēng)速為Um=25 m/s時各波浪譜函數(shù)的圖像Fig.7 Graphs of wave spectrums in wind speed Um=25 m/s

圖9 GZ曲線分段點位置對傾覆概率的影響Fig.9 Capsizing probability influenced by the position of GZ curve’s demarcation

1）GZ曲線分段線性近似法簡化了橫搖運動方程的非線性系數(shù)項，使得微分方程解析解的求解成為可能，且可以有效地計算船舶在橫風(fēng)橫浪下的傾覆概率，計算結(jié)果相比局部線性化方法偏于危險。

2）對于分段線性化的計算，暴露時間越長，傾覆概率的計算結(jié)果越危險。

3）不同的波浪譜對傾覆概率的計算影響是存在的，應(yīng)當(dāng)根據(jù)海區(qū)特征選取適當(dāng)?shù)牟ɡ俗V。

4）對某些具有凹凸變化的GZ曲線的船型（如CEHIPAR 2792船），分段線性化分界點的取法對其傾覆概率的評估影響較大，從提高穩(wěn)性衡準(zhǔn)安全裕度的角度，建議將其取在GZ曲線頂點位置處。