规则和随机横浪中破损船舶运动响应研究*

吴 浩 高志亮 苏 焱

(武汉理工大学高性能船舶技术教育部重点实验室1) 武汉 430063) (武汉理工大学交通学院2) 武汉 430063)

0 引 言

随着船舶尺度的增大，船舶发生碰撞或触礁的概率增加.2018年初，油船“桑吉”号与散货船“长峰水晶”号发生碰撞，“桑吉”号油船燃爆失火并最终沉没；同年，挪威护卫舰“海尔格·英斯塔”号与油轮“索拉”号相撞，事故发生后护卫舰虽采用了海岸搁浅、损管作业，以及钢索固定等补救措施，最终仍然沉没.对大量的船舶碰撞事故事后分析表明，船舶破舱进水后正确的应急决策能将人员伤亡、财产损失，以及环境污染降至最低.破损船舶在海浪中的运动行为十分复杂，对其准确判断决定了应急决策系统的正确性.

船模试验方法难以对破损船舶运动的影响因素进行系统分析，不利于深入探究破损船舶与流体相互作用的动力学机理.数值模拟方法为上述问题研究提供了可行的途径.早期，科研工作者将传统耐波性理论与基于准静态假设的破舱进水模型相结合，开展船舶在波浪中的运动研究[1-3].该方法采用势流理论计算船舶波浪力，基于伯努利方程计算破口处的流量，并假设舱室内水面始终保持水平.上述方法计算简单，但忽略了水体流经破口时的动态行为及舱室内部自由液面运动等瞬时因素，因此无法准确预报水流经破口和在舱内运动时的动态作用力.采用计算流体动力学(computational fluid dynamics，CFD)方法则可合理地计入破舱进水问题中流体的动态效应.

近年来，学者们采用基于求解纳维-斯托克斯方程的CFD方法开展了破损船舶在波浪中的运动研究.Gao等[4-5]采用势流理论求解船波干扰问题，并与CFD方法求解的破舱进水问题进行耦合，分析了破损滚装客船在横浪中的运动和受力.Sadat等[6]采用CFD方法模拟了破损客船在规则横浪中的横摇运动.Haro等[7]采用类似方法对同一艘船在破损状态下规则迎浪航行时的快速性问题进行了研究.上述研究分析了破损船舶在规则波中的运动，而破损事故通常在真实海况中发生，因此，开展破损船舶在随机波浪中的运动预报更加符合实际情况，同时也可为事故应急决策系统提供有用参考[8-9].

本文将基于求解雷诺平均纳维-斯托克斯(reynolds averaged navier stokes，RANS)方程的CFD方法与谱分析方法相结合，研究了破损舰船DTMB 5415在规则横浪和随机横浪中的运动响应.通过CFD方法建立数值波浪水池，进行了数值水池造波精度的验证.采用所建立的数值水池模拟了破损船舶在规则横浪中的运动，并将船舶横摇频率响应函数与试验结果进行对比验证.在此基础上，采用频谱分析方法预报了该艘破损船在四到六级海况的随机横浪中横摇运动响应.

1 数学模型

本文采用CFD商用软件STAR-CCM+对船舶运动流场进行模拟，将空气和水看作不可压缩流体，流体运动由连续性方程和RANS方程控制.在笛卡尔坐标系下，方程描述为

(1)

(2)

式中：t为时间；xi为坐标分量(i=1,2,3)；ui为xi方向的速度分量；ρ为流体密度；μ为黏性；I为单位矩阵；k为湍动能；p为压力；F为体积力；Si为自定义源项.

采用SSTk-ω湍流模型，方程描述如下.

(3)

(4)

式中：ω为耗散率；Γk和Γω为k和ω的有效扩散项；Pk和Pω为k和ω的产生项；fβ*为自由剪切修正系数；fβ为漩涡拉伸修正系数；ω1和k1为环境湍流值.上述项的计算和系数取值见STAR-CCM+用户手册.

本文采用流体体积(volume of fluid，VOF)方法对水和空气交界面进行捕捉，该方法通过计算网格单元中流体所占的体积分数(α)来捕捉流体的交界面，α满足以下方程.

(5)

船舶在流体中的运动满足

(6)

(7)

式中：vi和ωi分别为船体运动线速度和角速度；m为船体质量；Ji为船体转动惯量；fi和Mi分别为作用于船体上的合力和合力矩.

STAR-CCM+软件基于有限体积法对方程(1)～(5)进行离散，其中时间项采用二阶隐式格式，方程(5)对流项采用二阶高分辨率交界面捕捉(high-resolution interface capturing，HRIC)格式计算，其余方程中的对流项采用二阶迎风格式进行计算，扩散项则采用中心差分格式进行计算.采用半隐式速度压力耦合方程(semi-implicit method for pressure linked equation，SIMPLE)方法对流场的速度和压力进行耦合求解.

2 数值预报及分析

2.1 研究对象

Begovic等[10]对DTMB 5415船模破损情况下在波浪中的运动开展了试验研究.本文的研究对象为该船模(见图1)，其主要船型参数和破损舱室参数见表1～2.破损舱室为位于船中附近的两个舱室，两舱之间舱壁位于破口中间.数值模拟中，在破损舱室顶部布置了与试验相同的两个通气管，以便船舶在运动过程中舱室内的空气可以顺利排出.

图1 船型示意图

表1 DTMB 5415船型参数

表2 DTMB 5415破损舱室信息

2.2 数值波浪水池建立及验证

STAR-CCM+软件提供了基于Stokes波理论的边界造波方法，一阶Stokes波的流场水平速度(u1)、垂向速度(u3)和波面抬高(η)方程为

(8)

(9)

η=acos(Kx1-ω0t)

(10)

式中：a为波幅；ω0为波浪频率；K为波数；x1方向为波浪传播方向；x3为距静水面的垂直距离；d为水深.

采用上述边界造波方法生成一阶Stocks波，波长λ取3.2 m，波高取波长的1/50.图2为数值波浪水池的示意图，其中计算域长8λ、宽0.5 m，水深2.15 m，该水深与文献[10-11]试验水池深度一致，水面以上部分高2 m.水池左侧设为速度入口边界，右侧设为压力出口边界，前后两侧设为对称面边界，上下两侧设为不可滑移固壁边界.

图2 数值水池示意图

为了减少波浪在水池末端反射影响，在方程(2)中施加阻尼源项，源项作用范围为距离计算域末端2λ区域，源项表达式为

(11)

(12)

式中：xsd和xed分别为阻尼区的起始点和结束点.

数值水池的波面附近网格划分见图3，在波面区域内沿波高和波长方向进行网格加密；波高范围内网格均匀分布，并向水池顶部和底部逐渐变疏；波长方向水池前端6λ区域采用均匀网格布置形式，水池后端2λ区域网格逐渐变疏.此处模拟二维波浪传播，故在计算域宽度方向只设置了两层网格.

图3 数值水池网格划分示意图

首先对波浪传播模拟进行了网格依赖性分析，通过对波高和波长方向的网格调整生成了三种疏密不同的网格，分别记为Mesh1、Mesh2和Mesh3，计算时间步长取波浪周期的1/1 000，在距入流口4λ处记录波面高度，表3为计算结果与理论结果的比较.可以看出，随着网格加密，计算所得的波高和波周期误差降低.当网格达到Mesh2的分辨率时，误差已小于1%，故后面计算根据Mesh2的网格布置形式来生成网格.

表3 网格依赖性分析结果(距入口边界4λ处测点)

采用Mesh2的网格布置形式，通过选取三种时间步长，分别记为Timestep1、Timestep2和Timestep3，对计算模型进行了时间步长依赖性测试，时间步长依赖性分析结果见表4.随着时间步长变小，计算所得的波高和波周期误差降低；当时间步长小于波周期的1/800时，计算误差小于1%，所以后面计算时间步长取波周期的1/800.

表4 时间步长依赖性分析结果(距入口边界4λ处测点)

图4为时间步长为波周期的1/800时，距入口4λ处测得的波高历时曲线与理论值比较.图5为波传播20个周期后数值水池波面与理论值比较.可以看出，在工作区内，数值水池造波结果与理论结果基本一致；在消波区内，波浪朝着出口边界方向逐渐衰减为零.

图4 距入口4λ处计算波高与理论波高比较(Mesh2，Timestep2)

图5 波传播20周期时计算波高与理论波高比较(Mesh2，Timestep2)

2.3 破损船舶在规则横浪中运动模拟

破损船舶在规则横浪中运动模拟的计算区域及边界设置见图6，船体左右各取4λ，前后各取1倍船长，水深为2.15 m，水面以上高2 m.将整个计算域分为背景区域和重叠区域，在计算中背景区域保持静止，而重叠区域随船体一起运动.水面附近网格生成参考上述数值波浪水池Mesh2网格布置形式，舱室内网格尺寸为0.01 m.船体附近网格采用了局部加密的布置形式，远离该区域网格逐渐变稀疏，见图7.

图6 破损船舶在规则横浪中运动计算区域

图7 船体附近网格布置

在计算域末端2λ区域采用3.2中描述的阻尼消波处理.此外，为了消除波浪在计算域入口边界的二次反射，在方程(2)中源项加入力源项，力源作用区域为计算域前端2λ区域，力源项表达式为

(13)

采用上述模型模拟了破损船舶在不同波长的规则横浪中的运动，波高取波长的1/50，破口面向来波方向.计算中仅考虑船体的横摇、纵摇和垂荡三个自由度的运动，初始时刻舱内水面与静水面持平.表5为计算所得的船舶运动稳定后的横摇幅值与试验结果[12]比较，图8为横摇频率响应函数计算结果与试验结果比较.CFD计算结果与试验结果吻合较好，当波长较小时由于船舶运动幅度较小，数值计算结果与试验结果的相对误差较大，但横摇幅值绝对误差小于1°，说明本文计算所得结果比较可靠.

表5 计算横摇角与试验值比较

图8 横摇频率响应函数比较

图9为波浪周期为1.47 s时船舶的运动历时曲线图，该波浪周期与船舶破损情况下横摇固有周期相近.由图9可知，计算14 s后船舶的横摇运动趋于稳定，其横摇角幅值为16.67°；船舶的纵摇运动幅度较小，其幅值约为0.2°，平均纵倾角约为1.1°；船舶的垂荡运动幅值约为0.03 m.图10为船舶运动不同时刻舱内水体的分布情况，由图10可知，舱室内水体运动较为剧烈，出现自由液面爬升、翻卷和触碰舱顶等现象.

图9 船舶运动历时曲线

图10 舱内水体运动情况

2.4 破损船舶在随机横浪中运动预报

将上述CFD计算所得的船舶在规则横浪中的频率响应函数换算成实船尺度，在已知波浪谱密度的条件下，采用谱分析方法预报该实尺度破损船舶在随机横浪中的横摇运动.波浪谱采用国际拖曳水池会议(international towing tank conference，ITTC)建议的双参数谱，谱密度函数公式为

(14)

式中：H为三一波高；T为波浪特征周期.

表6为谱分析法得到的船舶横摇运动的统计值.由表6可知，对于四级和五级海况，该破损船舶的三一横摇幅值分别为1.47°和4.69°，对应的单幅有义值分别为1.4°和3.01°，均满足北欧合作研究计划建议的军船耐波性衡准指标横摇角单幅有义值小于4°的要求[13]，说明该船出现两舱破损情况时在四级和五级海况下横摇运动量较小；在这两种海况下，船舶的百一横摇幅值分别为2.46°和7.83°，表明当前破损条件下该船基本上不会出现较大的横摇角度和倾覆行为.对于六级海况，该破损船舶的三一横摇幅值为13.55°，对应的单幅有义值为5.06°，该值略大于北欧合作研究计划建议的衡准值；船舶的百一横摇幅值为22.63°，表明在六级海况下该船有可能出现较大的横摇角度，需要防范倾覆行为的发生.此外，船舶破损后在上述三种海况下的横摇运动的平均周期与船舶在完整状态下的横摇运动固有周期9.77 s较为接近，且均大于各自海况对应的波浪特征周期.

表6 破损船舶在随机横浪中的运动响应统计值

3 结 束 语

本文采用基于求解RANS方程的CFD方法对破损舰船DTMB 5415在规则横浪中的运动进行了时域模拟.研究中建立了数值波浪水池并对其进行了网格和时间步长依赖性分析；在此基础上模拟了破损船舶在规则横浪中的运动，计算所得的横摇运动响应与试验结果吻合较好，表明本文采用的计算模型是有效的.对该船破损情况下在规则横浪中发生横摇谐摇时的运动进行了分析，得出其横摇谐摇角幅值约为17°，平均纵倾角约为1°，垂荡运动的幅值约为0.03 m.随后，基于ITTC双参数谱，采用谱分析方法预报了该破损船在随机横浪中的横摇运动统计值.统计结果表明，船舶在四、五和六级海况下，单幅有义值分别约为1°，3°和5°，横摇运动平均幅度较小；在六级海况下，船舶百一横摇幅值约为23°，表明该船出现两舱破损情况时在该海况下会出现较大的横摇运动，需要注意其稳性相关的安全性问题.今后工作将在本研究基础上开展船舶在波浪中发生瞬间破损后的运动行为研究.