隨機風浪中破損船橫搖與癱船傾覆概率計算

賴國銳 毛筱菲 詹星宇

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

0 引 言

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

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

1 橫風橫浪中時域橫搖運動模型

隨機橫風橫浪中船舶癱船狀態的時域橫搖運動方程形式為

Mwind(t)+Mwaves(t)

(1)

1.1 非線性橫搖阻尼力矩

(2)

1.2 非線性復原力臂

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

1.3 風傾力矩

船舶受到的隨時間變化的風傾力矩表達式為

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

(3)

(4)

(5)

(6)

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

1.4 波浪力矩

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

(7)

(8)

(9)

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

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

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

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

(10)

當樣本容量n足夠大時，可認為其近似服從正態分布N(p,(1-p))，其置信區間為

(11)

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

3 計算算例分析

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

表1 DTMB5415主要船型參數

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

圖1 破損船舶示意圖

3.1 CFD求解橫搖阻尼系數

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

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

圖2 CFD模擬設置

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

圖3 模擬結果與實驗結果對比

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

圖4 橫搖消滅曲線

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

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

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

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

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

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

圖5 復原力臂曲線和有效波傾系數曲線

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

圖6 規則波中橫搖運動時歷曲線

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

表5 數值計算橫搖角與實驗結果比較

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

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)

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

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

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

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

3.4 隨機橫風橫浪中船舶傾覆概率計算

癱船狀態DTMB5415在隨機橫風橫浪聯合作用下的部分橫搖運動時歷見圖11，若在單次模擬時間內船舶橫搖角度超過所定義的實際傾覆角度50°，則認為船舶已經發生傾覆.

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

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

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

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

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

4 結 論

1) IMO第二代完整穩性中所推薦的池田法不適用于計算破損船的橫搖阻尼系數，破損船的橫搖阻尼系數應當通過模型實驗或者CFD模擬以獲得較為準確的結果.

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

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

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

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