基于模糊理論的船舶進水速度計算方法

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

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

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

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

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

1 進水速度評估中的模糊模型

1.1 模糊理論原理

Zadeh教授提出了隸屬函數概念，奠定了模糊理論的基礎。基于模糊理論發展起來的模糊綜合評價法是由汪莊培提出來的。該方法可基于模糊數學對受到多因素影響的事物做出一個總體評價。主要包括以下步驟[12]：

1)確定模糊評價等級。

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

2)權重確定。

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

3)建立隸屬度向量。

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

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

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

4)模糊算子確定。

得到隸屬度和權重后，需根據評價指標的特點選擇模糊算子。常見模糊算子如下：

主因素決定型：

主因素突出型：

不均衡平均型：

加權平均形：

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

5)綜合評估。

1.2 基于模糊理論的進水速度評估

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

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

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

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

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

2 數值仿真模擬

2.1 可行性分析

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

模型建立過程中最重要、最難的部分就是破口面的處理。研究CFD的各種邊界條件發現一種特殊邊界：interior邊界，此邊界允許流體通過，可很好體現破口這一情況。沿破口做一截面進行觀察。可看出破口處不存在任何邊界限制如圖2。

圖2 截面邊界Fig.2 Section boundary

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

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

設置艙室中心位置為水位監測點監測水位。分析不同時間的進水相位圖如圖4。

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

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

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

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

圖6 監測點水位Fig.6 Water level map of monitoring points

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

2.2 實驗方案

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

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

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

2.2.1 收斂性分析

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

表1 網格信息表Table 1 Grid information table

船艙吃水為4.5 m，設置破口處截面進行觀察。1、5、10 s不同網格尺寸的相位分布如圖7。

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

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

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

2.2.2 實驗方案

1)破口面積

為研究破口面積對進水速度的影響，選破口形狀均為正方形，中心位置均為(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)破口形狀

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

表2 不同形狀破口信息表

4)破口位置

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

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

3 單因素敏感性分析

逐一計算上述破損案例，監測進水艙內的質量流量，分析各個案例的進水速度。

3.1 破口面積

獲取進水速度，計算平均進水速度。v=m/tA。其中，m為質量流量，t為對應的時間，A為破口面積。

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

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

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

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

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

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

3.2 破口高度

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

計算各個時刻的平均進水速度如圖13。

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

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

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

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

3.3 破口形狀

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

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

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

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

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

3.4 破口位置

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

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

統計不同模型前3、5、7、9 s內的平均速度值，以模型序號為橫坐標、速度值為縱坐標，定性分析不同時間節點處進水速度隨破口位置的變化情況如圖18。

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

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

4 模糊理論評估進水速度算例

4.1 評價體系

本次評估模型如表4。

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

4.2 評語集及權重確定

1)第1層權重確定：

2)第2層破口形狀權重W1:

將統計的33個破口形狀進行分類，如表5。

表5 破口形狀統計表Table 5 Crevasse shape statistic table

以不同形狀破口在統計案例中出現的頻率來確定破口形狀的權重，則W1=[0.151 5 0.515 2 0.060 6 0.212 1 0.060 6]

3)第2層破口位置權重W2：

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

表6 評語集對應表Table 6 Corresponding table of comment set m/s

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

4.3 隸屬度及模糊算子確定

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

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

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

4.4 進水速度綜評估

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

故：

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

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

5 結論

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

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

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