基于CFD的渔船破舱进水时域模拟及破舱稳性和压强分布研究

于 博 文，张 维 英，陈 静，于 欣，于 洋，周 俊 秋

(大连海洋大学 航海与船舶工程学院，辽宁 大连 116023 )

0 引 言

据中国渔业互保协会统计，1994～2015年的渔船事故中，全损(无法修复)渔船2 996艘，部分(可修复)损失渔船76 454艘．事故原因方面，由人员操作不当造成的碰撞事故占比44.81%，触损事故占比39.12%，触礁事故占比49.23%，搁浅事故占比36.04%[1]．所以探讨渔船在破舱条件下的稳性，尤其是非对称进水过程的六自由度变化，可以判断渔船的分舱合理性和给出破舱应急操作的建议，以达到避免或减少渔船事故损失的目的，对渔船的安全性具有实用价值[2]．

不断发生的海难事件促使人们不断发展船舶破舱稳性的研究，在模型试验之外对于破舱进水的研究方法可以分为两种．准静态法：将整体分成破舱管涌、破损舱室内的水运动以及波浪运动下的船舶运动[3-9]．在这些研究中，基于伯努利方程水动力模型被用来确定通过开口的流速．假定涌入舱内的水是瞬间沉降到一个水位[3，5-7]，并作为块状物在特定路径表面上自由移动[8-9]，或者晃荡并且遵循浅水方程[4]．势流理论法：应用切片理论[3，4，6]或面元法[5，7-9]计算船-波相互作用动力学．计算机硬件的飞速发展，增强了数值模拟用于解决非线性问题的能力．破舱进水过程中水的涌入、冲击、飞溅等都是典型的非线性状况，通过数值模拟方法解决非线性问题的应用日趋成熟，对基于计算流体动力学方法的高保真度模拟的需求也在增加．

刘强等[10]基于Fluent软件用Navier-Stokes(N-S)方程的CFD求解器对进水以及二维舱室耦合运动进行模拟，从舱内空气对水进入舱内的影响分析，提出了减少外域计算量的方法，并探讨了空气压缩性对于进水过程的影响．卢俊尹[11]用模型试验法和一种雷诺平均的N-S(RANS)方程对船舶固定状态下不考虑外流体域的破舱进水问题进行了研究，两者结果吻合较好．Manderbacka等[12]为了增强水流和晃动的耦合效应，对集总质量模型进行了改进，使其能够考虑入流冲量和流量对分隔室内湍流的影响，并且具有计算效率高和实用性强的优点．Sadat-Hosseini等[13]使用CFD中的Ship-Iowa代码(RANS方程求解器)，通过水平设置方案和动态重叠网格技术，模拟了破浪船在横波海况中的淹没过程．Ming等[14]利用光滑粒子流体动力学(SPH)研究相对于入射波的损坏位置对船舱行为的影响．Cao等[15]在2019年建立了一个基于条带理论的振动模型，用于计算不同损伤和荷载条件下受损舱室的运动响应．Gao等[16]在2020年基于RANS方程结合流体体积法(VOF)、滑动界面和动力数值模拟，并采用分层技术处理网格更新，对护卫舰DTMB-5415在静水、规则波中的运动进行了模拟，理论计算和实验结果高度吻合．

CFD逐渐与实验流体力学一起成为船舶稳性研究中的重要手段，其动网格和重叠网格技术在对于非稳态运动物体的运动和强烈的非线性动力学的模拟中，有着很好的应用．

本文对一条养殖看护船DLHT 892 L288先用COMPASS进行稳性计算，再选取一个舱室并用UG NX建模，在考虑内外流体影响下，用建立的STAR-CCM+数值模拟方法计算船体和进水之间横摇的耦合运动，并对船模破舱进水进行时域模拟．通过对破损船和波浪的横向耦合，尽可能还原作业环境下渔船的稳性变化和压力分布，对船舶破舱后扶正及减损提供参考．

1 数学模型

马峥等[17]基于Fluent软件分别用Standardk-ε、RNGk-ε、SSTk-ω、Realizablek-ω和Standardk-ω湍流模型在同一条件下对不同船舶进行模拟并分析．对于方形系数较小、航速较高的船，Standardk-ε和SSTk-ω模型有相对较好的预测精度．SSTk-ω模型对方形系数较大、航速较低的船的总阻力也具有较高的预测精度．而RNGk-ε模型对剩余阻力具有较强的预测能力．综上所述，本文选取Standardk-ε模型．

1.1 湍流模型

1.1.1 不可压缩黏性流体控制方程 本文以三维不可压缩的黏性流体瞬态运动方程为理论基础，流体密度和黏性系数为常数．

质量守恒方程(连续性方程)为

(1)

式中：u、v、w分别为速度矢量v在x、y、z轴上的分量．

动量守恒方程(运动方程)为

(2)

式中：F为质量力，p为压强，μ为流体动力黏度．

有限体积法是在控制体积内对一般形式的控制微分方程的积分，即求解积分形式的守恒方程：

(3)

式中：φ为通用变量，V为控制体积，Γ为广义扩散系数，S为广义源项．

1.1.2 湍流方程 湍流的物理模型基础是三维非定常的，这使得模拟湍流的计算量非常大．前人为了在有限计算资源下模拟湍流，提出了大涡模拟(large eddy simulation，LES)和RANS两种方法．RANS法适用于本文．

连续方程

(4a)

动量方程

式中：u′i表示略去平均符号的雷诺平均速度分量，ρ为密度，u′j为脉动速度，σij为应力张量分量；i=1，2，3．

k-ε模型湍流控制方程如下：

湍流动能方程

Gb-ρε-YM+Sk

(5a)

扩散方程

式中：Gk为速度梯度产生的湍动能项；Gb为浮力产生的湍动能项；i，j=1，2，3分别表示x、y、z方向；YM为脉动扩张项，C1ε、C2ε、C3ε为经验常数；σk、σε分别为与湍动能k和耗散率相对应的Prandtl数；Sk和Sε为用户自定义的源项．

1.2 运动控制方程

船舶在海上运动时，把船体运动看成刚体运动，其平移运动通过高斯-赛德尔(GS)描述的线性方程求解：

(6)

船舶的旋转运动受以下固定于质心的相对坐标系(BS)中的动量方程控制：

(7)

式中：Jc是船舶相对于中心的惯性矩张量，并且相对于BS保持恒定；Ω是船的角速度矢量；M′c是相对于中心作用在船上力矩的合成矢量．

1.3 入射波的产生与消散

通过背景区域的边界条件规定入射波方向，入射波由斯托克斯(Stokes)波理论确定．

水平速度方程

u=Aωcos(K·x-ωt)ekz

(8)

垂直速度方程

w=Aωsin(K·x-ωt)ekz

(9)

表面高度方程

η=Acos(K·x-ωt)

(10)

式中：A为波幅值，ω为波频率，K为波矢量，k为波矢量的幅值，z为与平均水位的垂直距离．

为了减少波在反射边界附近的振荡，可以在选定边界附近添加阻尼衰减．阻尼将垂直阻力引入垂直运动，以避免其在边界反射．Choi等的方法[18]为垂直速度w方程添加了一个阻力项：

(11)

其中

(12)

式中：xsd为波阻尼(在x方向上传播)的起点，xed为波阻尼的终点(边界)，l为阻尼波长度，f1、f2和nd为阻尼模型的参数．

2 模型建立与计算域选取

2.1 实验船和破损舱室建模

在本研究中，养殖看护船DLHT 892 L288的主要参数如表1所示．

表1 养殖看护船DLHT 892 L288的主要参数Tab.1 Main parameters of breeding care boat DLHT 892 L288

依据型值表给出的船型参数，在建模软件UG-NX中，以尾垂线与BL的交点为坐标原点，船头为x轴正方向，左舷为y轴正方向，铅垂向上为z轴正方向建立基准笛卡儿坐标系．依据型线，采用1∶1进行网格曲面的建立，各曲面之间采用G1连续相切，如图1(a)所示；在同一个坐标系中，选取该船的工具舱(舱容为28.00 m3)并对该舱室进行建模，舷侧板厚为0.01 m，如图1(b)所示．

将模型导入STAR-CCM+中，并构建该模型破损舱室顶部的通风孔，其直径为1.50 m，保证船舱进水时，排除内部气压对进水的影响；其破损进水口的直径为0.30 m，如图2所示．

图2 受损船DLHT 892 L288的数值模拟模型Fig.2 Model of damaged boat DLHT 892 L288 in the numerical simulation

(a)船体

船舱的内表面与船外表面均为混合壁面，其切向速度固定，剪应力无滑移，壁面规格平滑．

为较好地捕捉进水流量特征，进水口的最小尺寸为0.005 m，最大尺寸为0.010 m，对舱内区域的空气和水体积采用各向异性加密的方法，其网格尺寸的z方向为0.050 m，x和y方向均为0.200 m．为合理表现出船外表面高曲率部分的曲面特征，表面曲率最大点数为200，其面网格最小尺寸为0.010 m，最大尺寸为0.400 m，划分结果如图3所示．

(a)局部加密

本文选取DLHT 892 L288的装载情况为空载到港，部分装载情况为No.1压载舱100%．No.2压载舱60%，船员100%，淡水10%，燃油10%，食品及备品10%．该装载下排水量为109.4 t，平均型吃水1.575 m．

通过COMPASS计算该船空载返港状况下船舶复原力臂曲线，如图4所示，横坐标为横倾角θ，纵坐标为静稳性力臂Ls．

图4 DLHT 892 L288空载返港状态下的复原力臂曲线Fig.4 DLHT 892 L288 restoring force arm curve under no-load return to port state

2.2 外计算域

2.2.1 外计算域划分 在本文中，监测破损船分别在静水、常规海况和危险海况作业下的六自由度变化．为更好地响应复杂运动，采用嵌套网格技术和VOF法用于破损船六自由度运动的网格更新．

将整个区域分为6个计算域，如图5所示．其中区域1为船舶运动的背景域；区域2为船舶运动的滑道域；区域3为嵌套网格域，其与船舶和舱室区域同步运动，产生相同的平移和旋转；区域4为波浪域；区域5是船舶本体域；区域6为破损舱室域；破损船的外计算域尺寸见表2．

图5 计算域划分Fig.5 Partition of the computational domain

表2 计算域的大小Tab.2 Size of computational domain

2.2.2 外计算域网格 区域1的网格基础尺寸为0.5 m，应用尺寸为2 m，并设置为各向同性．为更好地响应破损船运动，添加嵌套网格为区域3，基础尺寸为0.5 m，其面网格最小尺寸为0.05 m，最大尺寸为2 m，并设置为各向同性．在嵌套网格和背景网格之间添加一个区域2为滑道域，其尺寸为0.25 m，作为区域1和区域3的过渡层，使网格变化递进协调．

在网格控制中，体积增长率设为非常慢，面网格增长率为1.1，表面曲率为72，最小接近修复值为0.01．经以上设置，取该网格设置方案y向截面，如图6所示，各区域网格过渡较好．

经网格检测，在面网格中，高质量面网格数量为8 329 281，占比为100%；体网格中，较好质量体网格数量为18 464，占比为0.222%，高质量体网格数量为8 310 818，占比为99.778%．

2.2.3 波浪区网格 该船现作业区域为渤海范围，结合渤海海域波浪特征[19]，采用以下3种不同的波浪形态：Wave 1(静水)、Wave 2(常规)、Wave 3(恶劣)，分别对破损失速养殖看护船DLHT 892 L288进行稳性分析，其波浪参数见表3．

表3 数值模拟中入射波的参数Tab.3 Parameters of incident waves in the numerical simulation

为保证所加载波浪的准确性，采用二维溃坝方式对Wave 2和Wave 3进行了无风条件下波浪形态检验，网格尺寸的长度分别为对应波长的1/80，高度为对应波高的1/20；时间步Ts为0.01 s(时间步≤P/2.4n，其中P为波周期，n为用于划分波长的段数，为方便与破损船的同步流固耦合运动，取0.01 s)．

所加载的波浪运动的残差图如图7所示，其残差值ξ不大，趋于稳定，有较好的收敛性，该波浪网格方案合理可用．

(a)Wave 2

模拟波所呈现出来的形态分别如图8和9所示，其模拟波浪的形态与阻尼消波情况较好．

图8 阻尼条件下模拟Wave 2的波浪形态Fig.8 Wave shape of simulated Wave 2 under damping condition

图9 阻尼条件下模拟Wave 3的波浪形态Fig.9 Wave shape of simulated Wave 3 under damping condition

3 模拟结果与讨论

3.1 静水状态下破舱失速状态

Wave 1海况下的下模拟网格运动残差分析如图10所示．其残差值较小，趋势具有较好的稳定性，该网格划分具有较好的合理性并满足计算需要．

图10 Wave 1条件下的残差图Fig.10 Residual plot under Wave 1 condition

在Wave 1条件下，外域湍流由破口涌入破损舱室，在考虑进舱水晃荡影响下，其进水过程如图11所示．

(a)t=0.01 s

Wave 1条件下，破口侧的舷侧压力分布如图12所示．

(a)t=0.01 s

Wave 1条件下，相对于初始状态，受损船DLHT 892 L288的六自由度随时间变化如图13所示．

(a)横摇

Wave 1条件下的破口流量如图14所示．

图14 Wave 1条件下破口质量流量随时间变化Fig.14 The mass flow rate changes of the break under Wave 1 condition with time

在Wave 1条件下，在进水阶段的0.56 s时，从破口涌入的水喷射到舱壁，对舱壁产生冲击并开始形成新的自由液面，水流冲击和新的湍流运动使静止破损船产生横摇并在10.50 s达到最大，随后开始衰减．破口进水在0.81 s趋于稳定，舱内进水形成涡流，初期湍流的涌入对船舶的六自由度影响也趋于稳定，随后产生规律稳定的响应．破口舷侧的压力分布稳定．船尾为压强变化最大区域，在图12中，从初始的18 459 Pa在20 s内变化到19 357 Pa，变化幅度很小，对船体影响不大．破损船DLHT 892 L288在静水条件下的六自由度变化很小，横摇响应周期较长，舷侧压力没有较大变化，此状态下的船舶稳性较好．

3.2 Wave 2波浪状态下的破舱失速状态

Wave 2海况下的下模拟网格运动残差分析如图15所示．其残差值较小，趋势具有较好的稳定性，该网格划分具有较好的合理性并满足计算需要．

图15 Wave 2条件下的残差图Fig.15 Residual plot under Wave 2 condition

假定船初始位于Wave 2的波峰处，此时横波在船处浪位最高，其船随波浪产生的刚体运动变化如图16所示．Wave 2波浪载况下破口侧舷侧压力分布如图17所示．Wave 2工况下，破口涌入的海水响应比Wave 1条件下要强很多，故在此捕捉进水的瞬时状态，在考虑进舱水的晃荡影响下，进水过程的瞬时变化如图18所示．

(a)t=0.50 s

(a)t=0.50 s

(a)t=0.01 s

Wave 2条件下，相对于初始状态，受损船DLHT 892 L288的六自由度随时间变化如图19所示．

Wave 2条件下的破口质量流量如图20所示．

图20 Wave 2条件下破口质量流量随时间变化Fig.20 The mass flow rate changes of the break under Wave 2 condition with time

在Wave 2条件下，破舱管涌形成蘑菇状射流冲击舱壁，结合波浪作用，其艏摇在1.30 s达到6.8°，横摇在第一周期达到7.0°，并在持续冲击状态中于8.00 s前后达到15.0°的倾斜位置，随后逐渐收敛．纵荡以及横荡在波浪与破损船的耦合作用下产生了较大的变化．其压力最大处位于龙骨，随纵摇前后移动，最大压强变化区间为22 158～15 979 Pa，如图17所示，变化幅度相对平缓，无中拱中垂现象．此条件下破舱后的六自由度变化虽然较大，但船舶无倾覆趋势，最大横摇为最大横倾角一半，稳性较好，比较安全．

(a)横摇

3.3 Wave 3波浪状态下的破舱失速状态

Wave 3海况下的下模拟网格运动残差分析如图21所示．其残差值较小，趋势具有较好的稳定性，该网格划分具有较好的合理性并满足计算需要．

图21 Wave 3条件下的残差图Fig.21 Residual plot under Wave 3 condition

假定船初始位于Wave 3的波峰处，此时横波在船处浪位最高，其船随波浪刚体运动变化如图22所示．

Wave 3波浪载况下破口侧舷侧压力分布如图23所示．

在Wave 3海况下，破口涌入的海水响应比Wave 2海况下的要剧烈很多，故捕捉其短时间内进水的瞬时状况，在考虑进舱水的晃荡影响下，进水过程的瞬时变化如图24所示．

(a)t=0.20 s

(a)t=0.20 s

(a)t=0.01 s

Wave 3条件下，相对于初始状态，受损船DLHT 892 L288的六自由度随时间变化如图25所示．

Wave 3条件下的破口质量流量如图26所示．

(a)横摇

图26 Wave 3条件下破口质量流量随时间变化Fig.26 The mass flow rate changes of the break under Wave 3 condition with time

在Wave 3条件下，破口处初始质量流量高达42 000 kg/s，破舱管涌对舱壁产生剧烈冲击，新自由液面的晃荡再结合外域波浪的推动，使其横摇在1.10 s达到75°，纵摇在3.10 s达到12°，并且产生了较为严重的甲板上浪情况．虽然产生了巨大的横倾后趋于收敛，这得益于艏摇的快速转变，使得破损船舷侧受波浪冲击的相对面积大幅度减小，但如此夸张的横倾角已经远超该船的最大横倾角，依赖波浪作用才不至倾覆．由于出现的艏摇运动使该船从遭受横浪影响转为遭受纵浪影响，且在图23(b)中，龙骨中段最大压强为40 180.0 Pa，两端最小处为8 042.7 Pa，容易造成中拱现象．运动过程中，龙骨处最大压强变化区间为40 180.0～20 810.0 Pa，该变化剧烈且强度大，容易给船体结构带来一定损伤．此条件下破舱后的六自由度变化剧烈，并且伴随较大的甲板上浪和中拱现象，虽然不至于倾覆，但是极度危险．

4 结 语

沿海渔船作为高发事故群体，其破舱后失稳容易造成船体以外的各项损失．通过对比该船在Wave 1、Wave 2、Wave 3条件下的六自由度变化，可以发现在进水阶段的过渡相中，破损船会向破口侧倾斜，破舱管涌冲击舱壁使船反向侧倾，凭借复原力矩产生横摇运动，其运动剧烈程度与破舱管涌的强弱相关．在初始进水流量中Wave 1、Wave 2、Wave 3流量比约为1∶20∶40，结合风浪荷载产生的最大横摇比约为1∶10∶30，后者甚至达到了夸张的75°，按国际拖船大会的建议[20]，该船已经属于倾覆状态．幸运的是该破损舱室位于船头，通过观察舷侧压强分布变化，发现产生的艏摇运动在波浪的推动下使船头朝向改变为风浪前进方向，进而较大地减少了横摇运动的幅度．

综上所述，沿海渔船的静稳性对于作业时破舱进水后物理行为的参考意义不大．考虑船上人员对于进水的第一阶段难以及时反应，发现船舱破损时，进水阶段往往处于渐进相，此时可通过调转船头方向，使船尽可能处于纵波状态，以减小风浪和破舱管涌带来的影响．

对于渔船舱室，尽量设置在渔船中横剖线两侧，当发生破舱时，借助破舱进水的湍流冲击和波浪耦合作用，推动破损船尽快从横波状态变为纵波状态，以减小横摇幅度，减少损失．