鄭 宇，馬 寧，顧解忡

(上海交通大學(xué) 船舶與海洋建筑工程學(xué)院 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

基于CFD的船舶破艙進(jìn)水時(shí)域模擬

鄭 宇，馬 寧，顧解忡

(上海交通大學(xué) 船舶與海洋建筑工程學(xué)院 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

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

破艙；進(jìn)水；時(shí)域模擬；重疊網(wǎng)格

0 引 言

造成船舶海難的事故種類繁多，有碰撞、進(jìn)水、擱淺、火災(zāi)、沉沒(méi)等，海難給人員生命和財(cái)產(chǎn)帶來(lái)巨大的損失，某些情境例如油船破損還會(huì)給海洋環(huán)境帶來(lái)嚴(yán)重的破壞。因此，對(duì)船舶受損后的生存能力、浮態(tài)情況進(jìn)行研究具有重要的現(xiàn)實(shí)意義。

船舶的破損進(jìn)水過(guò)程一般分為3個(gè)階段：破損產(chǎn)生之后水從破口涌入船內(nèi)稱為瞬時(shí)進(jìn)水階段，水從內(nèi)部開口浸入未破損的艙室，稱為連續(xù)進(jìn)水（遞進(jìn)進(jìn)水）階段，最終船舶沒(méi)有傾覆或者在該階段沒(méi)有沉沒(méi)，將得到最后的穩(wěn)定階段。在瞬時(shí)進(jìn)水階段船舶可能達(dá)到比之后連續(xù)進(jìn)水和穩(wěn)定階段更大的橫傾角，造成船舶傾覆，這對(duì)進(jìn)水過(guò)程的瞬態(tài)時(shí)域模擬提出了要求。對(duì)船舶破艙進(jìn)水過(guò)程的時(shí)域模擬是解決船舶破損沉沒(méi)時(shí)間預(yù)報(bào)、可供救援時(shí)間預(yù)報(bào)、救援措施選擇等挑戰(zhàn)的有效方式。

對(duì)于船舶破艙穩(wěn)性與進(jìn)水的時(shí)域模擬研究已進(jìn)行20余年。發(fā)生在1987年的Herald of Free Enterprise客滾船傾覆事件和1994年Estonia ferry沉沒(méi)事件強(qiáng)調(diào)了船舶破艙穩(wěn)性的重要性，促使對(duì)于此問(wèn)題的研究加速。

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

本文對(duì)三維客滾船破艙進(jìn)水過(guò)程時(shí)域模擬的方法進(jìn)行探討。應(yīng)用STAR-CCM+計(jì)算軟件對(duì)三維客滾船計(jì)算模型進(jìn)行了破艙進(jìn)水的時(shí)域模擬。采用VOF法模擬進(jìn)水過(guò)程的自由液面，且首次采用重疊網(wǎng)格技術(shù)結(jié)合STAR-CCM+ DFBI （Dynamic Fluid Body Interaction）模塊實(shí)現(xiàn)對(duì)于三維客滾船實(shí)體在破損進(jìn)水過(guò)程中的6-DOF運(yùn)動(dòng)模擬。

1 數(shù)值模擬的理論基礎(chǔ)

1.1 控制方程

本文的理論基礎(chǔ)是三維不可壓縮的粘性流體瞬態(tài)運(yùn)動(dòng)方程，流體的密度和粘性系數(shù)為常數(shù)。

質(zhì)量守恒方程（連續(xù)性方程）：

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

動(dòng)量守恒方程（運(yùn)動(dòng)方程）：

式中：F為質(zhì)量力；p為壓強(qiáng)；μ為流體動(dòng)力粘度。

有限體積法是在控制體積內(nèi)對(duì)一般形式的控制微分方程的積分，即是求解積分形式的守恒方程。

式中：?為通用變量；V為控制體積；Γ為廣義擴(kuò)散系數(shù)；S為廣義源項(xiàng)。

1.2 VOF方法

VOF（Volume of fluid）方法的基本原理是通過(guò)研究網(wǎng)格單元中流體和網(wǎng)格體積比函數(shù)F來(lái)確定自由面，追蹤流體的變化，而非追蹤自由液面上質(zhì)點(diǎn)的運(yùn)動(dòng)。VOF方法可以處理自由面重入等強(qiáng)非線性現(xiàn)象，所需計(jì)算時(shí)間短、存儲(chǔ)量少。

如圖1所示，VOF方法中若體積比函數(shù)F=1，則說(shuō)明該單元全部為指定相流體所占據(jù)；若F=0，則該單元為無(wú)指定相流體單元；當(dāng)0＜F＜1時(shí)，則該單元稱為交界面單元。

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

1.3 重疊網(wǎng)格（Overset Grid）動(dòng)網(wǎng)格處理技術(shù)

Overset grid技術(shù)[8–9]是將復(fù)雜的流動(dòng)區(qū)域分成幾個(gè)幾何邊界較為簡(jiǎn)單的子區(qū)域，各個(gè)子區(qū)域之間的計(jì)算網(wǎng)格獨(dú)立生成，彼此存在重疊或覆蓋關(guān)系。流場(chǎng)信息通過(guò)插值在重疊區(qū)域的邊界上進(jìn)行交換和匹配。

重疊網(wǎng)格邏輯關(guān)系簡(jiǎn)單，對(duì)流場(chǎng)計(jì)算精度高、效率高、壁面粘性模擬能力強(qiáng)。重疊網(wǎng)格技術(shù)可以使得計(jì)算更容易執(zhí)行和自動(dòng)參數(shù)化：?jiǎn)谓M網(wǎng)格可以對(duì)應(yīng)多重配置，網(wǎng)格質(zhì)量不受位置和方向影響，邊界條件易于設(shè)置等；又有對(duì)復(fù)雜運(yùn)動(dòng)模擬的可能性：不必預(yù)先設(shè)定路徑，路徑可以相互交叉。

2 計(jì)算模型尺寸

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

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

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

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

破損艙室為船舶左側(cè)一個(gè)全部在水線之下的機(jī)艙R6.2。示意圖見圖3，相關(guān)參數(shù)見表2。

圖3 破損艙室Fig. 3 Damaged compartment

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

設(shè)計(jì)破口為長(zhǎng)度與艙室長(zhǎng)度一樣均為13.6 m（x=48 ～61.6 m），高度為0.96 m（z=4.44 ～5.4 m），位于艙室側(cè)壁上方的長(zhǎng)方形破口。

圖4 破口Fig. 4 Opennig

3 網(wǎng)格劃分與計(jì)算參數(shù)設(shè)定

3.1 網(wǎng)格劃分

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

圖5 整個(gè)計(jì)算域的網(wǎng)格配置概覽Fig. 5 Sketch of mesh arrangement in the whole computational domain

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

最終生成網(wǎng)格數(shù)量：背景網(wǎng)格區(qū)域444 448，重疊網(wǎng)格區(qū)域2 957 687（其中破口艙室內(nèi)及周圍外域加密網(wǎng)格數(shù)約147萬(wàn)）。

圖6 重疊網(wǎng)格區(qū)域內(nèi)部網(wǎng)格Fig. 6 Mesh in overset region

圖7 破口處網(wǎng)格細(xì)節(jié)Fig. 7 Detail of mesh around opening

圖8 為初始橫傾角=0°狀態(tài)下橫截面中重疊網(wǎng)格和背景網(wǎng)格的插值網(wǎng)格區(qū)域示意圖。

圖8 初始狀態(tài)插值區(qū)域Fig. 8 Interpolation region of initial condition

3.2 計(jì)算參數(shù)

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

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

求解采用一階時(shí)間離散，并根據(jù)最大流動(dòng)速度的預(yù)估和艙內(nèi)網(wǎng)格尺寸要求，取隱式計(jì)算時(shí)間步長(zhǎng)為0.005 s。

4 模擬計(jì)算與結(jié)果分析

4.1 計(jì)算結(jié)果

計(jì)算結(jié)果如圖9～圖11所示。

4.2 結(jié)果與準(zhǔn)靜態(tài)計(jì)算結(jié)果的比較分析

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

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

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

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

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

5 結(jié) 語(yǔ)

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

在與傳統(tǒng)準(zhǔn)靜態(tài)結(jié)果比較中發(fā)現(xiàn)，在船舶不對(duì)稱進(jìn)水過(guò)程中，瞬時(shí)橫傾角可能會(huì)達(dá)到最終穩(wěn)態(tài)橫傾角的近兩倍。在最終未沉沒(méi)的情況下，船舶進(jìn)水中后期階段以最終穩(wěn)態(tài)為平衡點(diǎn)作阻尼衰減橫搖。本研究結(jié)果為船舶破艙進(jìn)水過(guò)程的分析提供了一種新的技術(shù)路線，通過(guò)STAR-CCM+對(duì)破艙進(jìn)水時(shí)流體運(yùn)動(dòng)情況和船舶浮態(tài)進(jìn)行模擬分析，其結(jié)果可為破艙船舶的救援時(shí)間預(yù)報(bào)、救援計(jì)劃制定和新船設(shè)計(jì)方案提供數(shù)據(jù)參考依據(jù)，也為后續(xù)更復(fù)雜的船舶破艙問(wèn)題研究打下基礎(chǔ)。

[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]李佳. 船舶破艙浸水的橫搖運(yùn)動(dòng)時(shí)域計(jì)算及破艙穩(wěn)性研究[D]. 上海: 上海交通大學(xué), 2009.

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

[5]曹雪雁, 明付仁, 張阿漫. 船舶破損進(jìn)水特性的SPH三維模擬研究[C]//第十六屆中國(guó)海洋(岸)工程學(xué)術(shù)討論會(huì)論文集,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]劉強(qiáng), 段文洋. 艦船破艙進(jìn)水過(guò)程的時(shí)域模擬[J]. 艦船科學(xué)技術(shù), 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

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

鄭宇(1992 – )，男，碩士研究生，研究方向?yàn)榇捌婆撨M(jìn)水?dāng)?shù)值模擬。