王 濤, 張 維 英*, 胡 麗 芬, 于 博 文, 王 樹 祥, 魯 紀 元
( 1.大連海洋大學 航海與船舶工程學院, 遼寧 大連 116023;2.魯東大學 交通學院, 山東 煙臺 264025 )
海難事故多種多樣,例如著火或爆炸、進水、碰撞、擱淺、傾斜或傾覆、沉沒、失控及漂浮等[1].海難事故導致的破艙進水會造成嚴重的經(jīng)濟損失,眾多學者對船舶的破艙進水進行了研究.Spouge對一起海難事故進行了調(diào)查研究,考察了兩艘船在碰撞前后的運動、碰撞本身的機制以及隨后破損船舶快速傾斜和下沉的原因,觀察到“瞬時不對稱進水”現(xiàn)象,并解釋了破損船舶快速橫傾的原因,討論了這一現(xiàn)象的后果[2].
近年來,計算流體力學(CFD)的興起不僅促進了實驗研究和理論分析的發(fā)展,也為流動模型的簡化提供了更多依據(jù)[3].在船舶流體力學中,利用CFD進行船舶穩(wěn)性研究已經(jīng)形成趨勢,相比靜態(tài)穩(wěn)性,動態(tài)穩(wěn)性的研究更加復雜.Gao等觀察船舶破損過程,發(fā)現(xiàn)海水涌入受損艙室時,進艙水和受損船舶相互運動.為了解決這種相互運動狀態(tài)的數(shù)值模擬計算,Gao等采用基于自由表面捕獲技術的Navier-Stokes(N-S)求解模型開發(fā)了一種求解器,并驗證了求解器的準確性.最后,對滾裝渡船破損艙體進水現(xiàn)象進行了數(shù)值模擬,計算結(jié)果與實驗數(shù)據(jù)吻合較好[4].盧俊尹建立了破艙進水的數(shù)理模型,通過數(shù)理模型計算驗證了數(shù)值方法的正確性,并分析了影響破艙進水的要素,結(jié)果表明其數(shù)理模型方法可以模擬實際艙室進水[5].在研究過程中經(jīng)常要考慮水、油和空氣等自由液面的變化,尤其是船舶在航行過程中處于水、氣兩層流體中,CFD軟件能夠利用多種湍流模型和有限體積法模擬這種問題.劉強等對破艙進水時域模擬進行了探討,分析了空氣流對進水過程的影響和破口水流速度變化,提出了一種減少計算域的方法[6].
吳文鋒利用ANSYS軟件模擬了船舶碰撞過程,分析了撞擊參數(shù)對碰撞性能的影響,然后進行了物理模型實驗,利用FLUENT對雙殼油船在靜態(tài)水域碰撞后發(fā)生泄漏的情況進行三相流數(shù)值模擬,將模擬結(jié)果與物理實驗結(jié)果進行對比分析,驗證了FLUENT軟件對液貨泄漏模擬所得結(jié)果是可行的[7].在對動態(tài)穩(wěn)性的研究結(jié)果分析前,船舶在靜水中的穩(wěn)性研究相比于復雜海況中的研究更有必要.靜水中的船舶穩(wěn)性研究能夠給復雜工況下船舶穩(wěn)性研究提供基礎.李月萌等利用FLUENT對Ruponen的模型進行破艙進水時域模擬,驗證了數(shù)值模擬方法結(jié)果的準確性,并能觀察到破艙進水每個時間步長的進水過程變化,為破艙模擬分析提供了一種新的研究方法[8].
在船舶的破損位置、尺寸、進水量等因素確定的情況下,研究受損船舶的流體運動特性和自身動態(tài)行為可為有效減損提供參考,還可對船舶破艙穩(wěn)性進行計算與評估,為船舶設計提供重要依據(jù)[9].本文利用CFD軟件STAR-CCM+平臺,建立對一艘箱型駁船模型進行破艙進水兩相流的數(shù)值模擬模型,進行箱型駁船在靜水情況下的破艙進水時域模擬.將模擬結(jié)果與模型實驗結(jié)果進行比較,驗證STAR-CCM+對船舶破艙進水研究的可行性.然后根據(jù)實船數(shù)據(jù)利用NAPA建模,選取一個艙室,主要考慮船體和進艙水之間的橫搖耦合運動,利用建立的STAR-CCM+數(shù)值模擬方法,對實船模型進行破艙進水時域模擬,監(jiān)控靜水狀態(tài)下船體和水在六自由度下的耦合運動,并對模擬結(jié)果進行分析討論.
2009年馬崢等對船舶數(shù)值模擬中的湍流模型進行研究,利用FLUENT對三大主力船型即散貨船、油船和集裝箱船進行數(shù)值模擬,選取5種常用兩方程湍流模型,網(wǎng)格數(shù)量和劃分方法基本一致,實驗結(jié)果得到:在數(shù)值模擬選擇湍流模型時,對于集裝箱船一類,標準k-ε模型和SSTk-ω模型相對較好.對于油船、散貨船等,SSTk-ω模型具有較高的預測精度,而RNGk-ε模型對剩余阻力具有較強的預測能力[10].
本文以三維不可壓縮的黏性流體瞬態(tài)運動方程為理論基礎,假設流體密度和黏性系數(shù)為常數(shù).
質(zhì)量守恒方程(連續(xù)性方程)為
(1)
式中:u、v、w為速度矢量v沿著x、y、z軸3個方向的速度分量.
動量守恒方程(運動方程)為
(2)
式中:F為質(zhì)量力,p為壓強,μ為流體動力黏度,v為速度矢量.
有限體積法是在控制體積內(nèi)對一般形式的控制微分方程的積分,即是求解積分形式的守恒方程.
(3)
式中:φ為通用變量;V為控制體積;v為速度矢量;Γ為廣義擴散系數(shù);S為廣義源項.
動量守恒方程(N-S方程)為
(4)
式中:p為壓強,單位Pa;τxy、τxx、τxz等是黏性應力τ的分量,單位Pa;fx、fy、fz為x、y、z方向上的單位質(zhì)量力,單位m/s2.
k-ε湍流模型方程為
Gb-ρε-Ym+Sk
(5)
式中:Gk為速度梯度產(chǎn)生的湍動能項;Gb為浮力產(chǎn)生的湍動能項;i,j=1,2,3分別表示x、y、z3個方向;Ym為脈動擴張項;C1ε、C2ε、C3ε為經(jīng)驗常數(shù);σk、σε為與湍動能k和耗散率ε相對應的Prandtl數(shù);Sk、Sε為用戶自定義的源項.
采用基于有限元分析的STAR-CCM+軟件求解模型內(nèi)流場.不可壓縮黏性流體的控制方程由雷諾平均Navier-Stokes方程和連續(xù)性方程描述.為了封閉這組方程,采用了k-ε湍流模型,對于管內(nèi)的充分發(fā)展湍流,無論是穩(wěn)態(tài)還是非穩(wěn)態(tài)的流固耦合問題,k-ε湍流模型都是很適用的[11].
破損船舶進水過程時域模型的建立基于兩個條件:一是艙室流量平衡,即與艙室相連的所有破口的流量之和與艙室內(nèi)進水增量相等;二是破口處滿足伯努利方程[12].
計算域模型的網(wǎng)格劃分方法和網(wǎng)格質(zhì)量對數(shù)值求解的計算精度及模擬結(jié)果都具有非常大的影響,采用合理的網(wǎng)格控制參數(shù)和局部區(qū)域(關鍵主流區(qū)域)網(wǎng)格細化控制方法進行網(wǎng)格劃分,對減少網(wǎng)格數(shù)量、提高計算精度和求解效率具有非常重要的作用[13].
各個艙室均有通風口,不考慮封閉空間氣體壓強因素,劃分區(qū)域為船體(HULL)、進水艙室(HULL INSIDE)、隨體區(qū)域(OVERSET ZONE)、變形區(qū)域(TANK).設置氣液交界面(WAVE)、破口(INLET A),標準模型頂部設置壓力出口(INLET B)、混合壁面(INSIDE WALL),混合壁面切向速度固定,剪應力無滑移,壁面規(guī)整平滑.
本文首先采用Ruponen的實驗模型[14].模型來源于阿爾托大學和NAPA:ITTC基準研究的洪水模型測試(ITTC-Box),收錄于船舶與海洋運載器國際穩(wěn)性會議(STAB)文集,STAB是全球穩(wěn)定性和一般安全領域最具代表性的專業(yè)會議.模型主要參數(shù)見表1,箱型駁船標準模型如圖1所示,箱型駁船外形尺寸如圖2所示.
表1 標準模型主要參數(shù)
圖1 箱型駁船
圖2 箱型駁船外形尺寸
STAR-CCM+支持對幾何模型進行網(wǎng)格劃分,網(wǎng)格包括四面體網(wǎng)格、多面體網(wǎng)格等,與FLUENT相比,支持網(wǎng)格模型較多.此外,STAR-CCM+還支持從其他第三方網(wǎng)格軟件(ICEM CFD、Pointwise、Gridpro、Hypermesh等)直接導入體網(wǎng)格進行計算.
破艙進水模擬屬于模擬外部流,因此選擇切割體網(wǎng)格.切割體網(wǎng)格單元生成器能提供穩(wěn)定且有效的方法,針對簡單和復雜的網(wǎng)格生成問題生成高質(zhì)量的網(wǎng)格,同時適用于基于零部件的網(wǎng)格化(PBM)和基于區(qū)域的網(wǎng)格化(RBM).除了模型本身質(zhì)量的影響,面網(wǎng)格質(zhì)量最小值和時間步長的設置也大大影響了計算精度和計算時間.
為了得到不同網(wǎng)格和不同時間步長對模擬結(jié)果的影響,選取箱型駁船進水艙室(HULL INSIDE)進行收斂性分析.
進水艙室尺寸如圖3所示,不同網(wǎng)格比較情況如表2所示.由表2可以看出,面網(wǎng)格質(zhì)量最小值為0.05時,劃分網(wǎng)格質(zhì)量最好,網(wǎng)格數(shù)量也滿足計算要求.
圖3 進水艙室外形尺寸
表2 不同網(wǎng)格對比表
庫朗數(shù)能夠?qū)TAR-CCM+中的網(wǎng)格和時間步長聯(lián)系起來.其公式為
C=vΔt/Δx
(6)
式中:v為流體流速;Δt為用于VOF方法計算的時間步長;Δx為網(wǎng)格間距;庫朗數(shù)C<0.3.殘差收斂較好的可以適當放大庫朗數(shù).
面網(wǎng)格質(zhì)量最小值均為0.05,采用不同時間步長,模擬終止條件為30 s,則不同時間步長的模擬時間如表3所示.
表3 不同時間步長對比表
殘差分析是考察模型假設合理性的方法,通過對殘差及殘差圖的分析,可以判斷選取的時間步長是否滿足計算需要.
當時間步長設置為0.30 s時,湍動能隨著計算的增加呈上升趨勢,收斂較差.時間步長設置為0.02、0.05 s時,湍流耗散率和湍動能殘差基本相同,都較為平穩(wěn),略有向下趨勢,因此判斷收斂較好,但時間步長設置為0.02 s時,模擬時間較長.結(jié)合表3和圖4可以看出,時間步長設置為0.05 s,既能保證計算精度,同時模擬時間較少.
由于破艙進水模擬中存在自由液面,為保證計算精度,本文利用STAR-CCM+切割體網(wǎng)格生成器對實船模型進行網(wǎng)格劃分,對隨體區(qū)域和變形區(qū)域建立重疊網(wǎng)格.進水艙室及破口尺寸見表4.
(a) 時間步長0.30 s
(b) 時間步長0.05 s
(c) 時間步長0.02 s
圖4 殘差分析
Fig.4 Residuals analysis
表4 進水艙室及破口尺寸
選取表面重構(gòu),執(zhí)行曲率細化和接近值細化.自動表面修復設置面網(wǎng)格質(zhì)量最小值為0.05;選擇切割體網(wǎng)格單元生成器和棱柱層網(wǎng)格生成器,延伸函數(shù)為幾何級數(shù),填隙百分比為25%,最小厚度百分比為10%,層減少百分比為50%.面網(wǎng)格增長率1.3,自動修復最小接近值0.01.網(wǎng)格總數(shù)量為2 526 717,其中質(zhì)量較好的網(wǎng)格數(shù)量為2 507 240,占99.229%,如圖5所示.進水艙室網(wǎng)格劃分,設置表面重構(gòu),選取切割體網(wǎng)格單元生成器和棱柱層網(wǎng)格生成器.共生成網(wǎng)格684 196,其中質(zhì)量較好的網(wǎng)格數(shù)量為675 567,占98.739%,如圖6所示.變形區(qū)域網(wǎng)格劃分為877 160,其中質(zhì)量較好的網(wǎng)格數(shù)量為876 960,占99.977%,如圖7所示.隨體區(qū)域網(wǎng)格劃分為965 361,其中質(zhì)量較好的網(wǎng)格數(shù)量為954 713,占98.897%,如圖8所示.
選取常用k-ε湍流模型六自由度求解器,由于不考慮波浪在靜水狀態(tài)下進行模擬,所以VOF波設置為靜水VOF波.實際環(huán)境中,船底部為液體,上部為氣體,所以為氣液兩相模擬.選取歐拉多相流模型,水設置為恒密度,空氣部分也設置為恒密度,動力黏度和密度均設置為常數(shù).在STAR-CCM+中設置物理模型,模型參數(shù)選擇初始條件水和氣體復合湍流強度0.01,湍流速度1 m/s.設置求解器參數(shù):隱式不定常中時間步長設置為0.01 s,默認設置負載平衡、k-ε湍流、分離VOF模型等選項.
圖5 整體網(wǎng)格
圖6 艙室網(wǎng)格
圖7 船體外側(cè)變形區(qū)域網(wǎng)格劃分
圖8 隨體區(qū)域網(wǎng)格加密
利用STAR-CCM+模擬結(jié)果與實驗結(jié)果對比,結(jié)果高度吻合,數(shù)據(jù)基本準確.箱型駁船進水過程如圖9所示.0.40 s時,可以觀察到破口進水開始,進艙水在破口處形成水柱沖擊艙壁.6.60 s時,破口艙室的進水已經(jīng)流入中間艙室,中部艙室進水,部分水流入右側(cè)艙室.10.50 s時,3個艙室均有進艙水,艙室內(nèi)液面逐漸升高.29.95 s時,3個艙室達到水滿狀態(tài).
監(jiān)視破口進水艙室內(nèi)的自由液面高度,繪制液面高度隨時間變化曲線,并與文獻[8,14]結(jié)果進行比較,如圖10所示.
(a) t=0.40 s
(b) t=6.60 s
(c) t=10.50 s
(d) t=29.95 s
圖9 箱型駁船模擬結(jié)果
Fig.9 Simulation results of box barge
本文STAR-CCM+模擬結(jié)果與Ruponen的實驗結(jié)果[14]、FLUENT模擬結(jié)果對比[8],R21S艙在t=0~12.00 s時,由于模擬開始時部分進艙水噴射在R21S艙壁上產(chǎn)生濺射,液面高度較實驗結(jié)果與文獻結(jié)果略高;t=12.00 s以后,模擬結(jié)果逐漸趨于平穩(wěn),與實驗結(jié)果、文獻結(jié)果高度吻合.由于左側(cè)R21S艙一部分進艙水通過門直接噴射到R21艙室中,使監(jiān)測的R21艙平均液面高度增加,但隨著艙室內(nèi)進水量的增加,液面升高,誤差逐漸減小,最后監(jiān)測結(jié)果逐漸穩(wěn)定,與文獻結(jié)果相同.R21P艙進水無水柱噴射,進水狀態(tài)平穩(wěn),艙室液面高度變化與實驗結(jié)果、文獻結(jié)果對比高度吻合.
(a) R21S
(b) R21
(c) R21P
圖10 艙室液面高度隨時間變化圖
Fig.10 Diagram of cabin level change with time
模擬對比結(jié)果驗證了該模擬方法可行,利用STAR-CCM+可以模擬船舶破艙進水過程,模擬結(jié)果準確可靠且模擬方法較FLUENT模擬更方便,占用系統(tǒng)資源更少.
本文采用一艘522 kW漁政船的數(shù)據(jù),利用NAPA進行船舶建模.漁政船是指在漁業(yè)專屬水域執(zhí)行漁政任務,擔負海上漁政管理監(jiān)督、執(zhí)法的專業(yè)船只,也是漁政管理設施建設中重要的設施之一;主要用于漁場巡視并監(jiān)督、檢查漁船執(zhí)行國家漁業(yè)法規(guī)的情況,也可兼負漁業(yè)生產(chǎn)指揮、發(fā)布漁情和氣象通報以及海上醫(yī)療、海難救助等任務[15].因此漁政船所航行海域復雜多變,容易發(fā)生觸礁、剮蹭等事故.機艙破口進水將嚴重影響船舶航行性能.本文選取機艙進行破艙進水的時域模擬.
船體主要參數(shù)如表5所示.
表5 漁政船主要參數(shù)
建立漁政船模型,如圖11所示,艙室劃分如圖12所示.
選取機艙修改幾何,在右側(cè)做半徑0.30 m的圓形破口,艙室頂部做通風口,如圖13所示.由于模擬處于靜水環(huán)境中,漁政船無動力、無航速,水流速度為0 m/s,忽略進艙水晃蕩的影響,因此在模擬過程中,不考慮艙室內(nèi)設備影響,視機艙為空艙.
圖11 522 kW漁政船模型
圖12 艙室結(jié)構(gòu)劃分示意圖
圖13 機艙幾何模型
變形區(qū)域選取切割體網(wǎng)格單元生成器,打開表面重構(gòu)和自動表面生成器,面網(wǎng)格增長率1.3,自動修復最小接近值0.01;進水艙室選取切割體網(wǎng)格單元生成器、棱柱層網(wǎng)格生成器,打開表面重構(gòu)和自動表面生成器,面網(wǎng)格增長率1.3,自動修復最小接近值0.01,棱柱層數(shù)1,棱柱層延伸1.5,棱柱層總厚度絕對尺寸0.005 m.考慮運算時間,在重疊域之間建立重疊網(wǎng)格交界面.求解器設置為隱式不定常、六自由度求解器,六自由度運動、負載平衡、分離流、分離VOF波、k-ε湍流、k-ε湍流黏度.
總網(wǎng)格數(shù)量為699 301,檢查網(wǎng)格質(zhì)量,質(zhì)量較好的網(wǎng)格達到98.832%.
機艙網(wǎng)格劃分總網(wǎng)格數(shù)量為707 564,其中質(zhì)量較好的為699 301,占98.832%,如圖14所示;破口處網(wǎng)格加密如圖15所示.
圖14 機艙網(wǎng)格劃分
圖15 機艙破口網(wǎng)格細節(jié)
船體網(wǎng)格劃分如圖16所示,劃分總網(wǎng)格數(shù)量為373 556,其中質(zhì)量較好的為350 977,占93.956%.
圖16 船體網(wǎng)格劃分
變形區(qū)域網(wǎng)格劃分如圖17所示,劃分總數(shù)量為795 358,其中質(zhì)量較好的網(wǎng)格數(shù)量為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,重疊網(wǎng)格區(qū)域基本尺寸為56 m.設定0.05 s為時間步長,對于邊界條件的設置,頂部邊界設置為TOP,底部邊界設置為BOTTOM.
圖17 變形區(qū)域網(wǎng)格加密
DFBI設置六自由度體,啟用模型VOF波,考慮重力因素,兩層全y+壁面處理,精確壁面距離.選擇可實現(xiàn)的k-ε兩層模型,常用k-ε湍流模型,雷諾平均Navier-Stokes模型,三維隱式不定常模型.選取歐拉多相流模型,水設置為恒密度,空氣部分也設置為恒密度,動力黏度和密度均設置為常數(shù).右邊界破口設置為速度入口(INLET A),頂部通風口設置為壓力出口(INLET B).監(jiān)視器監(jiān)視模型六自由度的變化,并生成變化曲線.具體采用STAR-CCM+中的混合壁面函數(shù)功能,指定速度入口的函數(shù)值、壓力出口的函數(shù)值,控制計算中邊界處多相流分布和速度、壓力分布等.
與標準模型實驗結(jié)果進行對比,已經(jīng)證明STAR-CCM+在處理船舶穩(wěn)性研究中可以得到較為可信的結(jié)果.采用實船數(shù)據(jù)利用NAPA進行建模得到實船模型,在破艙進水的過程中,不考慮進艙水的晃蕩影響.進水過程瞬時變化如圖18所示.
模擬過程中可以觀察到在瞬時過程的進水前期,t=8.20 s時,破口處形成水柱噴射,沖擊在艙壁上,機艙內(nèi)產(chǎn)生自由液面,并在船艙底部形成兩處渦流,艙內(nèi)液面逐漸升高;t=80.00 s時,艙內(nèi)液面高度超過破口高度,艙內(nèi)進水不形成水柱噴射,但進艙水仍然有渦流.
監(jiān)測破口處流量隨時間變化,生成破口流量變化曲線,如圖19所示.可以看出在t=0.07 s時,破口流量為655.81 kg/s.在t=0.07~55.00 s時,由于船體橫搖影響,破口流量也產(chǎn)生規(guī)則波動,隨著進艙水的增加,橫搖逐漸趨于穩(wěn)定,破口流量也逐漸趨于平穩(wěn).t=100.00 s時升沉運動持續(xù)增大,破口到水面垂直距離增加,破口流量也隨之增加.對破口流量曲線進行積分,使用科研繪圖工具GraphPad Prism求出曲線下面積,機艙浸滿水時,艙內(nèi)水有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所示,船身向進水一側(cè)移動,橫蕩距離增加,船體形成較大的橫傾角.
第1次橫傾角最大值出現(xiàn)在t=0.66 s,偏移角度為1.29°,是艙室浸滿水時的瞬時橫傾角的16倍,橫搖周期約為5.080 s,比設計計算的橫搖周期5.199 s略?。?/p>
隨著進水時間增加,橫傾角逐漸減小,艙室內(nèi)液面升高,液面趨于平穩(wěn),機艙進水增加,船舶下沉,艏搖加快,機艙進水狀態(tài)趨于平穩(wěn).由于進艙水影響,縱搖逐漸減小,艙室浸滿水約160 s.
(a) 橫搖
(b) 橫蕩
(a) 縱搖
(b) 縱蕩
(a) 艏搖
(b) 升沉
本文利用CFD軟件STAR-CCM+平臺,運用動網(wǎng)格和重疊網(wǎng)格技術對靜水狀態(tài)下的實船模型進行破艙進水模擬.為了保證模擬方法的準確性,先對網(wǎng)格和時間步長收斂性進行分析,然后與標準箱型駁船實驗結(jié)果進行對比,驗證了STAR-CCM+應用于模擬船舶破艙進水的可行性,模擬結(jié)果準確可靠.STAR-CCM+的六自由度求解器、VOF波形模型、可實現(xiàn)的k-ε兩層模型、歐拉多相流模型均適用于船舶實驗的模擬.通過模擬實驗,可以觀察船艙進水的瞬時狀態(tài),通過對監(jiān)視器的六自由度變化曲線進行分析,能夠得到船舶破損時的穩(wěn)態(tài)瞬時數(shù)據(jù),依靠軟件強大的后處理功能,可將模擬結(jié)果生成圖表或視頻,轉(zhuǎn)為可視化數(shù)據(jù).
在確定了一艘漁政船破損位置、尺寸等因素的情況下,對漁政船進行破艙進水時域模擬得到漁政船每一時間步長破損艙室進水位置云圖和破艙進水破口流量及六自由度隨時間變化曲線.剛開始進水時船舶產(chǎn)生較大橫搖角,破口流量規(guī)則波動,隨著進水量的增加船體下沉,橫搖減弱并逐漸穩(wěn)定,破口流量增大,縱搖隨進水量增加逐漸減弱,模擬在艙室浸滿水后停止.
本文實驗結(jié)果可為研究復雜工況如多個艙室破損進水、大風大浪環(huán)境下的船舶破艙進水模擬提供基礎數(shù)據(jù),并為研究船舶破艙進水提供一種新思路和方法.其分析結(jié)果可為船體優(yōu)化、海難船舶救援提供參考.