基于模糊理论的船舶进水速度计算方法

2021-11-05边金宁张海鹏陈淼韩涛

哈尔滨工程大学学报 2021年10期

边金宁，张海鹏，陈淼，韩涛

(哈尔滨工程大学船舶工程学院，黑龙江哈尔滨 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。此值较采取主因素评估算子的值小。从安全性富裕方面考虑，在研究船舶破舱进水方面，选取中部的正方形破口有一定代表性，同时结果有一定的安全富裕。