Maxsurf Motions模块在船舶耐波性上的应用与验证

2023-03-15张大朋严谨赵博文朱克强

科学技术与工程 2023年4期

张大朋，严谨，赵博文，朱克强

(1.广东海洋大学船舶与海运学院，湛江 524088；2.浙江大学海洋学院，舟山 316021; 3.宁波大学海运学院，宁波 315211)

海浪是作用在海洋运动物体上的最主要的外力[1]。船舶作为海洋中最主要的运动物体，其在海浪中的运动问题，是船舶耐波性与适航性的基础。准确预报和正确了解船舶在风浪中的水动力性能，有利于设计和建造耐波性能优良的船舶，同时也有利于在船舶运营过程中采取适当措施以规避倾覆等海难事故的发生[2]。

随着高性能计算技术的迅速发展，船舶在波浪中的运动理论也越来越成熟[3]。各种理论和计算方法开始应用于船舶运动问题的研究，包括切片理论、高速细长体理论、三维频域理论和三维时域理论、非线性理论等[4-6]。船舶在波浪上运动的理论和计算方法从二维切片近似方法发展至三维边界元法、从零航速发展至有航速、从频域法发展至时域法，取得了一系列重要的进展[7-10]。这些方法都是应用势流理论，并且自由面条件进行了线性化近似。大多数方法都是应用面元法进行的。从频域计算转向时域计算是船舶运动的线性理论的另一个重要发展趋势，这也反映了工程设计中对一些非线性问题的需求。时域计算方法较频域法在处理一些带非线性问题上有明显的适用性。目前时域方法中一种是应用频域计算的结果，通过相应的变换确定辐射力的时延函数和波浪干扰力的脉冲响应函数。应用这些函数，通过对运动时历或来波的卷积计算相应的辐射力或干扰力。最后通过对时域微分方程积分解出运动。

计算机的高速发展为船舶在波浪中的运动响应预报问题提供了一条新的途径，计算流体动力学(computational fluid dynamics，CFD)方法[11-13]。应用CFD技术进行船舶与海洋工程结构物流动数值模拟和水动力性能预报的优越性和应用前景日益显现。耐波性中涉及的甲板上浪、船体砰击和载液船体的液舱晃荡等强非线性问题己经可以通过将求解非定常RANS方程进行研究，在一定程度上可以解决部分工程实际问题[14]。近年来，基于黏流理论的CFD软件如FLUENT、STAR-CCM+、OpenFOAM等得到了快速地发展，已经被广泛应用于船舶耐波性及水动力计算的各类问题，为众多学者的研究提供了强有力的计算工具。然而，目前来看理论分析和数值计算仍不能完全替代模型实验，最后通常还是需要用模型实验的结果来验证理论计算是否正确。此外，在非线性的运动计算分析、结果预报中还是有很多问题存在，需要进一步改进的理论部分还有很多。

目前，应用势流理论计算方法快速地预估某一船舶在指定海况中的运动特征，在设计方案的比较阶段用理论与计算相结合的手段对其耐波性能进行比较和选择，已经成为船舶设计流程中的常规环节[15]。Maxsurf Motions是一款基于势流理论的船体耐波性预测程序，目前在船体耐波性计算的研究上已得到较多的应用，但其计算精度仍需要进一步讨论。采用Maxsurf Motions程序对4种国际标准船模进行耐波性计算，验证Maxsurf Motions中的切片理论和面元法的可靠性，探讨Motions程序在船舶耐波性分析中的可行性，以期对Maxsurf Motions的应用起到一定的借鉴和指导作用。

1 基本理论

Maxsurf Motions程序有两种方法用于计算船舶在波浪中的运动响应，分别是切片理论和面元法。切片理论用于计算船舶的垂荡和纵摇响应。横摇响应通过线性横摇阻尼理论计算。面元法是一种一阶衍射/辐射水动力分析方法，其中使用了基于面元的恒定边界元法(boundary element method，BEM)，能够计算流体动力附加质量和阻尼、一阶波浪激振力和力矩、平均漂移力和力矩以及船舶湿表面上的压力。

1.1 切片理论

使用切片理论计算船舶垂荡和纵摇的耦合运动时，包含以下基本假设：船长远大于船宽和吃水，船宽远小于波长；船体为刚性，舷侧直立；航速适中且保持不变；相比于波幅，船体运动很小且呈线性；水深远大于波长；船体对入射波没有影响(Froude-Kriloff假设)。

船舶的垂荡、纵摇和横摇运动本质上都是振荡的，这是由于这些运动中均包含由浮力变化产生的恢复力。船舶在波浪作用下的运动可视为阻尼-弹簧-质量系统。

垂荡和纵摇的耦合运动中，垂荡方程为

C35η5=F3eiωt

(1)

纵摇方程为

(2)

当船舶在波浪中航行时，除了会产生与航行速度和水深相关的船行波外，还会产生另一种与入射波浪诱导船体运动相关的波。这两种波均会消耗船体的能量且互相叠加，因此船舶在波浪中航行时所消耗的能量要高于静水中，其中额外的能量损失即为所谓的波浪增阻。Maxsurf Motions程序的波浪增阻计算方法分别为：Gerritsma&Beukelman法、Salvesen法和Havelock法。

Gerritsma&Beukelman法所描述的波浪增阻RAW与船舶的相对垂直速度有关，计算公式为

(3)

Salvesen法中，波浪增阻的计算公式为

(4)

ib33(ξ)]}dξ

(5)

(6)

(7)

Havelock法中，波浪增阻的计算公式为

(8)

式(8)中：ε为运动与相应激振力或力矩的相位差；下角标3表示垂荡；下角标5表示纵摇运动。

1.2 面元法

面元法将物体表面离散，用一个平面或曲面代替原来的物面(称为面元)，在该面元上布置流动的奇点如源、涡、偶极子及其组合进行求解。面元法中引入一个旋转中心，如图1所示。物体遭遇波浪后通过求解绕旋转中心上的简谐运动，进而求解船体在波浪中的运动问题。

图1 面元法旋转中心

假定波高和陡度较小，同时流体不可压缩、无黏、无旋，因此可以使用线性波理论。流场用速度势的梯度来描述，速度势由Laplace方程控制，同时应满足适当的边界条件。对于简谐运动，速度势可表示为

Φ(x;t)=Re[Φ(x)e-iωt]

(9)

式(9)中：Re为雷诺数。

在整个流体运动域内，必须满足Laplace方程：

∇2Φ(x)=0,x∈Ω

(10)

自由水面边界条件为

(11)

物面边界条件为

(12)

海底边界条件为

(13)

式中：∇为拉普拉斯算子；Ω为实数集；g为重力加速度；n为海底的外法线方向；Un为速度；h为水深。

在线性理论的框架内，速度势Φ被细分为入射势ΦI、衍射势ΦD和辐射势Φj。

(14)

式(14)中:xj为第j个运动模式的复振幅。

入射波的速度势定义为

(15)

式(15)中：h为水深；β为波浪入射角。

假设浮体是刚性的，并且在平静的水中处于稳定的平衡状态。考虑作用在浮体表面上的水动力，平台运动方程为

(Cij+Kij)]=Fi

(16)

2 船舶耐波性计算与验证

2.1 Series-60船舶

Series-60是一艘不带球鼻艏的货船船型，是国际公认的用于验证的标准船模，有较丰富的实验数据和文献资料。Series-60船模的方形系数Cb=0.68，缩尺比λ=3.93。几何模型如图2所示。几何参数如表1所示。

图2 Series-60船舶几何模型

表1 Series-60船舶几何参数

验证工况为迎浪，分别计算两个航速Fn=0.198和Fn=0.297。将垂荡和纵摇运动的响应幅值算子(response amplitude operators，RAOs)与实验和SWAN(simulating waves near shore)程序形成对比，绘制在图3和图4中。其中，实验结果是由Gerritsma等[16]开展的实验，即G-B实验，MIT-SWAN是Sclavounos等[17]开发的一款船舶非线性时域分析程序。

图3 Fn=0.198的RAOs曲线

图4 Fn=0.297的RAOs曲线

当Fn=0.198时，垂荡响应中，Motions计算的垂荡幅值响应算子(response amplitude operator，RAO)精度不如MIT-SWAN，但成功预测到了处于共振峰值时的遭遇频率范围，Motions对峰值的预测误差约15%，SWAN的误差约2%。纵摇响应中，Motions对峰值的预测误差约-16%，SWAN的误差约8%。当Fn=0.297时，Motions对纵摇响应的计算精度要优于垂荡响应。Motions能合理地预测Series-60船舶的运动响应。当Series-60船处于垂荡和纵摇的共振频率范围时，此时船体运动最复杂，Motions的计算误差较大。当遭遇频率较低或较高时，Motions的计算结果与实验和SWAN程序较吻合。Motions的优势在于能够在短时间内(约几十秒)预测多个频率的RAOs，而相同的工况SWAN程序需要超过8 h的时间。

2.2 SA船舶

SA(SA指在特殊区域special area内航行的船舶)船舶是一艘不带有球鼻艏和方艉的快速货船，船型瘦长，符合切片理论的假设。SA船舶的几何模型如图5所示。几何参数如表2所示。

表2 SA船舶几何参数

图5 SA船舶几何模型

验证工况为迎浪，计算4个航速Fn=0.15、0.2、0.25、0.3，对应5.8、7.73、9.67、11.6 m/s。验证对象主要是垂荡和纵摇运动的RAOs曲线、运动相位滞后以及附加阻力，与Journee[18]在2001年公布的实验数据以及Geritsma[19]的实验数据形成对比。相位滞后是指垂荡和纵摇运动与相应激振力/力矩的相位差。

通过图6和图7可以看出，基于切片理论的SA船舶的运动响应与相位滞后总体上与实验基本吻合，但在某些频率点出现较大的误差。Fn=0.15时，纵摇运动的RAOs曲线比垂荡运动更贴合实验，而纵摇运动的相位滞后却与实验差距较大，尤其是当船长波长比大于1(也就是波长小于船长)时，实验的纵摇运动相位滞后与计算值相差一个负号。同样的情况发生在Fn=0.3。由于实验数据不足，实验结果存在相当大的分散性，且在两种航速下实验都没有很好地捕捉到响应的峰值，因此很难说明计算结果与实验结果不吻合的原因是切片理论的精度不准确。然而大部分计算结果显示，当波长较大时，SA船舶的垂荡和纵摇耦合运动响应与实验吻合良好，一定程度上可以证明对于SA船舶的运动计算，Motions程序在大部分工况下具有较强的精度。

L为船长；λ为缩尺比

图7 Fn=0.3的SA船体运动响应

Raw为附加阻力；ρ为水的密度；g为重力；ζ为波高；B为船宽；L为船舶水线长

观察图8可知，无论是低频段还是高频段，Havelock法的计算结果均比较符合实验情况，尤其是在航速较大的情况下，高频段的计算结果与实验数据十分贴合。Salvesen法明显不适用于SA船舶的附加阻力计算，无论是曲线趋势还是数据点，均与实验相差甚远。

值得注意的是，无论是哪一种计算方法，都无法做到在整个波频内与实验完全吻合。这其中有尺度效应的影响。SA船舶是实尺度，而实验是模型尺寸。在拖曳水池中精准地测量出船舶的附加阻力是十分困难的，且从模型尺度到实尺度的换算存在修正系数。SA船舶的修正系数到目前为止仍是未知数。综合以上原因，可以认为Motions模块中的附加阻力计算方法具备一定的计算精度，能够应用于船舶早期设计阶段和附加阻力的估算。

2.3 Wigley船舶

Wigley船型满足切片理论和Michelle兴波阻力积分公式的假设，型值由式(19)可得

(19)

式(19)中：Lpp为垂线间长；B为最大船宽；d为吃水；x、y、z为空间三维点坐标。

d/Lpp=0.062 5；B/Lpp=0.1。几何模型和参数分别如图9和表3所示。

表3 Wigley船舶几何参数

图9 Wigley船舶几何模型

验证工况为迎浪，计算3个航速Fn=0.2、0.3、0.4，对应1.09、1.64、2.19 m/s，验证对象主要是附加质量、阻尼系数和波浪力。

图10～图12为基于Motions模块得到的不同航速下Wigley船体的附加质量系数随遭遇频率的变化，并与Journée[20]的实验结果进行对比。Wigley船体在波浪中运动时，会促使周围流体的速度发生变化。由于流体存在惯性，流体将阻碍船体运动并反作用于船体上，其作用力称为附加惯性力，把该附加惯性力假想为在物体上附加了一部分质量，这部分质量称为附加质量。附加质量与物体本身的形状及运动方向有关。

由图10～图12可知，不同航速下Journee的实验结果均匀地分散在Motions曲线附近，说明基于切片理论的Motions模块在不同航速下计算船舶附加质量时均能得到与实验较为接近的结果。其中，中高频段的吻合度较好，部分低频段的吻合度较差，说明Motions模块对中高频的计算精度较高。3种航速下Motions计算的附加质量系数差距不大，说明航速效应对附加质量的计算影响较低。4个附加质量系数中，A33和在3个航速下的计算数据几乎一致，说明因垂荡产生的垂荡附加质量不随航速的变化而发生明显变化。A35和A53的结果相差一个负号，说明因纵摇产生的垂荡惯性力和因垂荡产生的纵摇惯性力大小相同，方向相反。A55是因纵摇产生的纵摇附加质量，航速的增大会使纵摇附加质量略微增大。4个附加质量系数随遭遇频率的变化是一致的：在中低频段，附加质量系数随遭遇频率的增大快速减小或增大，在高频段，附加质量系数曲线几乎趋于一条直线，说明此时的遭遇频率基本不会对附加质量造成影响。

V为船舶排水体积

图11 Fn=0.3时Wigley附加质量随遭遇频率的变化

图12 Fn=0.4时Wigley附加质量随遭遇频率的变化

图13～图15是基于Motions模块得到的不同航速下Wigley船体的阻尼系数随遭遇频率的变化，并与Journee的实验结果进行对比。在势流理论框架下已经假设流体无旋、无黏、不可压，因此这里的阻尼不是粘性阻尼，而是物体在波浪运动中所做的功。由于物体在含有自由表面的流体中运动，势必会兴波，而兴起波浪必然会带走能量，这个能量就是物体运动做的功，体现为物体受到的阻力，大小与物体运动的速度相关。

图13 Fn=0.2时Wigley阻尼系数随遭遇频率的变化

图14 Fn=0.3时Wigley阻尼系数随遭遇频率的变化

图15 Fn=0.4时Wigley阻尼系数随遭遇频率的变化

相比于附加质量，Motions对阻尼系数的计算精度明显降低。除了B33之外，其余3个阻尼系数均与计算值出现了较大的差别。对于B35和B55，中低频段的计算值偏大，高频段的计算值偏小，B53的现象与B35相反。这3个系数均与纵摇模态相关。而且也看出，随着航速的增大，阻尼系数与实验值的偏离度增大，航速效应对阻尼系数的影响要大于对附加质量的影响，尤其是对纵摇模态的阻尼系数的影响最为明显。

2.4 KVLCC2船舶

KVLCC2是一艘由韩国KRISO设计的现代油船，具有详细的试验数据。KVLCC2船模带有球鼻艏和U形船尾框架线，属于肥大型船。通过KVLCC2船舶验证面元法对斜浪中的运动响应与漂移力的计算情况。KVLCC2实物模型如图16[21]所示，几何参数如表4所示。

表4 KVLCC2船舶几何参数

图16 KVLCC2实物模型[21]

验证的工况共有2个，分别是迎浪和艏斜浪135°，航速为0。与切片方法不同的是，面元法的分析需要网格的支持。根据KVLCC2的船体特征，船首、球鼻艏和船尾的曲率较大，应该生成较小尺寸的网格以捕捉到船体曲面，平行中体处的曲率变化较小，可以将网格尺寸适当调大。

KVLCC2各曲面的网格尺寸如图17所示。网格划分结果如图18所示。图19为面元法对180°浪向下的KVLCC2运动响应计算结果，并与CFD结果以及三维势流理论结果[22]形成对比。对于垂荡运动，面元法的计算结果与CFD和三维势流理论吻合良好。误差较大的地方发生在1<λ/L<2的区间，也就是波长为1倍和2倍船长的工况。对于纵摇运动，面元法结果的总体趋势与CFD及三维势流理论一致。在长波工况下误差较大。

Mesh用于指定表面是否自动网格化，勾选表示该表面进行网格的划分；Type表示网格类型，目前只有Tri.一种，即三角形网格；Min.edge length用于指定三角形网格的最小边长；Max.edge length用于指定三角形网格的最大边长；Patch egde angle用于指定网格的接触角度；Node merge tolerance指相对容差，用于指示两个相邻的网格节点是否由于过于靠近而合并，该值应在0.001～0.05；Geometric tolerance指几何公差，用于指定网格边可以偏离与其关联的曲面的距离，该值应在0.001～0.05，较低的值将在高曲面曲率区域创建较小的三角形；Surface Mesh Parameters为表面网格参数；Surface Name为表面名称；Characteristic Leng为特征长度

Zero pt.为零点；AP为艉垂线；FP为艏垂线；Baseline为基线；CF为漂心；CG为重心；CB为浮心

图19 180°浪向下的响应幅值算子

对于在波浪中做摇荡运动的船舶来说，作用在KVLCC2的外力包括入射力、绕射力、辐射力和静水回复力。CFD结果基于黏流理论，考虑了黏性阻尼的影响。三维势流理论在求解边界积分方程时，用无航速的脉动源格林函数代替难以计算的移动源格林函数，在理论上与有航速船舶的水动力问题并不完全吻合。因此对前进速度效应的准确预测是格林函数的主要局限。面元法计算船舶在波浪中的运动响应时，将船体周围流场速度势分解为波浪入射势、自身绕射势、相邻船体产生的绕射势、自身辐射势和相邻船体运动产生的辐射势，通过源汇分布法进行求解。相比CFD方法，面元法能够大幅减少计算成本；相比三维势流理论，面元法通过对网格单元进行压力积分得到二阶漂移力和力矩，计算内容更加丰富。