邊金寧， 張海鵬， 陳淼， 韓濤

(哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001)

艦船由于破損進(jìn)水造成的安全性問題一直是行業(yè)屆關(guān)注的焦點[1]。研究進(jìn)水速度的影響因素對深入研究進(jìn)水機(jī)理、估算沉沒時間、評估人員疏散可用時間具有重要意義。船舶進(jìn)水速度的研究極其復(fù)雜，研究方法從理論推導(dǎo)轉(zhuǎn)向數(shù)值模擬。胡麗芬等[2]基于伯努利方程對進(jìn)水過程進(jìn)行建模模擬，得到速度緩慢增大再減小的變化趨勢。鄢凱等[3]針對多艙室建立了快速求解的解析算法。但理論推導(dǎo)受到多方面限制：伯努利方程只適用于小開口，真實海難中船舶開口基本均為大開口、該方法無法考慮破口信息的復(fù)雜性、水的隨機(jī)性蔓延等。隨著計算機(jī)技術(shù)的發(fā)展，科研人員開始利用CFD去解決流體問題[4]。Gao等[5]對橫浪中的船舶進(jìn)水進(jìn)行了數(shù)值模擬。鄭宇[6]、李月萌等[7]利用CFD技術(shù)模擬了客船、箱型多艙室的進(jìn)水過程。但上述學(xué)者僅驗證了CFD技術(shù)對船舶進(jìn)水模擬的可行性，沒有細(xì)致探討影響進(jìn)水速度的因素。船舶進(jìn)水是一個多因素的復(fù)雜問題，破口的位置、面積、形狀等都會影響進(jìn)水過程，不同因素之間還存在互相影響。

模糊理論處理模糊信息能力強(qiáng)、定性信息定量化研究方便[8]，容錯性高可將定性問題根據(jù)實際情況進(jìn)行綜合運(yùn)算，譚志榮等[9]利用模糊評價法對客船的安全隱患進(jìn)行評估，焦甲龍等[10]應(yīng)用多級模糊綜合評價模型對風(fēng)浪環(huán)境下的水面艦艇適應(yīng)性進(jìn)行評估。孫麗萍等[11]用模糊語言對原油艙的風(fēng)險進(jìn)行研究。模糊理論已經(jīng)被應(yīng)用在船舶風(fēng)險評估、安全評價等方面。其處理不確定性因素的能力已被學(xué)者們所證實。

本文以箱型住艙為研究對象，采用流體體積函數(shù)技術(shù)、明渠流動技術(shù)實現(xiàn)進(jìn)水過程的可視化，獲取進(jìn)水速度。首先對本次模擬所需的技術(shù)進(jìn)行可行性分析，從破口面積、高度、形狀、縱向位置4個方面對不同的模型進(jìn)行計算。定量分析破口面積、高度對進(jìn)水速度的影響。構(gòu)造基于模糊理論的評估模型，采用層次分析(analytic hierarchy process,AHP)方法確定權(quán)重，分析破口形狀及縱向位置對進(jìn)水速度的影響，最后采用不同模糊算子評估進(jìn)水速度。

1 進(jìn)水速度評估中的模糊模型

1.1 模糊理論原理

Zadeh教授提出了隸屬函數(shù)概念，奠定了模糊理論的基礎(chǔ)?；谀：碚摪l(fā)展起來的模糊綜合評價法是由汪莊培提出來的。該方法可基于模糊數(shù)學(xué)對受到多因素影響的事物做出一個總體評價。主要包括以下步驟[12]：

1)確定模糊評價等級。

通常模糊評語集設(shè)定為h={h1,h2,…，hm}。其中hi(i=1,2,…，m)代表不同的級別。

2)權(quán)重確定。

確定各個因素對上層因素的影響程度，權(quán)重的確定方法直接影響評價結(jié)果，目前比較流行的有：專家評議法、專家調(diào)查法、層次分析法等。

3)建立隸屬度向量。

模糊理論中的元素對于模糊集來說不僅僅是屬于和不屬于2種，隸屬度值只能在0～1取值。

μA:U→[0,1]u|→μA(u)

式中：μA(u)是集合A的隸屬度函數(shù)，μA(ui)中ui為i的的隸屬度，將隸屬度進(jìn)行歸一化處理即得到單因素的隸屬度向量ri。將所有指標(biāo)的隸屬度向量進(jìn)行組合可得到隸屬度矩陣R。

4)模糊算子確定。

得到隸屬度和權(quán)重后，需根據(jù)評價指標(biāo)的特點選擇模糊算子。常見模糊算子如下：

主因素決定型：

主因素突出型：

不均衡平均型：

加權(quán)平均形：

第1種算子只考慮主要因素的影響，忽視其余因素，會放大主要影響因素的影響效果。第2、3種算法對第1種算法進(jìn)行細(xì)分，考慮的因素更全面。第4種算子考慮了不同種因素的影響程度，對影響程度進(jìn)行了加權(quán)計算，評價結(jié)果更加的客觀合理。

5)綜合評估。

1.2 基于模糊理論的進(jìn)水速度評估

在實際海難發(fā)生前，破口的形狀、位置是不確定的。相同面積、高度時，不同形狀、位置的破口進(jìn)水速度有較大差距。破口位置及形狀2個因素很難進(jìn)行量化處理，本文將采用模糊綜合評價法，以進(jìn)水速度為研究目標(biāo)，將破口形狀和位置進(jìn)行量化處理，綜合各種不確定因素得到艙室的進(jìn)水速度。模糊綜合評價的邏輯如圖1。

1)AHP可用于船舶的多因素決策問題，同時可得到各因素對目標(biāo)層的影響權(quán)重[13]。參考AHP的層次結(jié)構(gòu)模型，確定第1層為綜合進(jìn)水速度，第2層為破口形狀、破口位置；第3層為不同的破口案例，為量化不同因素對進(jìn)水速度的影響程度，用速度值來確定評語集;

2) 利用AHP法的重要性原則確定第2層的權(quán)重，重要性原則將各因素的重要性進(jìn)行數(shù)值標(biāo)度(1～9)[14]；1代表二者對于研究目標(biāo)來說同等重要，數(shù)值越大，前者相比后者來說重要性越強(qiáng)。針對第3層權(quán)重：采用統(tǒng)計學(xué)理念，利用不同形式破口在真實海難上出現(xiàn)的頻率來近似代替概率;

3)利用CFD技術(shù)提供各種案例的進(jìn)水速度，參照評語集將速度值換算成隸屬度矩陣；

4)選取比較有代表性的主因素決定性算子和加權(quán)平均形算子進(jìn)行模糊計算。

2 數(shù)值仿真模擬

2.1 可行性分析

利用一簡單箱型艙室進(jìn)行可行性驗證。艙室尺寸為1 m×1 m×1 m，破口為正方形，尺寸為0.2 m×0.2 m，位于右舷側(cè)。模擬過程中必須考慮外界流域的作用。取計算外域為4 m×3 m×4 m。

模型建立過程中最重要、最難的部分就是破口面的處理。研究CFD的各種邊界條件發(fā)現(xiàn)一種特殊邊界：interior邊界，此邊界允許流體通過，可很好體現(xiàn)破口這一情況。沿破口做一截面進(jìn)行觀察?？煽闯銎瓶谔幉淮嬖谌魏芜吔缦拗迫鐖D2。

圖2 截面邊界Fig.2 Section boundary

對模型進(jìn)行網(wǎng)格劃分，導(dǎo)入fluent軟件：設(shè)置瞬態(tài)計算，重力方向為Y軸，利用VOF技術(shù)捕捉水氣交界面[15]，明渠流動技術(shù)保證計算過程中外界流場的穩(wěn)定性。利用Patch技術(shù)完成初始流場的填充，初始流場分布如圖3。

圖3 初始流場分布Fig.3 Initial flow distribution

設(shè)置艙室中心位置為水位監(jiān)測點監(jiān)測水位。分析不同時間的進(jìn)水相位圖如圖4。

圖4 截面相位分布Fig.4 The phase distribution of section

1 s時海水在內(nèi)外部壓力和自身重力作用下射向艙內(nèi)。破口上部產(chǎn)生一部分空白區(qū)，為研究此空白區(qū)，顯示1 s時破口處的空氣、水速度分布如圖5，圖中箭頭為速度矢量。將破口位置放大如圖5。此圖可看到破口上部分有艙內(nèi)空氣向外部逃逸。艙內(nèi)空氣被壓縮，被壓縮的空氣從破口處逃逸，將破口處的外界海水排開，在破口上端產(chǎn)生部分空白區(qū)域。此現(xiàn)象驗證了本模型可實現(xiàn)艙內(nèi)空氣逃逸的物理現(xiàn)象。

圖5 破口處速度分布Fig.5 Velocity distribution of crevasse

相位圖分析：1 s海水抵達(dá)艙內(nèi)底部，底部水開始向艙內(nèi)其他空間流動，當(dāng)海水流動至四周艙壁后艙壁將阻擋海水的繼續(xù)蔓延，艙壁將具有一定流度的海水阻擋回艙內(nèi)。相位圖可直觀觀察海水的分布，將海水隨機(jī)性蔓延這一過程進(jìn)行可視化。水位如圖6所示。

圖6 監(jiān)測點水位Fig.6 Water level map of monitoring points

由圖6可知，前40 s水位上升較快，隨后水位上升緩慢最后趨于穩(wěn)定。20 s之前水位有小幅波動。20 s之前進(jìn)水量較少，未形成相對平穩(wěn)的水面，艙內(nèi)水面波動較劇烈；隨著進(jìn)水的持續(xù)，艙內(nèi)形成較平緩水面，進(jìn)水對艙內(nèi)水面的波動較小，水位平緩上升。本模型進(jìn)水量為600 kg，進(jìn)水高度為0.6 m。而外界水位高度為0.8 m。進(jìn)水高度未達(dá)到外流域水位，艙內(nèi)存在空氣墊阻止了進(jìn)水過程。

2.2 實驗方案

進(jìn)水的影響因素主要分2方面；外部因素：海洋環(huán)境、海水流速；內(nèi)部因素：破口尺寸、位置、形狀、深度。不同的航道、時間，海洋環(huán)境都不同，本次主要研究內(nèi)部因素的影響。

進(jìn)水艙以外的船體對于本實驗的影響體現(xiàn)在吃水上。由于本次研究需大量實驗案例作為分析基礎(chǔ)，為減小計算量本次針對進(jìn)水艙進(jìn)行建模，通過設(shè)置吃水來體現(xiàn)船體對于進(jìn)水模型的影響。

以典型箱型住艙為例，艙室尺寸長×寬×高分別為9 m×6 m×5 m。破口為正方形，尺寸為1 m×1 m，在考慮計算量和計算精度的情況下確定計算外域為20 m×10 m×10 m。

2.2.1 收斂性分析

網(wǎng)格的尺寸決定計算精度，理論上網(wǎng)格尺寸越小，計算精度越好，但過多網(wǎng)格將增加計算時間。在計算前需要確定最佳網(wǎng)格尺寸。針對不同區(qū)域采取不同尺寸網(wǎng)格：進(jìn)水艙采用加密網(wǎng)格，設(shè)置3種網(wǎng)格進(jìn)行收斂性分析。網(wǎng)格數(shù)據(jù)如表1。

表1 網(wǎng)格信息表Table 1 Grid information table

船艙吃水為4.5 m，設(shè)置破口處截面進(jìn)行觀察。1、5、10 s不同網(wǎng)格尺寸的相位分布如圖7。

圖7 不同尺寸相位Fig.7 The phase of different size

將3種網(wǎng)格的垂向受力數(shù)據(jù)作圖如圖8。3種尺寸的網(wǎng)格均可實現(xiàn)進(jìn)水的模擬，但模擬效果有一定差距，0.1網(wǎng)格的相位分布圖最清晰。0.12網(wǎng)格的相位圖與0.1很接近。0.2網(wǎng)格的相位分布最模糊，垂向受力與另外2個案例差距較大。故中等網(wǎng)格模擬結(jié)果較好，計算精度、時間可接受。本次模擬取進(jìn)水艙內(nèi)網(wǎng)格尺寸為0.12。

圖8 垂向受力對比Fig.8 Vertical force contrast diagram

2.2.2 實驗方案

1)破口面積

為研究破口面積對進(jìn)水速度的影響，選破口形狀均為正方形，中心位置均為(0,0,3 m)，不同面積破口處的截面圖如圖9。

圖9 不同面積破口截面Fig.9 The crevasse section of different areas

2)破口中心位置

不同中心位置破口截面如圖10。

圖10 不同中心位置破口截面Fig.10 The crevasse section of different positions

3)破口形狀

為確定真實海難中的破損形狀，統(tǒng)計來自國家海事局上的碰撞案例(最新的60個案例)，發(fā)現(xiàn)碰撞后船舶會發(fā)生擦傷、破損、傾覆、斷裂，本文主要研究船體破損產(chǎn)生的破口形狀。統(tǒng)計產(chǎn)生破口的案例，統(tǒng)計出27例有效案例，這些案例中有的2艘船舶均產(chǎn)生破口、有的一艘船舶產(chǎn)生破口、有的一艘船舶產(chǎn)生多處破口，共計33個破口，破口可近似為以下形狀：圓形、縱向長條、垂向長條、斜正方形、正方形。故實驗的破口信息如表2。

表2 不同形狀破口信息表

4)破口位置

研究破口縱向位置對進(jìn)水速度的影響，破口面積為1，形狀為正方形，數(shù)據(jù)如表3。

表3 不同縱向位置破口信息表Table 3 Crevasse information table at different positions

3 單因素敏感性分析

逐一計算上述破損案例，監(jiān)測進(jìn)水艙內(nèi)的質(zhì)量流量，分析各個案例的進(jìn)水速度。

3.1 破口面積

獲取進(jìn)水速度，計算平均進(jìn)水速度。v=m/tA。其中，m為質(zhì)量流量，t為對應(yīng)的時間，A為破口面積。

做平均速度隨時間的變化曲線，如圖11。

圖11 不同面積條件下破口進(jìn)水速度變化Fig.11 The change of inlet velocity under different area conditions

由圖11可知模型一(0.5 m×0.5 m)進(jìn)水速度隨時間增大，1.8 s達(dá)到最大值，然后急劇下降，2.5 s開始繼續(xù)上升，上升至3.5 s后開始平穩(wěn)下降。模型2(1 m×1 m)的速度變化趨勢大致相同。模型3(1.5 m×1.5 m)前期變化趨勢與模型一相同，但后期進(jìn)水速度緩慢上升。2 m×2 m破口的速度變化有小幅波動，在11 s有較大波動，后期速度緩慢降低。由于此破口較大，進(jìn)水量較大，而網(wǎng)格尺寸仍為0.12，導(dǎo)致單位時間內(nèi)流過一個網(wǎng)格的數(shù)據(jù)量較大，造成計算的準(zhǔn)確性下降，產(chǎn)生較大誤差。

為定量得到破口面積與進(jìn)水速度的關(guān)系，統(tǒng)計前5、7、9、11、13、15 s內(nèi)的速度。以破口面積為橫坐標(biāo)、進(jìn)水速度為縱坐標(biāo)，做出不同時間的速度隨破口面積的變化曲線(2 m×2 m模型計算精度較差，不做統(tǒng)計)。如圖12。

圖12 不同面積條件下破口平均進(jìn)水速度變化Fig.12 The change of average inlet velocity under different area conditions

通過此圖可知：在0.25≤S≤2.25時，進(jìn)水速度隨著破口面積的增大呈線性的增長趨勢。

3.2 破口高度

為探討破口中心高度對進(jìn)水速度的影響，本文破口高度選取范圍從吃水位置至艙室底部選取5個位置，保持各個破口之間距離均為0.5 m，保證破口之間距離一致，減少破口距離變化不一致額外增加新因素對進(jìn)水速度造成影響。同時這樣分布破口可保證破口范圍較大，涵蓋水下的所有位置。

計算各個時刻的平均進(jìn)水速度如圖13。

圖13 不同高度條件下破口進(jìn)水速度變化Fig.13 The change of inlet velocity at different heights

前3個模型(0.5 m、0 m，-0.5 m)進(jìn)水速度趨勢及數(shù)值很接近。-1 m高的破口模型最大速度值發(fā)生增加，達(dá)到2.7 m/s。整體趨勢與前3個模型相似，初次上升時間提前到0.5 s，最終速度為2.5 m/s。-1.5 m案例趨勢與-1 m接近。最大值為4.8 m/s。最終數(shù)值與-1 m高案例結(jié)果一致。統(tǒng)計進(jìn)水速度值，并作圖14。

圖14 不同高度條件下破口平均進(jìn)水速度變化Fig.14 The change of average inlet velocity at different heights

前3 s的速度隨破口高度的降低而增大，當(dāng)-1≤H≤0.5時：5 s后速度值不隨高度而變化，-1.5 m高案例的進(jìn)水速度比其他高度的速度大。

3.3 破口形狀

計算各個時刻的平均進(jìn)水速度作圖如圖15。所有模型的進(jìn)水速度變化趨勢相同，速度值不同，進(jìn)水速度由高至低分別為：垂向長條破口、旋轉(zhuǎn)正方形、正方形、圓形破口、最小的為橫向長條破口。最終平穩(wěn)速度分別為3.35、2.23、2.06、1.94、1.31。

圖15 不同形狀破口進(jìn)水速度變化Fig.15 The change of inlet velocity in different breach shapes

圓形和正方形速度值很接近。相同面積時，破口垂向所占比例越大、進(jìn)水速度越大，統(tǒng)計不同模型前3、5、7、9的平均速度。以模型序號為橫坐標(biāo)、速度值為縱坐標(biāo)，分析不同時間節(jié)點處進(jìn)水速度隨破口形狀的變化曲線如圖16。

圖16 不同形狀破口平均進(jìn)水速度變化Fig.16 The change of average inlet velocity in different breach shapes

由圖可知，各曲線變化趨勢相同，各個時間點下模型3的速度最大，模型2速度值最低。

3.4 破口位置

計算各個時刻的平均進(jìn)水速度并作圖17。3種模型進(jìn)水速度變化趨勢一致，數(shù)值上存在微小差異，中部的整體速度值最高，其次是前部、后部，前部和后部的速度值在7 s后基本相同。

圖17 不同位置條件下破口進(jìn)水速度變化Fig.17 The change of inlet velocity at different positions

統(tǒng)計不同模型前3、5、7、9 s內(nèi)的平均速度值，以模型序號為橫坐標(biāo)、速度值為縱坐標(biāo)，定性分析不同時間節(jié)點處進(jìn)水速度隨破口位置的變化情況如圖18。

圖18 不同位置條件下破口平均進(jìn)水速度變化Fig.18 The change of average inlet velocity at different positions

5 s后速度值趨于穩(wěn)定，舯部速度最大，艏部和艉部模型速度值相近，呈對稱形狀分布。

4 模糊理論評估進(jìn)水速度算例

4.1 評價體系

本次評估模型如表4。

表4 評估因素等級表Table 4 Rating factor scale table

4.2 評語集及權(quán)重確定

1)第1層權(quán)重確定：

2)第2層破口形狀權(quán)重W1:

將統(tǒng)計的33個破口形狀進(jìn)行分類，如表5。

表5 破口形狀統(tǒng)計表Table 5 Crevasse shape statistic table

以不同形狀破口在統(tǒng)計案例中出現(xiàn)的頻率來確定破口形狀的權(quán)重，則W1=[0.151 5 0.515 2 0.060 6 0.212 1 0.060 6]

3)第2層破口位置權(quán)重W2：

SOLAS公約是一種基于統(tǒng)計學(xué)得到的公約，對于確定不同位置在真實海難中的破損概率有指導(dǎo)意義[17]。將住艙沿船長方向等分劃分為艉部、舯部、艏部3個區(qū)域。依據(jù)SOLAS公約計算破損概率的方法。得到艉部破損概率0.299 663 ，舯部為 0.265 993 ，艏部為0.299 663 。將破損概率進(jìn)行歸一化處理得到W2,W2=(0.346 5,0.307,0.346 5)評語集的間隔越小，評估的準(zhǔn)確性越高。由數(shù)值模擬可知速度值分布在C=[1.25 1.75 2.25 2.75 3.25]T，參考李克特量表將速度區(qū)間等分為5個等級，如表6所示。

表6 評語集對應(yīng)表Table 6 Corresponding table of comment set m/s

為將最終結(jié)果量化，取評語集分值的中間值C=[1.25 1.75 2.25 2.75 3.25]T代替本等級。

4.3 隸屬度及模糊算子確定

在破口面積為1 m2，高度為0 m情況下，評估出該住艙5 s后穩(wěn)定的進(jìn)水速度。

參考評語集將不同案例的速度值進(jìn)行換算，得到該因素的隸屬度。例如P11因素的隸屬度求解。V=2.02 m/s。對應(yīng)的隸屬度R11=[0 0.46 0.54 0 0]。將所有速度值轉(zhuǎn)換為對應(yīng)的隸屬度，得到破口形狀的隸屬度矩陣R1。同理得到R2。

模糊算子：本文選用主因素決定性算子M(∧,∨)和加權(quán)平均算子M(·,⊕)來評估綜合進(jìn)水速度。

4.4 進(jìn)水速度綜評估

1)考慮主要因素的影響作用，忽略小因素的影響，取模糊算子為主因素決定型算子M(∧,∨)：

故：

2)綜合考慮所有因素的影響作用，采用加權(quán)平均算子M(·,⊕)得到各層相對于評語集的評價向量：

在只考慮主要因素的情況下，進(jìn)水速度評估值為2.088 m/s、與位于中部的正方形破口速度2.07 m/s相近，在充分考慮所有因素后進(jìn)水速度評估值為1.88 m/s。此值較采取主因素評估算子的值小。從安全性富裕方面考慮，在研究船舶破艙進(jìn)水方面，選取中部的正方形破口有一定代表性，同時結(jié)果有一定的安全富裕。

5 結(jié)論

1)進(jìn)水速度隨破口面積線性增大；破口高度(-1.5～0.5 m)對進(jìn)水速度基本無影響；中部區(qū)域進(jìn)水速度稍大，首尾區(qū)域速度相同；不同破口形狀進(jìn)水速度差距較大，垂向占比大的破口進(jìn)水速度大。橫向長條破口發(fā)生的頻率最高。

2)基于模糊評價理論可以評估出艙室的進(jìn)水速度，選取主要因素算子，進(jìn)水速度與中部正方形破口速度值相近。選取加權(quán)平均算子后進(jìn)水速度為1.88 m/s。在船舶進(jìn)水機(jī)理的研究中，可選取中部的正方形破口進(jìn)行進(jìn)水的數(shù)值模擬。

3)本文對研究船舶破艙進(jìn)水過程中破口的選擇、影響船舶進(jìn)水速度因素的探討、基于多因素評估進(jìn)水速度有重大參考意義。