賴國銳 毛筱菲 詹星宇

(武漢理工大學(xué)船海與能源動力工程學(xué)院1) 武漢 430063) (武漢理工大學(xué)高性能艦船技術(shù)教育部重點實驗室2) 武漢 430063)

0 引 言

在隨機橫風(fēng)橫浪的共同作用下，失去動力的船舶會在產(chǎn)生橫蕩運動的同時橫搖至一個大角度，導(dǎo)致船上的開口進(jìn)水，造成船舶失穩(wěn)甚至傾覆.目前，IMO已經(jīng)制定了針對完整船舶癱船穩(wěn)性的第一層、第二層衡準(zhǔn)，并提出了癱船穩(wěn)性失效模式的直接評估指南.近年來，國內(nèi)外學(xué)者對癱船穩(wěn)性直接評估方法展開了大量研究.Belenky等[1]使用LAMP水動力學(xué)程序計算不規(guī)則橫浪中船舶的傾覆概率.Umeda等[2]研究了包含“橫蕩-垂蕩-橫搖-縱搖”的耦合數(shù)值模型，但缺乏對非定常風(fēng)的考慮.Kubo等[3]對Umeda的方法進(jìn)行了完善，建立了不規(guī)則橫風(fēng)和橫浪作用下的四自由度模型，并將計算結(jié)果與物理實驗結(jié)果進(jìn)行了比較.王新宇等[4]基于SDC1/INF.6的癱船薄弱性衡準(zhǔn)相關(guān)內(nèi)容，研究了癱船狀態(tài)的船舶因貨物移動而具有不同初始橫傾角時的傾覆概率變化情況.胡麗芬等[5]將破損船舶進(jìn)水過程的時域計算和傾覆概率計算相結(jié)合，探究了不同載況下破損船舶在進(jìn)水過程中穩(wěn)性高度和傾覆概率的關(guān)系.王田華等[6]建立了四自由度耦合數(shù)學(xué)模型及船舶水動力系數(shù)數(shù)據(jù)庫，運用非線性時域預(yù)報方法對船舶的整個傾覆過程進(jìn)行模擬計算.在2019年SDC 7-5會議中，IMO建議采用至少包含“橫蕩-垂蕩-橫搖-縱搖”四自由度耦合運動模型進(jìn)行船舶的運動模擬，但采用多自由度運動模型時為達(dá)到可信度要求需要計算的樣本量很大，導(dǎo)致多自由度耦合模型難以應(yīng)用.

文中在現(xiàn)有的癱船穩(wěn)性直接評估方法研究基礎(chǔ)上，根據(jù)橫搖運動時歷對船舶在隨機橫浪中的橫搖幅值進(jìn)行統(tǒng)計分析，對處于橫風(fēng)橫浪中不同破損形式船舶的橫搖運動進(jìn)行時域模擬，并運用蒙特卡洛方法計算破損船癱船狀態(tài)的傾覆概率.通過比較破損船與完整船的橫搖運動特性和傾覆概率，對船舶破損造成的穩(wěn)性損失進(jìn)行系統(tǒng)評估.

1 橫風(fēng)橫浪中時域橫搖運動模型

隨機橫風(fēng)橫浪中船舶癱船狀態(tài)的時域橫搖運動方程形式為

Mwind(t)+Mwaves(t)

(1)

1.1 非線性橫搖阻尼力矩

(2)

1.2 非線性復(fù)原力臂

船舶作大角度橫搖運動時，復(fù)原力臂曲線的非線性可能是導(dǎo)致其橫搖過大而傾覆的原因之一.此外，不同艙室破損形式會影響船舶復(fù)原力臂曲線的形狀，如艙室不對稱浸水會導(dǎo)致船舶復(fù)原力臂曲線不經(jīng)過原點.文中采用準(zhǔn)靜態(tài)模擬法計算完整船舶和破損船舶的復(fù)原力臂曲線，并采用高次多項式擬合用于后續(xù)船舶橫搖運動的模擬.

1.3 風(fēng)傾力矩

船舶受到的隨時間變化的風(fēng)傾力矩表達(dá)式為

χ(ωn))·AL·HC·cos2φ(t)·{1-sinφ(t)}

(3)

(4)

(5)

(6)

式中：Nw為規(guī)則波數(shù)量；ωn為波浪圓頻率；εn為分布在[0,2π]的隨機相位；HS為有義波高.

1.4 波浪力矩

波浪力矩考慮復(fù)原擾動力矩，其計算公式為

(7)

(8)

(9)

式中：k為波數(shù)；δω為波浪頻率間隔；HS為有義波高；TZ為平均過零周期.

2 破損船舶癱船傾覆概率計算方法

船舶破艙中最普遍的進(jìn)水情況為第三類艙室破損進(jìn)水.該類艙室破損后，艙的頂蓋在水線以上，艙內(nèi)外的海水相連通且水面保持一致.對于第三類艙室破損進(jìn)水，宜采用損失浮力法計算破損船舶的浮態(tài)以及穩(wěn)性.艙室進(jìn)水使船舶損失部分浮力，船體將會下沉以獲得補償浮力，同時產(chǎn)生一定的傾斜，當(dāng)船體下沉和傾斜后的浮心與重心重新共垂線時，船舶將處于新的平衡狀態(tài)，并認(rèn)為船舶的排水量和重心位置保持不變.

不考慮進(jìn)水過程的影響，采用四階龍格庫塔法在時域內(nèi)求解破損船舶的單自由度橫搖運動方程，并運用蒙特卡洛方法計算船舶在不同風(fēng)浪環(huán)境條件下的傾覆概率.對于船舶癱船傾覆概率的計算，可以通過大量的單自由度時域模擬得到足夠的樣本容量n，統(tǒng)計時域模擬中船舶發(fā)生傾覆的次數(shù)nc來計算船舶的傾覆概率p，即有：

(10)

當(dāng)樣本容量n足夠大時，可認(rèn)為其近似服從正態(tài)分布N(p,(1-p))，其置信區(qū)間為

(11)

IMO SDC 7-5文件的直接評估指南中建議α′取0.05，即置信度為95%，由正態(tài)分布表可以查得zα′/2的數(shù)值.為使置信度達(dá)到95%，同時誤差值不超過5%，根據(jù)統(tǒng)計學(xué)中隨機抽樣樣本容量設(shè)計方法[7]，結(jié)合后續(xù)計算所取的環(huán)境條件，文中的船舶模擬計算的樣本數(shù)取為1 000次.在數(shù)值模擬中，若船舶橫搖角超過實際傾覆角度，則認(rèn)為船舶發(fā)生傾覆.IMO將實際傾覆角度定義為船舶進(jìn)水角、穩(wěn)性消失角和50°的較小者，文中取實際傾覆角為50°，每次時域模擬的時間為1 h，若傾覆發(fā)生，則該次時域模擬結(jié)束.

3 計算算例分析

文中研究對象為DTMB5415，其主要船型參數(shù)見表1，兩種破損形式分別為第三類艙室破損的對稱浸水和不對稱浸水(僅右舷浸水)形式，破損艙信息[8]見表2.實海域中船舶航行時需安裝舭龍骨以減小船舶橫搖角，舭龍骨長47.3 m、寬0.65 m.圖1為破損船舶示意圖.

表1 DTMB5415主要船型參數(shù)

表2 DTMB5415破損艙信息 單位：m

圖1 破損船舶示意圖

3.1 CFD求解橫搖阻尼系數(shù)

采用CFD軟件STAR-CCM+模擬DTMB5415的橫搖自由衰減運動.數(shù)值模擬中縮尺比為1∶51，將船模置于靜水中，給予一定的初始橫傾角，并記錄船模橫搖角隨時間變化的曲線.

由于船模作自由衰減運動時的橫搖幅值較大，需采用重疊網(wǎng)格設(shè)置.將計算區(qū)域分為背景區(qū)域和重疊區(qū)域兩大部分.背景區(qū)域中的出口邊界類型設(shè)置為壓力出口，入口、底部、頂部以及左右兩側(cè)面均設(shè)置為速度進(jìn)口，見圖2a).重疊區(qū)域中的船體表面設(shè)置為壁面，其余邊界的類型設(shè)置為重疊網(wǎng)格.在網(wǎng)格加密方面，對自由液面進(jìn)行加密，以提高流體特征的分辨率，同時在大于重疊區(qū)域的背景區(qū)域中應(yīng)進(jìn)行圓柱型的局部加密，設(shè)置加密尺寸與重疊區(qū)域的網(wǎng)格尺寸相近，以防止重疊網(wǎng)格在隨船體橫搖發(fā)生旋轉(zhuǎn)運動時與背景網(wǎng)格之間出現(xiàn)插值錯誤，見圖2b).在重疊區(qū)域中，對于船體變化特征較明顯的區(qū)域進(jìn)行網(wǎng)格加密，以細(xì)化幾何.

圖2 CFD模擬設(shè)置

以DTMB5415完整船模(無舭龍骨)為對象，設(shè)置初始橫傾角為28°，通過改變網(wǎng)格數(shù)量(疏網(wǎng)格、中等網(wǎng)格、密網(wǎng)格數(shù)量分別為106萬、222萬、493萬)以及時間步長(分別取0.001 s和0.002 s)進(jìn)行網(wǎng)格及時間步長的敏感性分析.采用CFD進(jìn)行模擬的結(jié)果與實驗結(jié)果相比，見圖3，其橫搖幅值誤差和橫搖周期誤差均小于3%，說明CFD模擬船模的橫搖自由衰減運動具有較高的準(zhǔn)確性，且可以看出不同的網(wǎng)格數(shù)量以及時間步長對船模的橫搖運動模擬結(jié)果的影響小.綜合考慮計算資源，選取中等網(wǎng)格方案、時間步長為0.002 s作為后續(xù)計算的設(shè)置方案.

圖3 模擬結(jié)果與實驗結(jié)果對比

采用該設(shè)置方案分別模擬DTMB5415(無舭龍骨)的完整狀態(tài)、對稱浸水形式、不對稱浸水形式(初始橫傾角分別為28°、19°、0°)的橫搖自由衰減運動，得到橫搖消滅曲線見圖4.將根據(jù)橫搖消滅曲線計算得到的實船的線性橫搖阻尼系數(shù)α(1/s)、平方橫搖阻尼系數(shù)β(1/rad)和池田法計算結(jié)果進(jìn)行比較，見表3.

圖4 橫搖消滅曲線

表3 CFD和池田法計算的橫搖阻尼系數(shù)對比

由表3可知:采用CFD方法和池田法分別計算完整船舶的橫搖阻尼系數(shù)時，兩者的計算結(jié)果比較接近，說明采用池田法計算完整船的橫搖阻尼系數(shù)是可行的.但當(dāng)采用池田法計算破損船的橫搖阻尼系數(shù)時，計算結(jié)果要明顯小于CFD方法得到的結(jié)果，這是由于池田法使用船舶的船型參數(shù)計算橫搖阻尼系數(shù)，沒有計及破損艙室內(nèi)液體振蕩產(chǎn)生的水艙阻尼對船舶橫搖阻尼的影響[9].因此，目前一些文獻(xiàn)[10-11]采用池田法計算破損船舶的橫搖阻尼，其計算結(jié)果會與實際結(jié)果有較大誤差，應(yīng)當(dāng)通過模型試驗或者CFD模擬以獲得較為準(zhǔn)確結(jié)果.

實際海域中船舶安裝舭龍骨可以有效減少船舶的橫搖，采用CFD方法繼續(xù)模擬DTMB5415(帶舭龍骨)的橫搖自由衰減運動，計算得到船舶的橫搖阻尼系數(shù)和橫搖周期見表4.

表4 DTMB5415(帶舭龍骨)橫搖阻尼系數(shù)和橫搖周期

3.2 規(guī)則波中船舶橫搖運動模擬及驗證

計算得到的DTMB5415復(fù)原力臂曲線和采用“等效船舶”法計算的有效波傾系數(shù)曲線見圖5.由于第三類艙室破損后存在自由液面的影響，兩種破損船舶的最大復(fù)原力臂、穩(wěn)性消失角,以及靜穩(wěn)性曲線下的面積相比于完整狀態(tài)均減小，即船舶破損后的穩(wěn)性比完整狀態(tài)差，且不對稱浸水船舶由于有較大的初始橫傾角，其穩(wěn)性要比對稱浸水船舶的穩(wěn)性更差.此外，破損船的有效波傾系數(shù)曲線比完整船略有下降，原因是船舶破損后對波浪的繞射效應(yīng)相比完整狀態(tài)有所減弱.

圖5 復(fù)原力臂曲線和有效波傾系數(shù)曲線

編制程序模擬規(guī)則波中DTMB5415的單自由度橫搖運動，并將時域橫搖運動結(jié)果與實驗結(jié)果進(jìn)行比較，以驗證程序的可靠性.計算與實驗均不考慮風(fēng)的影響，兩者的波浪條件一致，波幅取為波長的1/100，以對稱浸水形式的DTMB5415破損船(無舭龍骨)為例，模擬破損船在規(guī)則波中的橫搖運動時歷曲線見圖6.

圖6 規(guī)則波中橫搖運動時歷曲線

各波浪條件下數(shù)值模擬與實驗[12-13]得到的橫搖角幅值結(jié)果見表5.

表5 數(shù)值計算橫搖角與實驗結(jié)果比較

由表5可知：程序計算的橫搖角幅值相對于實驗結(jié)果的誤差在可接受的范圍內(nèi)，說明程序計算具有可靠性.造成誤差的主要原因是程序采用單自由度橫搖運動模型，沒有考慮多自由度運動的耦合作用.

3.3 隨機橫浪中船舶橫搖運動模擬與分析

采用程序模擬隨機橫浪中船舶的時域橫搖運動，以對稱浸水、不對稱浸水形式破損船(均帶舭龍骨)在隨機橫浪(HS=3.5 m、9.5 m，TZ=8.5 s)中的橫搖運動為例，運動時歷見圖7～8.

圖7 隨機橫浪中破損船橫搖運動時歷(HS=3.5 m,TZ=8.5 s)

圖8 隨機橫浪中破損船橫搖運動時歷(HS=9.5 m,TZ=8.5 s)

對于線性系統(tǒng)，若其輸入為平穩(wěn)高斯隨機過程，則其輸出仍為平穩(wěn)高斯隨機過程，而船舶大角度橫搖時存在非線性，其輸出將不再符合上述規(guī)律.統(tǒng)計兩種隨機橫浪(HS分別為3.5，9.5 m，TZ均為8.5 s)時歷曲線(時長為1 h)中的一系列波幅，以及相對應(yīng)的DTMB5415(帶舭龍骨)橫搖運動時歷曲線中的一系列橫搖幅值(峰值與谷值之差)，分別繪制波幅分布和橫搖幅值分布直方圖，見圖9～10，直方圖中實折線表示瑞利分布.

圖9 波幅分布和橫搖幅值分布直方圖(HS=3.5 m,TZ=8.5 s)

圖10 波幅分布和橫搖幅值分布直方圖(HS=9.5 m,TZ=8.5 s)

如果平穩(wěn)隨機過程的瞬時值服從高斯分布，則其幅值服從瑞利分布.兩種隨機橫浪的幅值近似服從瑞利分布.當(dāng)船舶在有義波高為3.5 m、過零周期為8.5 s的隨機橫浪中作橫搖運動時，由于橫搖角度較小，此時船舶的復(fù)原力矩與橫傾角、阻尼力矩與橫搖角速度近似成線性關(guān)系，船舶系統(tǒng)仍為線性系統(tǒng)，其橫搖幅值也近似服從瑞利分布.當(dāng)船舶在有義波高為9.5 m、過零周期為8.5 s的隨機橫浪中作大角度橫搖運動時，由于船舶的復(fù)原力矩非線性和阻尼力矩非線性，此時船舶系統(tǒng)為非線性系統(tǒng)，輸出具有非線性，因此其橫搖響應(yīng)不再服從高斯分布，橫搖幅值不再服從瑞利分布，且船舶的橫搖幅值分布會整體向右偏移.

3.4 隨機橫風(fēng)橫浪中船舶傾覆概率計算

癱船狀態(tài)DTMB5415在隨機橫風(fēng)橫浪聯(lián)合作用下的部分橫搖運動時歷見圖11，若在單次模擬時間內(nèi)船舶橫搖角度超過所定義的實際傾覆角度50°，則認(rèn)為船舶已經(jīng)發(fā)生傾覆.

圖11 隨機橫風(fēng)橫浪中橫搖運動時歷(HS=9.5 m,TZ=8.5 s)

運用蒙特卡洛方法對完整狀態(tài)和兩種破損形式的DTMB5415(帶舭龍骨)在多個海況中的癱船傾覆概率進(jìn)行計算，其結(jié)果見圖12.

由圖12可知:相同海況中DTMB5415完整船的傾覆概率最小，對稱浸水形式的傾覆概率次之，不對稱浸水形式的傾覆概率最大.當(dāng)海浪的過零周期處于7.5 s附近時，此時海浪的譜峰周期處于10.56 s附近，海浪主成分波的譜峰周期接近船舶的橫搖固有周期，船舶進(jìn)入橫搖的臨界區(qū)域，產(chǎn)生十分嚴(yán)重的橫搖，三者癱船時的傾覆概率均迅速增加.而當(dāng)海浪主成分波的周期遠(yuǎn)離船舶的橫搖固有周期時，船舶發(fā)生大角度橫搖的可能性迅速減小，此時船舶的癱船傾覆概率較小.

DTMB5415癱船傾覆概率隨有義波高的變化趨勢見圖13，圖中各海況的過零周期根據(jù)北大西洋海浪譜中各有義波高對應(yīng)的最高出現(xiàn)頻率的過零周期進(jìn)行選擇.由圖13可知:當(dāng)有義波高在8.5～12.5 m，隨著有義波高的繼續(xù)增大，三者癱船時的傾覆概率均迅速增加.且在該有義波高范圍內(nèi)，相同海況中DTMB5415完整船的癱船穩(wěn)性最好，對稱浸水形式的癱船穩(wěn)性次之，不對稱浸水形式的癱船穩(wěn)性最差.

圖13 癱船狀態(tài)DTMB5415傾覆概率隨有義波高變化趨勢

4 結(jié) 論

1) IMO第二代完整穩(wěn)性中所推薦的池田法不適用于計算破損船的橫搖阻尼系數(shù)，破損船的橫搖阻尼系數(shù)應(yīng)當(dāng)通過模型實驗或者CFD模擬以獲得較為準(zhǔn)確的結(jié)果.

2) 將數(shù)值計算得到的船舶在規(guī)則波中的橫搖角幅值與實驗結(jié)果進(jìn)行對比，驗證了單自由度數(shù)值計算方法的可靠性.

3) 通過統(tǒng)計隨機橫浪中完整船和破損船的橫搖幅值，分析得出船舶在小角度橫搖時，由于船舶的復(fù)原力矩與橫搖角、阻尼力矩與橫搖角速度近似成線性關(guān)系，其橫搖幅值也近似服從瑞利分布.而當(dāng)船舶發(fā)生大角度橫搖時，由于復(fù)原力矩和阻尼力矩的非線性，船舶系統(tǒng)為非線性系統(tǒng)，其輸出具有非線性，橫搖幅值不再服從瑞利分布.

4) 相同海況中，兩種破損形式的船舶相較于完整船舶更容易發(fā)生癱船傾覆現(xiàn)象，且破損船舶不對稱浸水形式比對稱浸水形式更危險.處于癱船狀態(tài)的船舶在遭受譜峰周期與自身橫搖固有周期相近的波浪作用時，船舶進(jìn)入橫搖的臨界區(qū)域，產(chǎn)生十分嚴(yán)重的橫搖，發(fā)生傾覆的風(fēng)險顯著增加.

在計算DTMB5415癱船狀態(tài)的傾覆概率過程中，文中對有效波傾系數(shù)的計算采用SDC 7-INF.2文件中的“等效船舶”方法進(jìn)行估算，并對船舶破損區(qū)域進(jìn)行了簡化處理.在后續(xù)研究中，應(yīng)采用實驗方法確定有效波傾系數(shù)，以獲得更準(zhǔn)確的結(jié)果.同時，程序計算采用單自由度橫搖運動模型以快速獲得癱船傾覆概率計算結(jié)果，并未考慮多自由度運動之間的耦合作用，其計算結(jié)果與實際結(jié)果之間存在誤差，今后應(yīng)采用多自由度運動模型以獲得更接近船舶實際運動的結(jié)果.