船舶波浪增阻及运动响应的黏性数值计算研究

陆泽华， 董国祥， 陈伟民

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

0 引 言

船舶大多数情况下都是在波浪环境中航行的，尤其是远洋运输船舶，在航行的过程中会遭遇各种海况，海上波浪环境是影响船舶运营效能的最主要外界环境。因此，船舶在波浪中的航行性能才能更真实地反映船舶的运营效能，针对实际运营过程中遭遇的波浪环境进行船舶性能研究，将成为提高船舶能源利用效率、减少船舶能耗的重要技术途径，对提高船舶远洋航行安全也具有重要意义。

波浪中船舶阻力增加的主要研究手段包括模型试验和数值计算。对于短波中船舶的阻力增加，模型试验存在物理水池的造波机能力限制短波环境的模拟和短波具有稳定性较差等两方面困难，在水池中衰减与耗散比较严重，使得水池中较难模拟出高品质的短波，综合其他各因素的影响，导致短波中船舶阻力增加模型试验的不确定度较大。计算流体动力学(Computational Fluid Dynamics, CFD)作为一种建立在黏流理论基础上的研究方法，与势流理论相比，CFD数值计算可考虑到流体黏性的影响，且不需势流理论中的诸多假设，更符合实际的物理过程。船舶在波浪中航行时通常伴随着明显的波浪破碎及翻卷等非线性现象，势流理论无法对此类非线性现象进行研究。近年来，运用数值水池技术模拟高频波浪环境下船舶周围流场的研究正在发展。

国外，WILSON等[1]以RANS方程为控制方程，对全约束的DTMB5415船模和Wigley船模在规则波顶浪情况下的水动力进行了非定常模拟研究。SATO等[2]应用有限体积法，对Wigley和S60船模在迎浪情况下进行运动模拟分析，通过标记密度函数法求得水动力及力矩系数，进而求解船舶六自由度运动。HOCHBAUM等[3]分别研究全约束船模和自由船模在迎浪下的波场和阻力，利用有限体积法离散N-S方程和连续性方程，应用SIMPLE方法进行速度与压力的耦合。DENG等[4]利用ISIS-CFD求解器计算KVLCC2在波浪中受到的阻力，应用有限体积法离散非定常的RANS方程，近似的SIMPLE方法进行速度与压力耦合，EASM模型考虑湍流黏性的影响。

国内，齐鹏等[5]较早建立了三维数值波浪水池，通过模拟物理推板造波机的运动，成功模拟出正向入射波浪和30°浪向角的斜射波浪。周勤俊等[6]将动量方程引入入射波场，提出消除波浪二次反射的源造波-消波技术。林兆伟等[7]建立了二维数值波浪水池，模拟计算了甲板上浪问题。朱德祥等[8]模拟柔性造波板运动产生入射波,并利用人工阻尼消波，模拟了非线性规则波和不规则波。吴乘胜等[9]模拟有航速船模的辐射问题，计算得到顶浪中船模的运动和水动力性能。董志等[10]利用商业软件FLUENT及其二次开发功能，提出动边界造波和多孔介质消波等多种数值方法，建立出可有效模拟弱非线性波的数值波浪水槽。

1 数值波浪水池模型

1.1 数学模型

对数值波浪水池的建立中,将波浪流动作为两相(空气和水)流动处理，自由液面即为水和空气的交界面，使用VOF(Volume Of Fluid)方法进行处理。假设水是不可压缩流体，由于空气的速度远小于音速，可作为不可压缩流体处理。因此，数学模型的控制方程包括连续性方程、动量方程、体积分数方程、湍流模型的k方程和ε方程。

1.2 计算域建立及网格划分

建立数值波浪水池的最终目标是计算目标船波浪中的性能，计算域的设定主要考虑与波浪计算的匹配性，水池前端位于目标船艏部上游约1.2倍船长(LPP)处，尾端位于目标船艉部后下游约2LPP处，其中包含约1倍LPP的消波区，具体尺度见表1。

计算中使用的网格为Trim网格，由STAR-CCM+自带工具划分，网格划分需满足波浪数值模拟的精度，每个波长范围内网格单元数为60～100个，波高范围内网格单元数为10～20个，整个计算域宽度范围内网格单元数为100个。

表1 计算区域范围及网格划分

1.3 边界条件及初始条件

数值模拟中边界条件的具体设置为：

(1) 在入口边界上，根据目标船的前进速度和入射规则波条件，设定入口流体流动速度；

(2) 出口边界设置为压力出口，由于出口距船模足够远，波面经消波后基本平静，压力分布设置为静水压力；

(3) 计算域顶部和底部边界设置为速度入口条件，与入口边界条件相同，以更好地模拟波浪前进，减少波浪衰减；

(4) 在对称面上，满足对称条件。

以规则入射波的均匀流场作为数值计算的初始条件。

2 静水阻力数值计算

2.1 目标船型参数

与模型试验一致，缩尺比取为27.3，船模水线长为5.317 m,目标船19 000 DWT多用途船的实船与船模具体参数见表2。

表2 球艏型19 000 DWT多用途船主尺度参数 m

2.2 目标船静水阻力计算

船舶的波浪增阻是船舶在波浪中受到的总阻力与静水阻力的差值，需计算目标船的静水阻力。为简化计算模型，不考虑舵的影响。

取上述计算域范围进行19 000 DWT多用途船裸船体的静水阻力计算，在船首、船尾及开尔文波区域进行网格加密，网格总数约140万个(半域)。入口和出口边界条件设置与上述相同，计算域顶部和底部为固壁边界，对称面满足对称边界条件，流体初始速度取实船航行速度为11.5 kn，对应船模速度为1.132 m/s。时间步长取Δt=0.01 s，计算过程中监测船模受到的阻力。静水阻力时历曲线和最终流场波形分布见图1。

a) 船模(半船)静水阻力计算时历曲线 b) 船模周围流场波形分布

图1 船模静水阻力计算时历曲线及流场波形分布图

由图1中可知，船模周围流场处于随时间振荡衰减的过程，最终船模阻力稳定地处于小幅波动状态，取最终稳定阶段的阻力平均值作为船模静水阻力，计算结果为17.024 N，与试验值17.591 N间存在的误差为3.22%。其中，模型试验为带舵试验，舵的影响会对静水阻力计算结果产生一定的误差。

3 数值造波技术研究

数值波浪水池采用STAR-CCM+软件中的VOF Wave模块来模拟水池造波，该模块直接通过RANS求解器分别控制速度入口及压力出口参数，来达到模拟线性波(一阶波)的造波和消波。在线性波浪理论造波中，波浪参数满足

ω2=gktanhkh

(1)

U=aωsin(kx-ωt)ekz

(2)

式(1)和式(2)中：ω为波浪圆频率；k为波数；h为水深；a为波幅；U为速度入口边界处的流体速度。

为使波浪模拟工况尽量涵盖波浪增阻及运动响应计算中的工况，分别选取波长为2.5 m，5.0 m，8.0 m,对应波高为4 cm，5 cm，6 cm的线性波进行模拟，以规则入射波的均匀流场作为初始条件，计算时间t=10 s，取x=6 m处为监测点，监测该点波高随计算时间的变化。波长为2.5 m、波高为4 cm线性波，在初始化后和计算10 s后的自由液面波形和波高监测见图2。3种波浪工况下，在x=6 m处监测的波高时历曲线见图3。

由图3可知，利用VOF Wave模块造波具有一定的精度，波浪在传播的过程中没有出现明显衰减，满足计算目标船波浪增阻及运动响应中对波浪精度的要求。

4 目标船波浪增阻及运动响应数值计算

4.1 网格划分

参照表1建立计算域并进行网格划分，并对船舶艏部和艉部网格适当加密，舯部网格较为稀疏；船舶表面网格采用边界层数5层，其中第一层网格尺度根据y+值确定(y+值为50～100)；自由液面加密区域包括xy方向的整个计算域范围及z方向的1.2倍波高范围。根据以上原则，计算域网格总数约为250万(半域)，主要集中在自由液面加密区域；船后下游区域主要为消波区，为提高计算效率，网格较为稀疏(见图4)。

c) 艏艉部网格加密及船体表面边界层示意图

图4 目标船波浪增阻及运动响应数值计算网格划分示意图

4.2 波浪增阻及运动响应数值计算

数值计算的工况根据模型试验工况选取，速度取11.5 kn，波长船长比取0.2～1.8共计11个工况(见表3)。各工况下的波高设定参照试验波高。此外，为保证目标船波浪中的运动在各时间步长内稳定迭代，计算过程中对各时间步长的迭代次数略有增加，由静水数值计算的5次增加为10次，时间步长取0.001 s，计算过程中监测阻力、垂荡、纵摇及波高。波长船长比λ/Lpp=0.8工况下阻力、垂荡、纵摇、波高时历监测曲线及自由液面波形图见图5。

表3 模型速度及波浪参数

4.3 计算结果分析

取稳定变化段的阻力瞬时值进行平均，得到目标船在各工况下的波浪阻力值；取运动稳定变化段的垂荡、纵摇和波高幅值进行平均，得到目标船在各工况下的垂荡、纵摇和波高值，并按式(3)～式(5)无因次化后，与试验值比较(见表4和图6)。

(3)

(4)

(5)

工况增阻系数垂荡系数纵摇系数计算值试验值计算值试验值计算值试验值0.20.9530.7680.0160.0140.0040.0040.30.6510.5940.0190.0140.0070.0080.40.6420.5490.0110.0170.0080.0100.50.5670.5690.0280.0320.0120.0230.80.9480.8150.0910.0940.1420.1211.01.5331.6130.6160.4860.4910.4191.11.7201.8330.8020.7130.6320.5791.21.8341.9430.8950.8670.7690.7831.41.2981.1740.8430.7891.0060.9411.60.6560.6730.8550.7791.0580.9481.80.4190.4850.8600.8151.0681.003

由图6可知波浪增阻及运动响应黏流数值计算结果与试验值较为吻合，建立的数值计算模型有较高精度，计算结果具有参考价值。

5 结 语

基于软件STAR-CCM+，采用VOF方法捕捉自由面，通过求解RANS方程和运动方程，进行数值造波技术研究，进而研究船舶在波浪中迎浪航行时的数值模拟。分别对目标船在静水和规则波中进行计算，重点计算不同波长与船长比工况下目标船的波浪增阻及运动响应，计算结果与试验值较为吻合，表明全黏流数值计算模型能准确预报出增阻和运动响应的变化趋势。研究工作为波浪中航行船舶的波浪增阻和运动响应研究提供了一种有效的技术手段，可为波浪中船舶快速性预报和船型优化提供技术支撑。

[1] WILSON R, PATERSONE, STERN F. Unsteady RANS CFD for Naval Combatants in Waves[C]∥Proceeding of the 22nd Symposium on Naval Hydrodynamics, 1998： 532-549.

[2] SATO Y， MIYATA H, SATO T. CFD Simulation of 3-Dimentional Motion of a Ship in Waves: Application to an Advancing Ship in Regular Heading Waves[J]. Journal of Marine Science and Technology， 1999，4(3)： 108-116.

[3] HOCHBAUM A C， VOGT M. Towards the Simulation of Seakeeping and Manoeuvring Based on the Computation of the Free Surface Viscous Ship Flow[C]∥Proceeding of the 24th Symposium on Naval Hydrodynamics, 2002：268-281.

[4] DENG G B, QUEUTEY Y， VISONNEAU M. RANS Prediction of the KVLCC2 Tanker in Head Waves[J]. Journal of Hydrodynamics， 2010，22(5)： 476-481.

[5] 齐鹏，王永学. 三维数值波浪水池技术与应用[J]. 大连理工学报，2003,43(6)：825-830.

[6] 周勤俊，王本龙，兰亚梅， 等. 海堤越浪的数值模拟[J]. 力学季刊， 2005,26(4)：629-633.

[7] 林兆伟，朱仁传，缪国平. 甲板上浪问题的二维数值模拟[J]. 船舶力学，2009,1(1)：1-8.

[8] 朱德祥，沈泓萃，洪方文， 等. 船模数值水池框架及其研究基础[J]. 水动力研究与进展，2008,1(1)：24-32.

[9] 吴乘胜，朱德祥，顾民. 数值波浪水池中船舶顶浪运动模拟研究[J]. 船舶力学，2008,5(2)：692-696.

[10] 董志，詹杰民. 基于VOF方法的数值波浪水槽以及造波、消波方法研究[J]. 水动力学研究与进展，2009,24(1)：15-21.

[11] 陆泽华，李建鹏，李传庆，等. 波浪自航与非自航试验中的船舶运动对比研究[C]//中国造船工程学会.船舶力学学术委员会第八次全体会议论文集.北京：船舶力学学校委员会，2014:157-162.