基于CFD的船舶破艙進水時域模擬

鄭 宇，馬 寧，顧解忡

(上海交通大學 船舶與海洋建筑工程學院 海洋工程國家重點實驗室 高新船舶與深海開發裝備協同創新中心，上海 200240)

基于CFD的船舶破艙進水時域模擬

鄭 宇，馬 寧，顧解忡

(上海交通大學 船舶與海洋建筑工程學院 海洋工程國家重點實驗室 高新船舶與深海開發裝備協同創新中心，上海 200240)

船舶因破損進水導致危險情況甚至傾覆是船舶安全性研究中的重要問題。船舶在破損后進水的過程中，可能存在比最終狀態更加危險的中間狀態。本文以STAR-CCM+軟件為研究工具，實現了對某客滾船瞬時非對稱進水的動態模擬。針對自由漂浮的客滾船破損情景，論文利用VOF方法模擬水氣自由液面，并首次利用重疊網格（Overset Grid）動網技術結合DFBI（Dynamic Fluid Body Interaction）六自由度求解器處理船舶破艙進水過程中的船舶升沉、橫搖等運動。模擬結果可以得到艙內液面和船舶浮態隨進水時間的變化情況，觀察到水柱射流、水花飛濺等瞬時現象；將最終浮態結果與傳統準靜態方法相比較，吻合較好。

破艙；進水；時域模擬；重疊網格

0 引 言

造成船舶海難的事故種類繁多，有碰撞、進水、擱淺、火災、沉沒等，海難給人員生命和財產帶來巨大的損失，某些情境例如油船破損還會給海洋環境帶來嚴重的破壞。因此，對船舶受損后的生存能力、浮態情況進行研究具有重要的現實意義。

船舶的破損進水過程一般分為3個階段：破損產生之后水從破口涌入船內稱為瞬時進水階段，水從內部開口浸入未破損的艙室，稱為連續進水（遞進進水）階段，最終船舶沒有傾覆或者在該階段沒有沉沒，將得到最后的穩定階段。在瞬時進水階段船舶可能達到比之后連續進水和穩定階段更大的橫傾角，造成船舶傾覆，這對進水過程的瞬態時域模擬提出了要求。對船舶破艙進水過程的時域模擬是解決船舶破損沉沒時間預報、可供救援時間預報、救援措施選擇等挑戰的有效方式。

對于船舶破艙穩性與進水的時域模擬研究已進行20余年。發生在1987年的Herald of Free Enterprise客滾船傾覆事件和1994年Estonia ferry沉沒事件強調了船舶破艙穩性的重要性，促使對于此問題的研究加速。

Spouge[1]對于European Gateway客滾船的傾覆進行了分析探討，首次提出了“對稱艙室的瞬時非對稱進水”。Santos and Guedes[2]應用了6自由度時域模擬了滾裝船的瞬時進水過程，在進水艙室中用大的流動障礙物模擬瞬時不對稱進水，整個過程假定液面水平。Zhiliang Gao等[3]用一種N-S方程求解器結合VOF方法模擬了固定狀態下客滾船進水，并用該方法測試了二維和三維的潰壩模擬問題。李佳[4]采用伯努利方程的準靜態法求解了自由漂浮狀態下不考慮外流體域的船舶進水過程。曹雪雁[5]采用SPH法對于三維船舶破艙進水特性進行了模擬研究。劉強[6]基于Fluent平臺局部網格重構動網格技術和用戶自定義函數（UDF）進行了對進水以及二維艙室耦合運動的模擬，并探討了空氣壓縮性對于進水過程的影響。上述方法均未對三維實船對象在自由漂浮狀態下的破艙流動特性、船體與水的橫搖耦合運動進行時域模擬研究。

本文對三維客滾船破艙進水過程時域模擬的方法進行探討。應用STAR-CCM+計算軟件對三維客滾船計算模型進行了破艙進水的時域模擬。采用VOF法模擬進水過程的自由液面，且首次采用重疊網格技術結合STAR-CCM+ DFBI （Dynamic Fluid Body Interaction）模塊實現對于三維客滾船實體在破損進水過程中的6-DOF運動模擬。

1 數值模擬的理論基礎

1.1 控制方程

本文的理論基礎是三維不可壓縮的粘性流體瞬態運動方程，流體的密度和粘性系數為常數。

質量守恒方程（連續性方程）：

式中：u，v，w為速度矢量v沿著x，y，z軸3個方向的速度分量。

動量守恒方程（運動方程）：

式中：F為質量力；p為壓強；μ為流體動力粘度。

有限體積法是在控制體積內對一般形式的控制微分方程的積分，即是求解積分形式的守恒方程。

式中：?為通用變量；V為控制體積；Γ為廣義擴散系數；S為廣義源項。

1.2 VOF方法

VOF（Volume of fluid）方法的基本原理是通過研究網格單元中流體和網格體積比函數F來確定自由面，追蹤流體的變化，而非追蹤自由液面上質點的運動。VOF方法可以處理自由面重入等強非線性現象，所需計算時間短、存儲量少。

如圖1所示，VOF方法中若體積比函數F=1，則說明該單元全部為指定相流體所占據；若F=0，則該單元為無指定相流體單元；當0＜F＜1時，則該單元稱為交界面單元。

圖1 VOF方法[7]Fig. 1 Volume of fluid method[7]

1.3 重疊網格（Overset Grid）動網格處理技術

Overset grid技術[8–9]是將復雜的流動區域分成幾個幾何邊界較為簡單的子區域，各個子區域之間的計算網格獨立生成，彼此存在重疊或覆蓋關系。流場信息通過插值在重疊區域的邊界上進行交換和匹配。

重疊網格邏輯關系簡單，對流場計算精度高、效率高、壁面粘性模擬能力強。重疊網格技術可以使得計算更容易執行和自動參數化：單組網格可以對應多重配置，網格質量不受位置和方向影響，邊界條件易于設置等；又有對復雜運動模擬的可能性：不必預先設定路徑，路徑可以相互交叉。

2 計算模型尺寸

本文模擬對象為某1 600客的客滾船，船體主要參數如表1所示。

建立完整船體幾何模型見圖2。計算坐標系的原點位于船體尾垂線、中心線與基線的交點，X軸向首為正，Y軸向左舷為正，Z軸向上為正。

表1 1 600 RoRo主尺度信息Tab. 1 Principal dimension information for 1 600 RoRo

圖2 完整船體幾何模型Fig. 2 Geometry model of intact ship

破損艙室為船舶左側一個全部在水線之下的機艙R6.2。示意圖見圖3，相關參數見表2。

圖3 破損艙室Fig. 3 Damaged compartment

表2 破損艙室信息Tab. 2 Information for damaged compartment

設計破口為長度與艙室長度一樣均為13.6 m（x=48 ～61.6 m），高度為0.96 m（z=4.44 ～5.4 m），位于艙室側壁上方的長方形破口。

圖4 破口Fig. 4 Opennig

3 網格劃分與計算參數設定

3.1 網格劃分

將建立的幾何面網格信息導入STAR-CCM+后，需要對其進行面的檢查、修補，例如去除自由邊、穿刺面等，進行面網格重構以滿足后續體網格生成和計算的模型質量要求。此次計算采用質量較高的切割體網格（Trimmed Cell），整個計算域即背景網格區域，如圖5所示，計算域內作橫、縱兩剖面（x=55 m和y=10 m）以展示計算區域內部的體網格情況。

圖5 整個計算域的網格配置概覽Fig. 5 Sketch of mesh arrangement in the whole computational domain

數值模擬的計算精度和計算量與網格密切相關，因此生成一套合適的網格用于計算顯得十分重要。此處根據計算區域的重要性，對于不同區域建立疏密有別的網格。背景網格基本尺寸（Base Size）為15 m，重疊網格區域網格基本尺寸為3 m。對于背景網格區域，在近重疊網格區域處、近水線面處作了網格加密；對于重疊網格區域，如圖6和圖7所示，在船體周圍、近水線面處和破艙內及周圍區域作了不同程度的網格控制加密，其中破口周圍區域的六面體網格尺寸為0.1 m。

最終生成網格數量：背景網格區域444 448，重疊網格區域2 957 687（其中破口艙室內及周圍外域加密網格數約147萬）。

圖6 重疊網格區域內部網格Fig. 6 Mesh in overset region

圖7 破口處網格細節Fig. 7 Detail of mesh around opening

圖8 為初始橫傾角=0°狀態下橫截面中重疊網格和背景網格的插值網格區域示意圖。

圖8 初始狀態插值區域Fig. 8 Interpolation region of initial condition

3.2 計算參數

計算采用隱式非定常（Implicit Unsteady）下的歐拉多相流（Eulerian Multiphase）計算模型，選用應用較多的k-epsilon湍流模型，船體在運動中視為不可變形剛體，不計空氣壓縮性。

對于邊界條件的設置，背景網格前、后邊界面設置為壁面（wall），上、下、左邊界面設置為速度入口（Velocity Inlet），右邊界面設置為壓力出口（Pressure Outlet）。具體采用STAR-CCM+中的場函數（Field Function）功能，指定速度入口的Volume Fraction函數、速度值，壓力出口的Volume Fraction函數、壓力值，控制計算中邊界處多相流分布和速度、壓力值分布。對于多相流的初始化設置，采用類似的自定義場函數方法進行。客滾船初始吃水5.323 m，無橫傾。

求解采用一階時間離散，并根據最大流動速度的預估和艙內網格尺寸要求，取隱式計算時間步長為0.005 s。

4 模擬計算與結果分析

4.1 計算結果

計算結果如圖9～圖11所示。

4.2 結果與準靜態計算結果的比較分析

NAPA軟件穩性模塊的計算精度得到了各大船級社和設計公司的肯定，將其作為驗證工具，結果可信。NAPA計算船舶浮態的方法是假定艙內液面均為水平的準靜態方法。由NAPA穩性模塊計算得到此客滾船在R6.2單艙室破損情況下最終船舶橫傾角為3.4°，吃水增加0.039 m。

如數值模擬結果圖9所示，在瞬時進水過程中，可以看到破口水柱射流、進水沖擊艙底并濺開、艙內自由液面的波動等瞬態現象；在進水后期，由于進水速度減小，液面趨于平穩，呈近似準靜態過程。如圖10和圖11所示，船舶受到的橫傾力矩在初始進水時候非線性波動很大，隨后以0為平衡位置呈阻尼衰減正弦趨勢。最大瞬時橫傾角在第2次波谷值t=25 s處出現，大小為6.12°，是最終穩態橫傾角的1.8倍。隨著橫搖時間增加，船體圍繞橫搖角–3.4°作阻尼衰減橫搖。初始2個波峰之間時間間隔為14.18 s，2個波谷之間時間間隔為13.43 s，整個計算過程橫搖周期值約為13.9 s，略大于按照杜埃爾公式估算的完整船舶固有橫搖周期12.65 s。

圖9 各時刻自由液面示意圖Fig. 9 Sketch of free surface of different time

圖10 船舶橫傾力矩時歷曲線Fig. 10 Rolling moment curve

圖11 船舶橫傾角時歷曲線Fig. 11 Rolling angle curve

5 結 語

本文以STAR-CCM+為軟件平臺，首次運用重疊網格動網技術在自由漂浮的三維客滾船實體瞬時非對稱破艙進水時域模擬中。在模擬研究中，為了提高數值運算的效率且保證計算精度，對重疊區域和背景區域進行疏密有別的網格劃分；根據破艙進水過程的物理特點，對邊界條件進行合理的設置，選擇隱式非定常求解和0.005 s的時間步長。

在與傳統準靜態結果比較中發現，在船舶不對稱進水過程中，瞬時橫傾角可能會達到最終穩態橫傾角的近兩倍。在最終未沉沒的情況下，船舶進水中后期階段以最終穩態為平衡點作阻尼衰減橫搖。本研究結果為船舶破艙進水過程的分析提供了一種新的技術路線，通過STAR-CCM+對破艙進水時流體運動情況和船舶浮態進行模擬分析，其結果可為破艙船舶的救援時間預報、救援計劃制定和新船設計方案提供數據參考依據，也為后續更復雜的船舶破艙問題研究打下基礎。

[1]SPOUGE J R. The technical investigation of the sinking of the Ro-Ro ferry EUROPEAN GATEWAY [J]. Transactions of RINA, 1985, 127: 49–72.

[2]SANTOS T A, GUEDES SOARES C. The influence of obstructions on the transient asymmetric flooding of roro ships[C]// 7thInternational Conference on stability of ships and ocean vehicles.

[3]GAO Zhi-liang, VASSALOS D, GAO Qiu-xin. Numerical simulation of water flooding into a damaged vessel's compartment by the volume of fluid method[J]. Ocean Engineering, 2010, 37: 1428–1442.

[4]李佳. 船舶破艙浸水的橫搖運動時域計算及破艙穩性研究[D]. 上海: 上海交通大學, 2009.

LI Jia. Time domain calculation of damaged ship floating and study of damage stability [D]. Shanghai: Shanghai Jiaotong University, 2009.

[5]曹雪雁, 明付仁, 張阿漫. 船舶破損進水特性的SPH三維模擬研究[C]//第十六屆中國海洋(岸)工程學術討論會論文集,2013.

CAO Xue-yan, MING Fu-ren, ZHANG A-man. The SPH 3D simulation study of ship damage and flooding character[C]//The 16thChinese Ocean Engineering Seminar Collectred Papers,2013.

[6]劉強, 段文洋. 艦船破艙進水過程的時域模擬[J]. 艦船科學技術, 2012, 3(34): 45–49.

LIU Qiang, DUAN Wen-yang. Time-domain simulation of the flooding process of damaged warships[J]. Ship Science and Technology, 2012, 3(34): 45–49.

[7]FERZIGER J H, PERIC M. Computational methods for fluid dynamics[M]. Berlin: Springer, 1996.

[8]ZHANG Z, ZHAO F, LI B. Numerical calculation of viscous free-surface flow about ship hull[J]. Journal of Ship Mechanics,2002, 6(6): 10–7.

[9]SUHS N E, ROGERS S E, DIETZ W E. PEGASUS 5: An automated pre-processor for overset-gnd CFD[R]. AIAA-2002–3186, 2002.

Time-domain simulation of the flooding of a damaged roro ship based on CFD

ZHENG Yu, MA Ning, GU Xie-chong
(State Key Laboratory of Ocean Engineering, Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiaotong University, Shanghai 200240, China)

Damage and flooding of ships will lead to danger even capsizing so that is a important problem of ship safety research. Ship may reach a more dangerous situation in mid flooding progress than final condition. In this paper, we used STAR-CCM+ as research tool to achieved the time-domain simulation of the triansient asymmetry flooding progress of a RoRo ship.Direct at free-floating damaged RoRo ship, we used VOF method to simulate the water-air free surface, and firstly used the Overset Grid dynamic mesh method with DFBI six degree of freedom solver to deal with the ship motion like heaving and rolling in ship damage and flooing problem. In the simulation result, we can get the the liquid level in compartment and ship floating conditions over time, and observe transient phenomenons like jet flow and splashing. The final floating condition result agreed well with the traditional quasi-static method result.

damage；flooding；time-domain simulation；overset

U661.3

A

1672 – 7649(2017)10 – 0029 – 05

10.3404/j.issn.1672 – 7649.2017.10.005

2016 – 12 – 23；

2017 – 02 – 13

教育部重大專項資助項目(GKZY010004)

鄭宇(1992 – )，男，碩士研究生，研究方向為船舶破艙進水數值模擬。