基于STAR-CCM+的船舶破艙進水模擬

王 濤， 張 維 英*， 胡 麗 芬， 于 博 文， 王 樹 祥， 魯 紀 元

( 1.大連海洋大學 航海與船舶工程學院， 遼寧 大連 116023；2.魯東大學 交通學院， 山東 煙臺 264025 )

0 引 言

海難事故多種多樣，例如著火或爆炸、進水、碰撞、擱淺、傾斜或傾覆、沉沒、失控及漂浮等[1]．海難事故導致的破艙進水會造成嚴重的經濟損失，眾多學者對船舶的破艙進水進行了研究．Spouge對一起海難事故進行了調查研究，考察了兩艘船在碰撞前后的運動、碰撞本身的機制以及隨后破損船舶快速傾斜和下沉的原因，觀察到“瞬時不對稱進水”現象，并解釋了破損船舶快速橫傾的原因，討論了這一現象的后果[2]．

近年來，計算流體力學(CFD)的興起不僅促進了實驗研究和理論分析的發展，也為流動模型的簡化提供了更多依據[3]．在船舶流體力學中，利用CFD進行船舶穩性研究已經形成趨勢，相比靜態穩性，動態穩性的研究更加復雜．Gao等觀察船舶破損過程，發現海水涌入受損艙室時，進艙水和受損船舶相互運動．為了解決這種相互運動狀態的數值模擬計算，Gao等采用基于自由表面捕獲技術的Navier-Stokes(N-S)求解模型開發了一種求解器，并驗證了求解器的準確性．最后，對滾裝渡船破損艙體進水現象進行了數值模擬，計算結果與實驗數據吻合較好[4]．盧俊尹建立了破艙進水的數理模型，通過數理模型計算驗證了數值方法的正確性，并分析了影響破艙進水的要素，結果表明其數理模型方法可以模擬實際艙室進水[5]．在研究過程中經常要考慮水、油和空氣等自由液面的變化，尤其是船舶在航行過程中處于水、氣兩層流體中，CFD軟件能夠利用多種湍流模型和有限體積法模擬這種問題．劉強等對破艙進水時域模擬進行了探討，分析了空氣流對進水過程的影響和破口水流速度變化，提出了一種減少計算域的方法[6]．

吳文鋒利用ANSYS軟件模擬了船舶碰撞過程，分析了撞擊參數對碰撞性能的影響，然后進行了物理模型實驗，利用FLUENT對雙殼油船在靜態水域碰撞后發生泄漏的情況進行三相流數值模擬，將模擬結果與物理實驗結果進行對比分析，驗證了FLUENT軟件對液貨泄漏模擬所得結果是可行的[7]．在對動態穩性的研究結果分析前，船舶在靜水中的穩性研究相比于復雜海況中的研究更有必要．靜水中的船舶穩性研究能夠給復雜工況下船舶穩性研究提供基礎．李月萌等利用FLUENT對Ruponen的模型進行破艙進水時域模擬，驗證了數值模擬方法結果的準確性，并能觀察到破艙進水每個時間步長的進水過程變化，為破艙模擬分析提供了一種新的研究方法[8]．

在船舶的破損位置、尺寸、進水量等因素確定的情況下，研究受損船舶的流體運動特性和自身動態行為可為有效減損提供參考，還可對船舶破艙穩性進行計算與評估，為船舶設計提供重要依據[9]．本文利用CFD軟件STAR-CCM+平臺，建立對一艘箱型駁船模型進行破艙進水兩相流的數值模擬模型，進行箱型駁船在靜水情況下的破艙進水時域模擬．將模擬結果與模型實驗結果進行比較，驗證STAR-CCM+對船舶破艙進水研究的可行性．然后根據實船數據利用NAPA建模，選取一個艙室，主要考慮船體和進艙水之間的橫搖耦合運動，利用建立的STAR-CCM+數值模擬方法，對實船模型進行破艙進水時域模擬，監控靜水狀態下船體和水在六自由度下的耦合運動，并對模擬結果進行分析討論．

1 數學模型

2009年馬崢等對船舶數值模擬中的湍流模型進行研究，利用FLUENT對三大主力船型即散貨船、油船和集裝箱船進行數值模擬，選取5種常用兩方程湍流模型，網格數量和劃分方法基本一致，實驗結果得到：在數值模擬選擇湍流模型時，對于集裝箱船一類，標準k-ε模型和SSTk-ω模型相對較好．對于油船、散貨船等，SSTk-ω模型具有較高的預測精度，而RNGk-ε模型對剩余阻力具有較強的預測能力[10]．

1.1 控制方程

本文以三維不可壓縮的黏性流體瞬態運動方程為理論基礎，假設流體密度和黏性系數為常數．

質量守恒方程(連續性方程)為

(1)

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

動量守恒方程(運動方程)為

(2)

式中：F為質量力，p為壓強，μ為流體動力黏度，v為速度矢量．

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

(3)

式中：φ為通用變量；V為控制體積；v為速度矢量；Γ為廣義擴散系數；S為廣義源項．

1.2 數值算法

動量守恒方程(N-S方程)為

(4)

式中：p為壓強，單位Pa；τxy、τxx、τxz等是黏性應力τ的分量，單位Pa；fx、fy、fz為x、y、z方向上的單位質量力，單位m/s2．

k-ε湍流模型方程為

Gb-ρε-Ym+Sk

(5)

式中：Gk為速度梯度產生的湍動能項；Gb為浮力產生的湍動能項；i，j=1，2，3分別表示x、y、z3個方向；Ym為脈動擴張項；C1ε、C2ε、C3ε為經驗常數；σk、σε為與湍動能k和耗散率ε相對應的Prandtl數；Sk、Sε為用戶自定義的源項．

采用基于有限元分析的STAR-CCM+軟件求解模型內流場．不可壓縮黏性流體的控制方程由雷諾平均Navier-Stokes方程和連續性方程描述．為了封閉這組方程，采用了k-ε湍流模型，對于管內的充分發展湍流，無論是穩態還是非穩態的流固耦合問題，k-ε湍流模型都是很適用的[11]．

2 STAR-CCM+模擬要點

2.1 數值模擬要點

破損船舶進水過程時域模型的建立基于兩個條件：一是艙室流量平衡，即與艙室相連的所有破口的流量之和與艙室內進水增量相等；二是破口處滿足伯努利方程[12]．

計算域模型的網格劃分方法和網格質量對數值求解的計算精度及模擬結果都具有非常大的影響，采用合理的網格控制參數和局部區域(關鍵主流區域)網格細化控制方法進行網格劃分，對減少網格數量、提高計算精度和求解效率具有非常重要的作用[13]．

2.2 邊界條件設置

各個艙室均有通風口，不考慮封閉空間氣體壓強因素，劃分區域為船體(HULL)、進水艙室(HULL INSIDE)、隨體區域(OVERSET ZONE)、變形區域(TANK)．設置氣液交界面(WAVE)、破口(INLET A)，標準模型頂部設置壓力出口(INLET B)、混合壁面(INSIDE WALL)，混合壁面切向速度固定，剪應力無滑移，壁面規整平滑．

3 標準模型破艙模擬研究

3.1 標準模型數據與模型尺寸

本文首先采用Ruponen的實驗模型[14]．模型來源于阿爾托大學和NAPA：ITTC基準研究的洪水模型測試(ITTC-Box)，收錄于船舶與海洋運載器國際穩性會議(STAB)文集，STAB是全球穩定性和一般安全領域最具代表性的專業會議．模型主要參數見表1，箱型駁船標準模型如圖1所示，箱型駁船外形尺寸如圖2所示．

表1 標準模型主要參數

圖1 箱型駁船

圖2 箱型駁船外形尺寸

3.2 收斂性分析

STAR-CCM+支持對幾何模型進行網格劃分，網格包括四面體網格、多面體網格等，與FLUENT相比，支持網格模型較多．此外，STAR-CCM+還支持從其他第三方網格軟件(ICEM CFD、Pointwise、Gridpro、Hypermesh等)直接導入體網格進行計算．

破艙進水模擬屬于模擬外部流，因此選擇切割體網格．切割體網格單元生成器能提供穩定且有效的方法，針對簡單和復雜的網格生成問題生成高質量的網格，同時適用于基于零部件的網格化(PBM)和基于區域的網格化(RBM)．除了模型本身質量的影響，面網格質量最小值和時間步長的設置也大大影響了計算精度和計算時間．

為了得到不同網格和不同時間步長對模擬結果的影響，選取箱型駁船進水艙室(HULL INSIDE)進行收斂性分析．

進水艙室尺寸如圖3所示，不同網格比較情況如表2所示．由表2可以看出，面網格質量最小值為0.05時，劃分網格質量最好，網格數量也滿足計算要求．

圖3 進水艙室外形尺寸

表2 不同網格對比表

庫朗數能夠將STAR-CCM+中的網格和時間步長聯系起來．其公式為

C=vΔt/Δx

(6)

式中：v為流體流速；Δt為用于VOF方法計算的時間步長；Δx為網格間距；庫朗數C<0.3．殘差收斂較好的可以適當放大庫朗數．

面網格質量最小值均為0.05，采用不同時間步長，模擬終止條件為30 s，則不同時間步長的模擬時間如表3所示．

表3 不同時間步長對比表

殘差分析是考察模型假設合理性的方法，通過對殘差及殘差圖的分析，可以判斷選取的時間步長是否滿足計算需要．

當時間步長設置為0.30 s時，湍動能隨著計算的增加呈上升趨勢，收斂較差．時間步長設置為0.02、0.05 s時，湍流耗散率和湍動能殘差基本相同，都較為平穩，略有向下趨勢，因此判斷收斂較好，但時間步長設置為0.02 s時，模擬時間較長．結合表3和圖4可以看出，時間步長設置為0.05 s，既能保證計算精度，同時模擬時間較少．

3.3 箱型駁船網格建立

由于破艙進水模擬中存在自由液面，為保證計算精度，本文利用STAR-CCM+切割體網格生成器對實船模型進行網格劃分，對隨體區域和變形區域建立重疊網格．進水艙室及破口尺寸見表4．

(a) 時間步長0.30 s

(b) 時間步長0.05 s

(c) 時間步長0.02 s

圖4 殘差分析

Fig.4 Residuals analysis

表4 進水艙室及破口尺寸

選取表面重構，執行曲率細化和接近值細化．自動表面修復設置面網格質量最小值為0.05；選擇切割體網格單元生成器和棱柱層網格生成器，延伸函數為幾何級數，填隙百分比為25%，最小厚度百分比為10%，層減少百分比為50%．面網格增長率1.3，自動修復最小接近值0.01．網格總數量為2 526 717，其中質量較好的網格數量為2 507 240，占99.229%，如圖5所示．進水艙室網格劃分，設置表面重構，選取切割體網格單元生成器和棱柱層網格生成器．共生成網格684 196，其中質量較好的網格數量為675 567，占98.739%，如圖6所示．變形區域網格劃分為877 160，其中質量較好的網格數量為876 960，占99.977%，如圖7所示．隨體區域網格劃分為965 361，其中質量較好的網格數量為954 713，占98.897%，如圖8所示．

選取常用k-ε湍流模型六自由度求解器，由于不考慮波浪在靜水狀態下進行模擬，所以VOF波設置為靜水VOF波．實際環境中，船底部為液體，上部為氣體，所以為氣液兩相模擬．選取歐拉多相流模型，水設置為恒密度，空氣部分也設置為恒密度，動力黏度和密度均設置為常數．在STAR-CCM+中設置物理模型，模型參數選擇初始條件水和氣體復合湍流強度0.01，湍流速度1 m/s．設置求解器參數：隱式不定常中時間步長設置為0.01 s，默認設置負載平衡、k-ε湍流、分離VOF模型等選項．

圖5 整體網格

圖6 艙室網格

圖7 船體外側變形區域網格劃分

圖8 隨體區域網格加密

3.4 標準模型進水過程及結果分析

利用STAR-CCM+模擬結果與實驗結果對比，結果高度吻合，數據基本準確．箱型駁船進水過程如圖9所示．0.40 s時，可以觀察到破口進水開始，進艙水在破口處形成水柱沖擊艙壁．6.60 s時，破口艙室的進水已經流入中間艙室，中部艙室進水，部分水流入右側艙室．10.50 s時，3個艙室均有進艙水，艙室內液面逐漸升高．29.95 s時，3個艙室達到水滿狀態．

監視破口進水艙室內的自由液面高度，繪制液面高度隨時間變化曲線，并與文獻[8，14]結果進行比較，如圖10所示．

(a) t=0.40 s

(b) t=6.60 s

(c) t=10.50 s

(d) t=29.95 s

圖9 箱型駁船模擬結果

Fig.9 Simulation results of box barge

本文STAR-CCM+模擬結果與Ruponen的實驗結果[14]、FLUENT模擬結果對比[8]，R21S艙在t=0～12.00 s時，由于模擬開始時部分進艙水噴射在R21S艙壁上產生濺射，液面高度較實驗結果與文獻結果略高；t=12.00 s以后，模擬結果逐漸趨于平穩，與實驗結果、文獻結果高度吻合．由于左側R21S艙一部分進艙水通過門直接噴射到R21艙室中，使監測的R21艙平均液面高度增加，但隨著艙室內進水量的增加，液面升高，誤差逐漸減小，最后監測結果逐漸穩定，與文獻結果相同．R21P艙進水無水柱噴射，進水狀態平穩，艙室液面高度變化與實驗結果、文獻結果對比高度吻合．

(a) R21S

(b) R21

(c) R21P

圖10 艙室液面高度隨時間變化圖

Fig.10 Diagram of cabin level change with time

模擬對比結果驗證了該模擬方法可行，利用STAR-CCM+可以模擬船舶破艙進水過程，模擬結果準確可靠且模擬方法較FLUENT模擬更方便，占用系統資源更少．

4 漁政船破艙進水模擬研究

4.1 NAPA建模

本文采用一艘522 kW漁政船的數據，利用NAPA進行船舶建模．漁政船是指在漁業專屬水域執行漁政任務，擔負海上漁政管理監督、執法的專業船只，也是漁政管理設施建設中重要的設施之一；主要用于漁場巡視并監督、檢查漁船執行國家漁業法規的情況，也可兼負漁業生產指揮、發布漁情和氣象通報以及海上醫療、海難救助等任務[15]．因此漁政船所航行海域復雜多變，容易發生觸礁、剮蹭等事故．機艙破口進水將嚴重影響船舶航行性能．本文選取機艙進行破艙進水的時域模擬．

船體主要參數如表5所示．

表5 漁政船主要參數

建立漁政船模型，如圖11所示，艙室劃分如圖12所示．

選取機艙修改幾何，在右側做半徑0.30 m的圓形破口，艙室頂部做通風口，如圖13所示．由于模擬處于靜水環境中，漁政船無動力、無航速，水流速度為0 m/s，忽略進艙水晃蕩的影響，因此在模擬過程中，不考慮艙室內設備影響，視機艙為空艙．

圖11 522 kW漁政船模型

圖12 艙室結構劃分示意圖

圖13 機艙幾何模型

4.2 STAR-CCM+網格劃分

變形區域選取切割體網格單元生成器，打開表面重構和自動表面生成器，面網格增長率1.3，自動修復最小接近值0.01；進水艙室選取切割體網格單元生成器、棱柱層網格生成器，打開表面重構和自動表面生成器，面網格增長率1.3，自動修復最小接近值0.01，棱柱層數1，棱柱層延伸1.5，棱柱層總厚度絕對尺寸0.005 m．考慮運算時間，在重疊域之間建立重疊網格交界面．求解器設置為隱式不定常、六自由度求解器，六自由度運動、負載平衡、分離流、分離VOF波、k-ε湍流、k-ε湍流黏度．

總網格數量為699 301，檢查網格質量，質量較好的網格達到98.832%．

機艙網格劃分總網格數量為707 564，其中質量較好的為699 301，占98.832%，如圖14所示；破口處網格加密如圖15所示．

圖14 機艙網格劃分

圖15 機艙破口網格細節

船體網格劃分如圖16所示，劃分總網格數量為373 556，其中質量較好的為350 977，占93.956%．

圖16 船體網格劃分

變形區域網格劃分如圖17所示，劃分總數量為795 358，其中質量較好的網格數量為795 166，占99.976%．

船艙尺寸9 m×6 m×2.8 m，艙壁厚度0.1 m，船艙吃水2.15 m，破口為半徑0.30 m的圓形口，外域尺寸150 m×100 m×40 m，重疊網格區域基本尺寸為56 m．設定0.05 s為時間步長，對于邊界條件的設置，頂部邊界設置為TOP，底部邊界設置為BOTTOM．

圖17 變形區域網格加密

DFBI設置六自由度體，啟用模型VOF波，考慮重力因素，兩層全y+壁面處理，精確壁面距離．選擇可實現的k-ε兩層模型，常用k-ε湍流模型，雷諾平均Navier-Stokes模型，三維隱式不定常模型．選取歐拉多相流模型，水設置為恒密度，空氣部分也設置為恒密度，動力黏度和密度均設置為常數．右邊界破口設置為速度入口(INLET A)，頂部通風口設置為壓力出口(INLET B)．監視器監視模型六自由度的變化，并生成變化曲線．具體采用STAR-CCM+中的混合壁面函數功能，指定速度入口的函數值、壓力出口的函數值，控制計算中邊界處多相流分布和速度、壓力分布等．

4.3 結果分析討論

與標準模型實驗結果進行對比，已經證明STAR-CCM+在處理船舶穩性研究中可以得到較為可信的結果．采用實船數據利用NAPA進行建模得到實船模型，在破艙進水的過程中，不考慮進艙水的晃蕩影響．進水過程瞬時變化如圖18所示．

模擬過程中可以觀察到在瞬時過程的進水前期，t=8.20 s時，破口處形成水柱噴射，沖擊在艙壁上，機艙內產生自由液面，并在船艙底部形成兩處渦流，艙內液面逐漸升高；t=80.00 s時，艙內液面高度超過破口高度，艙內進水不形成水柱噴射，但進艙水仍然有渦流．

監測破口處流量隨時間變化，生成破口流量變化曲線，如圖19所示．可以看出在t=0.07 s時，破口流量為655.81 kg/s．在t=0.07～55.00 s時，由于船體橫搖影響，破口流量也產生規則波動，隨著進艙水的增加，橫搖逐漸趨于穩定，破口流量也逐漸趨于平穩．t=100.00 s時升沉運動持續增大，破口到水面垂直距離增加，破口流量也隨之增加．對破口流量曲線進行積分，使用科研繪圖工具GraphPad Prism求出曲線下面積，機艙浸滿水時，艙內水有94 850 kg．

(a) 船體t=1.00 s

(b) 機艙t=1.00 s

(c) 船體t=8.20 s

(d) 機艙t=8.20 s

(e) t=80.00 s

(f) t=100.00 s

(g) t=125.00 s

(h) t=155.00 s

圖18 破艙進水過程圖

Fig.18 Water intake process of damaged cabin

圖19 破口流量圖

船體六自由度變化如圖20～22所示，船身向進水一側移動，橫蕩距離增加，船體形成較大的橫傾角．

第1次橫傾角最大值出現在t=0.66 s，偏移角度為1.29°，是艙室浸滿水時的瞬時橫傾角的16倍，橫搖周期約為5.080 s，比設計計算的橫搖周期5.199 s略小．

隨著進水時間增加，橫傾角逐漸減小，艙室內液面升高，液面趨于平穩，機艙進水增加，船舶下沉，艏搖加快，機艙進水狀態趨于平穩．由于進艙水影響，縱搖逐漸減小，艙室浸滿水約160 s．

(a) 橫搖

(b) 橫蕩

(a) 縱搖

(b) 縱蕩

(a) 艏搖

(b) 升沉

5 結 語

本文利用CFD軟件STAR-CCM+平臺，運用動網格和重疊網格技術對靜水狀態下的實船模型進行破艙進水模擬．為了保證模擬方法的準確性，先對網格和時間步長收斂性進行分析，然后與標準箱型駁船實驗結果進行對比，驗證了STAR-CCM+應用于模擬船舶破艙進水的可行性，模擬結果準確可靠．STAR-CCM+的六自由度求解器、VOF波形模型、可實現的k-ε兩層模型、歐拉多相流模型均適用于船舶實驗的模擬．通過模擬實驗，可以觀察船艙進水的瞬時狀態，通過對監視器的六自由度變化曲線進行分析，能夠得到船舶破損時的穩態瞬時數據，依靠軟件強大的后處理功能，可將模擬結果生成圖表或視頻，轉為可視化數據．

在確定了一艘漁政船破損位置、尺寸等因素的情況下，對漁政船進行破艙進水時域模擬得到漁政船每一時間步長破損艙室進水位置云圖和破艙進水破口流量及六自由度隨時間變化曲線．剛開始進水時船舶產生較大橫搖角，破口流量規則波動，隨著進水量的增加船體下沉，橫搖減弱并逐漸穩定，破口流量增大，縱搖隨進水量增加逐漸減弱，模擬在艙室浸滿水后停止．

本文實驗結果可為研究復雜工況如多個艙室破損進水、大風大浪環境下的船舶破艙進水模擬提供基礎數據，并為研究船舶破艙進水提供一種新思路和方法．其分析結果可為船體優化、海難船舶救援提供參考．