基于STAR-CCM+的船舶破舱进水模拟

王 涛， 张 维 英*， 胡 丽 芬， 于 博 文， 王 树 祥， 鲁 纪 元

( 1.大连海洋大学 航海与船舶工程学院， 辽宁 大连 116023；2.鲁东大学 交通学院， 山东 烟台 264025 )

0 引 言

海难事故多种多样，例如着火或爆炸、进水、碰撞、搁浅、倾斜或倾覆、沉没、失控及漂浮等[1]．海难事故导致的破舱进水会造成严重的经济损失，众多学者对船舶的破舱进水进行了研究．Spouge对一起海难事故进行了调查研究，考察了两艘船在碰撞前后的运动、碰撞本身的机制以及随后破损船舶快速倾斜和下沉的原因，观察到“瞬时不对称进水”现象，并解释了破损船舶快速横倾的原因，讨论了这一现象的后果[2]．

近年来，计算流体力学(CFD)的兴起不仅促进了实验研究和理论分析的发展，也为流动模型的简化提供了更多依据[3]．在船舶流体力学中，利用CFD进行船舶稳性研究已经形成趋势，相比静态稳性，动态稳性的研究更加复杂．Gao等观察船舶破损过程，发现海水涌入受损舱室时，进舱水和受损船舶相互运动．为了解决这种相互运动状态的数值模拟计算，Gao等采用基于自由表面捕获技术的Navier-Stokes(N-S)求解模型开发了一种求解器，并验证了求解器的准确性．最后，对滚装渡船破损舱体进水现象进行了数值模拟，计算结果与实验数据吻合较好[4]．卢俊尹建立了破舱进水的数理模型，通过数理模型计算验证了数值方法的正确性，并分析了影响破舱进水的要素，结果表明其数理模型方法可以模拟实际舱室进水[5]．在研究过程中经常要考虑水、油和空气等自由液面的变化，尤其是船舶在航行过程中处于水、气两层流体中，CFD软件能够利用多种湍流模型和有限体积法模拟这种问题．刘强等对破舱进水时域模拟进行了探讨，分析了空气流对进水过程的影响和破口水流速度变化，提出了一种减少计算域的方法[6]．

吴文锋利用ANSYS软件模拟了船舶碰撞过程，分析了撞击参数对碰撞性能的影响，然后进行了物理模型实验，利用FLUENT对双壳油船在静态水域碰撞后发生泄漏的情况进行三相流数值模拟，将模拟结果与物理实验结果进行对比分析，验证了FLUENT软件对液货泄漏模拟所得结果是可行的[7]．在对动态稳性的研究结果分析前，船舶在静水中的稳性研究相比于复杂海况中的研究更有必要．静水中的船舶稳性研究能够给复杂工况下船舶稳性研究提供基础．李月萌等利用FLUENT对Ruponen的模型进行破舱进水时域模拟，验证了数值模拟方法结果的准确性，并能观察到破舱进水每个时间步长的进水过程变化，为破舱模拟分析提供了一种新的研究方法[8]．

在船舶的破损位置、尺寸、进水量等因素确定的情况下，研究受损船舶的流体运动特性和自身动态行为可为有效减损提供参考，还可对船舶破舱稳性进行计算与评估，为船舶设计提供重要依据[9]．本文利用CFD软件STAR-CCM+平台，建立对一艘箱型驳船模型进行破舱进水两相流的数值模拟模型，进行箱型驳船在静水情况下的破舱进水时域模拟．将模拟结果与模型实验结果进行比较，验证STAR-CCM+对船舶破舱进水研究的可行性．然后根据实船数据利用NAPA建模，选取一个舱室，主要考虑船体和进舱水之间的横摇耦合运动，利用建立的STAR-CCM+数值模拟方法，对实船模型进行破舱进水时域模拟，监控静水状态下船体和水在六自由度下的耦合运动，并对模拟结果进行分析讨论．

1 数学模型

2009年马峥等对船舶数值模拟中的湍流模型进行研究，利用FLUENT对三大主力船型即散货船、油船和集装箱船进行数值模拟，选取5种常用两方程湍流模型，网格数量和划分方法基本一致，实验结果得到：在数值模拟选择湍流模型时，对于集装箱船一类，标准k-ε模型和SSTk-ω模型相对较好．对于油船、散货船等，SSTk-ω模型具有较高的预测精度，而RNGk-ε模型对剩余阻力具有较强的预测能力[10]．

1.1 控制方程

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

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

(1)

式中：u、v、w为速度矢量v沿着x、y、z轴3个方向的速度分量．

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

(2)

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

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

(3)

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

1.2 数值算法

动量守恒方程(N-S方程)为

(4)

式中：p为压强，单位Pa；τxy、τxx、τxz等是黏性应力τ的分量，单位Pa；fx、fy、fz为x、y、z方向上的单位质量力，单位m/s2．

k-ε湍流模型方程为

Gb-ρε-Ym+Sk

(5)

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

采用基于有限元分析的STAR-CCM+软件求解模型内流场．不可压缩黏性流体的控制方程由雷诺平均Navier-Stokes方程和连续性方程描述．为了封闭这组方程，采用了k-ε湍流模型，对于管内的充分发展湍流，无论是稳态还是非稳态的流固耦合问题，k-ε湍流模型都是很适用的[11]．

2 STAR-CCM+模拟要点

2.1 数值模拟要点

破损船舶进水过程时域模型的建立基于两个条件：一是舱室流量平衡，即与舱室相连的所有破口的流量之和与舱室内进水增量相等；二是破口处满足伯努利方程[12]．

计算域模型的网格划分方法和网格质量对数值求解的计算精度及模拟结果都具有非常大的影响，采用合理的网格控制参数和局部区域(关键主流区域)网格细化控制方法进行网格划分，对减少网格数量、提高计算精度和求解效率具有非常重要的作用[13]．

2.2 边界条件设置

各个舱室均有通风口，不考虑封闭空间气体压强因素，划分区域为船体(HULL)、进水舱室(HULL INSIDE)、随体区域(OVERSET ZONE)、变形区域(TANK)．设置气液交界面(WAVE)、破口(INLET A)，标准模型顶部设置压力出口(INLET B)、混合壁面(INSIDE WALL)，混合壁面切向速度固定，剪应力无滑移，壁面规整平滑．

3 标准模型破舱模拟研究

3.1 标准模型数据与模型尺寸

本文首先采用Ruponen的实验模型[14]．模型来源于阿尔托大学和NAPA：ITTC基准研究的洪水模型测试(ITTC-Box)，收录于船舶与海洋运载器国际稳性会议(STAB)文集，STAB是全球稳定性和一般安全领域最具代表性的专业会议．模型主要参数见表1，箱型驳船标准模型如图1所示，箱型驳船外形尺寸如图2所示．

表1 标准模型主要参数

图1 箱型驳船

图2 箱型驳船外形尺寸

3.2 收敛性分析

STAR-CCM+支持对几何模型进行网格划分，网格包括四面体网格、多面体网格等，与FLUENT相比，支持网格模型较多．此外，STAR-CCM+还支持从其他第三方网格软件(ICEM CFD、Pointwise、Gridpro、Hypermesh等)直接导入体网格进行计算．

破舱进水模拟属于模拟外部流，因此选择切割体网格．切割体网格单元生成器能提供稳定且有效的方法，针对简单和复杂的网格生成问题生成高质量的网格，同时适用于基于零部件的网格化(PBM)和基于区域的网格化(RBM)．除了模型本身质量的影响，面网格质量最小值和时间步长的设置也大大影响了计算精度和计算时间．

为了得到不同网格和不同时间步长对模拟结果的影响，选取箱型驳船进水舱室(HULL INSIDE)进行收敛性分析．

进水舱室尺寸如图3所示，不同网格比较情况如表2所示．由表2可以看出，面网格质量最小值为0.05时，划分网格质量最好，网格数量也满足计算要求．

图3 进水舱室外形尺寸

表2 不同网格对比表

库朗数能够将STAR-CCM+中的网格和时间步长联系起来．其公式为

C=vΔt/Δx

(6)

式中：v为流体流速；Δt为用于VOF方法计算的时间步长；Δx为网格间距；库朗数C<0.3．残差收敛较好的可以适当放大库朗数．

面网格质量最小值均为0.05，采用不同时间步长，模拟终止条件为30 s，则不同时间步长的模拟时间如表3所示．

表3 不同时间步长对比表

残差分析是考察模型假设合理性的方法，通过对残差及残差图的分析，可以判断选取的时间步长是否满足计算需要．

当时间步长设置为0.30 s时，湍动能随着计算的增加呈上升趋势，收敛较差．时间步长设置为0.02、0.05 s时，湍流耗散率和湍动能残差基本相同，都较为平稳，略有向下趋势，因此判断收敛较好，但时间步长设置为0.02 s时，模拟时间较长．结合表3和图4可以看出，时间步长设置为0.05 s，既能保证计算精度，同时模拟时间较少．

3.3 箱型驳船网格建立

由于破舱进水模拟中存在自由液面，为保证计算精度，本文利用STAR-CCM+切割体网格生成器对实船模型进行网格划分，对随体区域和变形区域建立重叠网格．进水舱室及破口尺寸见表4．

(a) 时间步长0.30 s

(b) 时间步长0.05 s

(c) 时间步长0.02 s

图4 残差分析

Fig.4 Residuals analysis

表4 进水舱室及破口尺寸

选取表面重构，执行曲率细化和接近值细化．自动表面修复设置面网格质量最小值为0.05；选择切割体网格单元生成器和棱柱层网格生成器，延伸函数为几何级数，填隙百分比为25%，最小厚度百分比为10%，层减少百分比为50%．面网格增长率1.3，自动修复最小接近值0.01．网格总数量为2 526 717，其中质量较好的网格数量为2 507 240，占99.229%，如图5所示．进水舱室网格划分，设置表面重构，选取切割体网格单元生成器和棱柱层网格生成器．共生成网格684 196，其中质量较好的网格数量为675 567，占98.739%，如图6所示．变形区域网格划分为877 160，其中质量较好的网格数量为876 960，占99.977%，如图7所示．随体区域网格划分为965 361，其中质量较好的网格数量为954 713，占98.897%，如图8所示．

选取常用k-ε湍流模型六自由度求解器，由于不考虑波浪在静水状态下进行模拟，所以VOF波设置为静水VOF波．实际环境中，船底部为液体，上部为气体，所以为气液两相模拟．选取欧拉多相流模型，水设置为恒密度，空气部分也设置为恒密度，动力黏度和密度均设置为常数．在STAR-CCM+中设置物理模型，模型参数选择初始条件水和气体复合湍流强度0.01，湍流速度1 m/s．设置求解器参数：隐式不定常中时间步长设置为0.01 s，默认设置负载平衡、k-ε湍流、分离VOF模型等选项．

图5 整体网格

图6 舱室网格

图7 船体外侧变形区域网格划分

图8 随体区域网格加密

3.4 标准模型进水过程及结果分析

利用STAR-CCM+模拟结果与实验结果对比，结果高度吻合，数据基本准确．箱型驳船进水过程如图9所示．0.40 s时，可以观察到破口进水开始，进舱水在破口处形成水柱冲击舱壁．6.60 s时，破口舱室的进水已经流入中间舱室，中部舱室进水，部分水流入右侧舱室．10.50 s时，3个舱室均有进舱水，舱室内液面逐渐升高．29.95 s时，3个舱室达到水满状态．

监视破口进水舱室内的自由液面高度，绘制液面高度随时间变化曲线，并与文献[8，14]结果进行比较，如图10所示．

(a) t=0.40 s

(b) t=6.60 s

(c) t=10.50 s

(d) t=29.95 s

图9 箱型驳船模拟结果

Fig.9 Simulation results of box barge

本文STAR-CCM+模拟结果与Ruponen的实验结果[14]、FLUENT模拟结果对比[8]，R21S舱在t=0～12.00 s时，由于模拟开始时部分进舱水喷射在R21S舱壁上产生溅射，液面高度较实验结果与文献结果略高；t=12.00 s以后，模拟结果逐渐趋于平稳，与实验结果、文献结果高度吻合．由于左侧R21S舱一部分进舱水通过门直接喷射到R21舱室中，使监测的R21舱平均液面高度增加，但随着舱室内进水量的增加，液面升高，误差逐渐减小，最后监测结果逐渐稳定，与文献结果相同．R21P舱进水无水柱喷射，进水状态平稳，舱室液面高度变化与实验结果、文献结果对比高度吻合．

(a) R21S

(b) R21

(c) R21P

图10 舱室液面高度随时间变化图

Fig.10 Diagram of cabin level change with time

模拟对比结果验证了该模拟方法可行，利用STAR-CCM+可以模拟船舶破舱进水过程，模拟结果准确可靠且模拟方法较FLUENT模拟更方便，占用系统资源更少．

4 渔政船破舱进水模拟研究

4.1 NAPA建模

本文采用一艘522 kW渔政船的数据，利用NAPA进行船舶建模．渔政船是指在渔业专属水域执行渔政任务，担负海上渔政管理监督、执法的专业船只，也是渔政管理设施建设中重要的设施之一；主要用于渔场巡视并监督、检查渔船执行国家渔业法规的情况，也可兼负渔业生产指挥、发布渔情和气象通报以及海上医疗、海难救助等任务[15]．因此渔政船所航行海域复杂多变，容易发生触礁、剐蹭等事故．机舱破口进水将严重影响船舶航行性能．本文选取机舱进行破舱进水的时域模拟．

船体主要参数如表5所示．

表5 渔政船主要参数

建立渔政船模型，如图11所示，舱室划分如图12所示．

选取机舱修改几何，在右侧做半径0.30 m的圆形破口，舱室顶部做通风口，如图13所示．由于模拟处于静水环境中，渔政船无动力、无航速，水流速度为0 m/s，忽略进舱水晃荡的影响，因此在模拟过程中，不考虑舱室内设备影响，视机舱为空舱．

图11 522 kW渔政船模型

图12 舱室结构划分示意图

图13 机舱几何模型

4.2 STAR-CCM+网格划分

变形区域选取切割体网格单元生成器，打开表面重构和自动表面生成器，面网格增长率1.3，自动修复最小接近值0.01；进水舱室选取切割体网格单元生成器、棱柱层网格生成器，打开表面重构和自动表面生成器，面网格增长率1.3，自动修复最小接近值0.01，棱柱层数1，棱柱层延伸1.5，棱柱层总厚度绝对尺寸0.005 m．考虑运算时间，在重叠域之间建立重叠网格交界面．求解器设置为隐式不定常、六自由度求解器，六自由度运动、负载平衡、分离流、分离VOF波、k-ε湍流、k-ε湍流黏度．

总网格数量为699 301，检查网格质量，质量较好的网格达到98.832%．

机舱网格划分总网格数量为707 564，其中质量较好的为699 301，占98.832%，如图14所示；破口处网格加密如图15所示．

图14 机舱网格划分

图15 机舱破口网格细节

船体网格划分如图16所示，划分总网格数量为373 556，其中质量较好的为350 977，占93.956%．

图16 船体网格划分

变形区域网格划分如图17所示，划分总数量为795 358，其中质量较好的网格数量为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，重叠网格区域基本尺寸为56 m．设定0.05 s为时间步长，对于边界条件的设置，顶部边界设置为TOP，底部边界设置为BOTTOM．

图17 变形区域网格加密

DFBI设置六自由度体，启用模型VOF波，考虑重力因素，两层全y+壁面处理，精确壁面距离．选择可实现的k-ε两层模型，常用k-ε湍流模型，雷诺平均Navier-Stokes模型，三维隐式不定常模型．选取欧拉多相流模型，水设置为恒密度，空气部分也设置为恒密度，动力黏度和密度均设置为常数．右边界破口设置为速度入口(INLET A)，顶部通风口设置为压力出口(INLET B)．监视器监视模型六自由度的变化，并生成变化曲线．具体采用STAR-CCM+中的混合壁面函数功能，指定速度入口的函数值、压力出口的函数值，控制计算中边界处多相流分布和速度、压力分布等．

4.3 结果分析讨论

与标准模型实验结果进行对比，已经证明STAR-CCM+在处理船舶稳性研究中可以得到较为可信的结果．采用实船数据利用NAPA进行建模得到实船模型，在破舱进水的过程中，不考虑进舱水的晃荡影响．进水过程瞬时变化如图18所示．

模拟过程中可以观察到在瞬时过程的进水前期，t=8.20 s时，破口处形成水柱喷射，冲击在舱壁上，机舱内产生自由液面，并在船舱底部形成两处涡流，舱内液面逐渐升高；t=80.00 s时，舱内液面高度超过破口高度，舱内进水不形成水柱喷射，但进舱水仍然有涡流．

监测破口处流量随时间变化，生成破口流量变化曲线，如图19所示．可以看出在t=0.07 s时，破口流量为655.81 kg/s．在t=0.07～55.00 s时，由于船体横摇影响，破口流量也产生规则波动，随着进舱水的增加，横摇逐渐趋于稳定，破口流量也逐渐趋于平稳．t=100.00 s时升沉运动持续增大，破口到水面垂直距离增加，破口流量也随之增加．对破口流量曲线进行积分，使用科研绘图工具GraphPad Prism求出曲线下面积，机舱浸满水时，舱内水有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所示，船身向进水一侧移动，横荡距离增加，船体形成较大的横倾角．

第1次横倾角最大值出现在t=0.66 s，偏移角度为1.29°，是舱室浸满水时的瞬时横倾角的16倍，横摇周期约为5.080 s，比设计计算的横摇周期5.199 s略小．

随着进水时间增加，横倾角逐渐减小，舱室内液面升高，液面趋于平稳，机舱进水增加，船舶下沉，艏摇加快，机舱进水状态趋于平稳．由于进舱水影响，纵摇逐渐减小，舱室浸满水约160 s．

(a) 横摇

(b) 横荡

(a) 纵摇

(b) 纵荡

(a) 艏摇

(b) 升沉

5 结 语

本文利用CFD软件STAR-CCM+平台，运用动网格和重叠网格技术对静水状态下的实船模型进行破舱进水模拟．为了保证模拟方法的准确性，先对网格和时间步长收敛性进行分析，然后与标准箱型驳船实验结果进行对比，验证了STAR-CCM+应用于模拟船舶破舱进水的可行性，模拟结果准确可靠．STAR-CCM+的六自由度求解器、VOF波形模型、可实现的k-ε两层模型、欧拉多相流模型均适用于船舶实验的模拟．通过模拟实验，可以观察船舱进水的瞬时状态，通过对监视器的六自由度变化曲线进行分析，能够得到船舶破损时的稳态瞬时数据，依靠软件强大的后处理功能，可将模拟结果生成图表或视频，转为可视化数据．

在确定了一艘渔政船破损位置、尺寸等因素的情况下，对渔政船进行破舱进水时域模拟得到渔政船每一时间步长破损舱室进水位置云图和破舱进水破口流量及六自由度随时间变化曲线．刚开始进水时船舶产生较大横摇角，破口流量规则波动，随着进水量的增加船体下沉，横摇减弱并逐渐稳定，破口流量增大，纵摇随进水量增加逐渐减弱，模拟在舱室浸满水后停止．

本文实验结果可为研究复杂工况如多个舱室破损进水、大风大浪环境下的船舶破舱进水模拟提供基础数据，并为研究船舶破舱进水提供一种新思路和方法．其分析结果可为船体优化、海难船舶救援提供参考．