基于数值计算的船舶波浪增阻和运动响应预报方法论证

陆泽华， 曹 旭， 李传庆

(上海船舶运输科学研究所有限公司 航运技术与安全国家重点实验室， 上海 200135)

0 引 言

船舶在海上航行时，受风、浪等因素的影响，在功率不变的情况下，除了顺风浪航行，其他状态下的航速都比静水状态下的航速小，该现象称为“船舶失速”。船舶失速研究是目前航运业最关心的课题之一，失速特性是定量分析、优选航运路线的主要依据，失速预报精度直接影响最佳航线的设计质量。准确预报船舶失速情况的关键在于精确预报船舶波浪增阻，目前关于船舶波浪增阻和运动响应的预报暂无标准化方法，通常采用船模试验和数值计算2种方法实现，其中：

1) 数值计算方法主要包括势流二维切片法、势流三维面元法和黏流数值计算法等3种；

2) 船模试验即在波浪中以定常速度拖动船模并测出其受到的平均阻力，该值与在静水中以相同速度拖动船模测得的阻力的差值为波浪增阻。

势流理论分析主要有二维计算和三维计算2种方式，忽略了流体的黏性效应。二维切片理论假定船体细长，可沿船长方向将船体分成若干段，各段的截面形状相同，对于各截面来说，流体流动可近似认为是二维的。但是，船舶实际上是三维物体，三维面元法更适合用来求解其水动力问题，反映其船型特征。三维面元法将控制微分方程转换为边界上的积分方程进行求解,通常不易获得解析解。

黏流数值计算法考虑流体的黏性效应，采用一系列离散点的变量值集合代替时域和体域连续的物理量，并通过建立反映该离散点上的场变量之间关系的代数方程组，求得场变量的近似解。

上述计算方法对是否考虑水黏性的影响及对船体表面采用的离散方法都不统一，且都根据实际情况有一定的假设条件，这些因素均会导致船舶波浪增阻计算结果出现差异。

本文以一艘19 000载重吨多用途船为研究对象，分别采用上述方法计算其波浪增阻和运动响应，结合国际拖曳水池会议(International Towing Tank Conference，ITTC)双参数谱，预报该船在实际海况下的波浪增阻和运动响应，并将所得结果与物模试验结果相对比，为这些方法的实际工程应用提供参考。

1 国内外研究进展

1) 在船模试验研究方面：为提高波浪增阻船模试验预报精度，近年来一些学者对船模试验方法和测量设备进行了改进。例如：UENO[1]开发出了一种新的试验技术，并开展了船舶在斜浪中航行时的波浪增阻试验研究；GUO等[2]对短波增阻进行了详细的试验研究，开展了波浪校核，并采取了消除噪声等措施；封培元等[3]采用自由自航模型对某超大型集装箱船进行了全浪向波浪增阻测试。

2) 在势流理论分析研究方面：SALVESEN等[4]基于细长体和高频低速假设提出了STF(Salvesen-Tuck-Faltinsen)法；LEGOVIC等[5]提出了混合计算方法，短波部分采用Faltinsen渐近公式计算，长波部分采用Salvesen方法计算，过渡部分由二者的计算结果相加得到；戴遗山等[6]提出了由近场公式直接推导远场公式的方法；刘亮[7]采用切片理论和辐射能法计算了运动增阻，采用日本经验公式计算了绕射增阻；李传庆[8]采用STF切片理论计算了船舶运动响应，采用扩展积分方程的方法求解了二维辐射问题；吴思莹等[9]基于势流理论对大型船舶最小装机功率进行了波浪增阻预报，提出采用三维方法进行预报，效率较高且所得结果准确。

3) 在黏流数值计算研究方面：RHEE等[10]采用泊松方程投影法对速度与压力进行耦合，研究了船模波场和阻力；DENG等11]采用ISIS-CFD求解器研究了KVLCC2(KRISO Very Large Crude Oil Carrier 2)船模在波浪中受到的阻力；董志等[12]采用商业软件FLUENT建立了可模拟弱非线性波的数值波浪水池；李纳等[13]基于商业软件STAR-CCM+对远洋船开展了波浪增阻数值计算，并对球鼻艏进行了选型分析。

2 目标船的波浪增阻和运动响应计算

2.1 目标船主尺度参数

本文以19 000载重吨多用途船为研究对象，其主尺度参数见表1。

表1 目标船主尺度参数

2.2 势流二维切片法

本文沿纵向将目标船等间距划分为40个分段，建立二维切片模型，见图1。采用基于二维切片理论的波浪增阻数值计算程序[8]计算目标船在11.5 kn航速下的波浪增阻和运动响应。采用STF切片法计算船舶运动响应；采用辐射能法计算船舶运动增阻，采用FALTINSEN等提出的渐近公式计算波浪反射增阻，二者之和为目标船在波浪中的阻力增加值。

图1 目标船二维切片模型

2.2.1 基于二维切片法的船舶运动响应计算

假设船体响应线性简谐，即不考虑大幅横摇和甲板上浪等非线性问题，船体六自由度耦合运动方程为

(1)

考虑x处剖面的二维辐射问题，假设流体理想、不可压缩且无旋，则辐射势可分解为3种基本运动之和，即

(2)

式(2)中：Φ(y,z,t)为船舶辐射势；φj(y,z)为单位复速度的j模态振荡运动引起的辐射速度势；ω为物体简谐运动频率，在有航速的情况下等于遭遇频率。

各速度势分量满足以下边值条件，即

(3)

式(3)中：Nj为二维单位外法向速度矢量；k为波数;ν=ω2/g,g为重力加速度。

2.2.2 波浪增阻计算

2.2.2.1 基于辐射能法的运动增阻计算

考虑纵向规则来波，船舶在1个遭遇周期Te内受到的辐射能P为

(4)

式(4)中：L为船长；xb为纵向坐标；vZ为相对垂向速度；b33,2D为有航速时的截面阻尼系数。b33,2D和vZ的计算公式为

(5)

(6)

考虑1个遭遇周期内的辐射能量的平均值，船舶运动增阻RaωR为

(7)

2.2.2.2 基于渐近公式法的反射增阻计算

Faltinsen渐近理论考虑波浪遇到船体发生反射(绕射效应)时作用在船体上的力，该方法假设船体水下部分是垂直的，在波浪中不产生振荡运动。当来流方向沿船体水线方向时，短波部分的绕射问题即为入射波作用于有限长垂直壁面上的问题，结合能量动量守恒定律，迎浪情况下船体单位长度受到的反射增阻可表示为

(8)

由此，船舶反射增阻为

(9)

综上，采用二维切片法计算的船舶波浪增阻为

Raω=RaωR+RaωD

(10)

2.3 势流三维面元法

本文将目标船划分为四边形面元，共1 632个(全船)，建立三维面元模型，见图2。利用基于三维格林函数的波浪增阻数值计算程序[14]计算目标船在11.5 kn航速下的波浪增阻和运动响应。采用移动脉动源格林函数法计算有航速船舶的水动力系数，从而得到船舶在1个遭遇周期内的辐射能量，进而采用辐射能法计算船舶在波浪中的辐射增阻。

图2 目标船三维面元模型

2.3.1 基于三维法的船舶运动计算

三维频域内的船舶六自由度运动方程可表示为

(11)

假设流体理想、不可压缩且无旋，引入微幅波假定，即入射波和船舶运动引起的辐射波及由于船体存在的绕射波的波幅均为小量，则一阶速度势Φ(x,y,z,t)满足边值条件

(12)

对于频域问题，Φ(x,y,z,t)=φ(x,y,z)·e-iωt，空间速度势φ(x,y,z)可分解为φ=φI+φR+φD，其中：φI为入射波速度势；φR为辐射势；φD为绕射势。采用源偶混合分布法求解速度势的边界积分方程,表达式为

(13)

(14)

式(13)和式(14)中：G(P,Q)为满足自由面条件、水底条件和远方辐射条件的格林函数；P(x,y,z)为场点坐标；Q(ξ,η,ζ)为源点坐标。

采用移动脉动源格林函数求解边界积分方程，表达式为

(15)

(16)

(17)

(18)

2.3.2 波浪增阻计算

2.3.2.1 基于辐射能法的运动增阻计算

从辐射能法的角度考虑，在航速为U0、波速为C的情况下，船后某一位置处1个遭遇周期内波浪的辐射能量与船舶阻力平均增值所做的功相等，即

E=RaωR(U0+C)Te=RaωRλ

(19)

式(19)中：Te为遭遇周期；λ为入射波波长；E为1个波浪遭遇周期内船舶辐射能量，可表示为各模态辐射力做功之和。

(20)

将式(13)代入式(12)中，得到船舶在波浪中的运动增阻为

(21)

式(21)中：XjI和XjR分别为j模态运动幅值的实部和虚部;XkR和XkI分别为k模态运动幅值的实部和虚部；RAM为与附加质量相关的增阻；RDP为与阻尼系数相关的增阻。

当无航速时，水动力系数Ajk=Akj，即RAM=0，则船舶运动增阻为

(22)

2.3.2.2 基于Salvensen法的反射增阻计算

采用Salvensen法公式中与绕射有关的项对阻力增加进行修正，根据Newman的推导，带自由面物体所受二阶力的计算公式为

(23)

式(23)中:φB为扰动速度势;φ0为入射速度势;SB为自由面。

基于弱散射假定，仅考虑绕射势的影响，船舶在波浪中航行时的绕射增阻可表示为

(24)

同样，由式(10)可得基于三维面元法的波浪增阻计算结果。

2.4 黏流数值计算法

本文在目标船的艏部上游1.2倍船长至艉部下游约2倍船长范围内建立计算域，采用Trim网格，利用商业软件STAR-CCM+中自带的工具进行网格划分，各波长范围内有60～100个网格单元，各波高范围内有10～20个网格单元，整个计算域宽度范围内有100个网格单元，在自由液面处进行网格局部加密，计算域网格总数约为250万个(半域)。船模表面网格划分图见图3。

图3 船模表面网格划分图

采用VOF(Volume Of Fluid)方法处理自由液面，求解RANS方程。湍流模型取k方程和ε方程，入口边界给定入口流体流动速度，出口边界设置为压力出口，计算域顶部和底部边界设为速度入口条件，时间步长取0.001 s，在各时间步长内迭代10次。参照目标船波浪增阻试验工况，计算11.5 kn航速下的波浪增阻和运动响应[15]。

2.5 各方法的计算结果分析

按式(20)～式(22)对采用各方法得到的计算结果进行无因次化处理，得到波浪增阻频率响应函数KΔR(ω)、垂荡频率响应函数KZ(ω)和纵摇频率响应函数Kθ(ω)。采用各方法计算得到的频率响应结果与模型试验结果对比见图4。

(25)

(26)

(27)

式(25)～式(27)中：RW为波浪中阻力；R0为静水中阻力；Za为垂荡幅值；θa为纵摇幅值；ζa为波幅；k为波数；B为船宽；L为船长；ρ为水的密度；g为重力加速度。

a) 波浪增阻频率响应曲线

3 不规则波中的波浪增阻和运动响应预报

采用ITTC双参数谱，结合11.5 kn航速下目标船的波浪增阻和运动响应计算结果，预报不同工况下的波浪增阻和运动幅值，并将所得结果与模型试验结果相对比，结果见表2;其中:H为浪高;T为波浪周期。采用不同方法得到的目标船波浪增阻和运动响应预报结果对比见图5。

表2 各海况下采用不同方法得到的目标船波浪增阻和运动响应预报结果

a) 波浪增阻预报结果

从图5中可看出，采用各方法得到的目标船在实际海况下的波浪增阻和垂荡预报结果均较为接近，采用黏流数值计算方法得到的预报结果与模型试验结果的吻合度最高，势流三维面元法次之，采用势流二维切片法得到的纵摇预报结果与模型试验结果的偏差相对较大，这是由于二维切片法假设船体细长，忽略了船型特征，而本文选取的目标船的船型特征正是影响其垂荡和纵摇等运动响应的重要因素。

4 结 语

本文以19 000载重吨多用途船为研究对象，开展了船舶波浪增阻和运动响应标准化预报研究，分别采用势流二维切片法、势流三维面元法和黏流数值计算法计算了船舶的波浪增阻和运动响应，并结合ITTC双参数谱预报了目标船在不同海况下的波浪增阻和运动响应。研究发现，采用势流理论计算耗时较短，采用黏流数值计算法计算耗时较长且在某种程度上受计算能力的限制。通过对计算结果进行对比发现：采用各方法得到的波浪增阻和垂荡预报结果比较接近，均在工程允许的误差范围内；势流二维切片法未考虑船型特征，对本文选取的船型的纵摇预报结果与模型试验结果的偏差相对较大。

因此，在进行船型优化等前期设计工程应用时，为节约试验时间和成本，可采用势流三维面元法进行评估；在进行船舶性能验证和航线规划等工程应用时，为提升预报精准度，建议采用黏流数值计算法或通过模型试验进行预报。