船舶在波浪上的运动与阻力增加一直是船舶与海洋工程领域研究的热点，准确预报船舶在波浪上的运动是船舶在水上航行、作业安全性和舒适性的基础，船舶波浪增阻在船型优化乃至航线优化问题中是经济性的重要参数指标。

目前，关于波浪中航行船舶运动与增阻研究主要有船模试验(experimental fluid mechanics, EFD)、基于势流理论的理论计算^[1]和考虑黏性的计算流体力学(CFD)数值模拟等方法。基于势流理论的方法主要有切片法、三维面元法等，相对简单且计算效率高，但势流理论没有考虑粘性的影响，且对诸如船体大幅度运动、波浪破碎等强非线性因素难以适用。而船模试验往往需要进行费时费力的系列模型试验，成本较高，试验频率也有限制，由于仪器设备和测试技术等方面的原因，目前还难以准确地监测出船体周围复杂流场的信息。

CFD方法由于考虑了粘性能对非线性问题作较为准确的数值模拟，随着计算机计算能力的提高日渐成为求解船舶与海洋工程水动力问题的一种重要手段，近年来已经取得了一系列重要的进展。Carrica等^[2-3]利用重叠网格采用单相 level-set方法捕捉自由面，考虑了纵摇和垂荡两个自由度，对船模DTMB 5512在波浪中的自由运动的进行了数值模拟。GUO等^[4-5]研究了KVLCC2在短波中的增阻，验证了Faltinsen渐进公式预报船舶在短波中增阻的准确性，同时还对KVLCC2的耐波性做了研究，并分析了波浪对螺旋桨盘面处流场的影响。Sadat-Hosseini等^[6]考虑纵荡的影响，数值预报超大型油轮KVLCC2在不同波长迎浪工况下的波浪增阻和运动响应，结果表明纵荡对纵摇和波浪增阻几乎没有影响。国内基于CFD方法关于船舶在波浪上运动与增阻的研究也取得了很多成果。方昭昭等^[7-9]基于CFD方法建立了数值波浪水池，使用Fluent软件结合动网格技术对不同航速顶浪航行的Wigley-III型船模的水动力系数、垂荡、纵摇、波浪增阻等进行了数值计算，结果与实验结果^[10]吻合较好。沈志荣等^[11-12]利用基于开源代码OpenFOAM平台开发计算分析了Wigley-III 型、DTMB5415、S175等船型在迎浪中的运动响应以及波浪增阻。

CFD模拟船体运动需要使用动网格技术，涉及到网格的再生与变形，若船体运动响应较大，网格畸变容易引发数值计算发散或是计算结果误差大。重叠网格方法对于复杂曲面离散以及物体大幅运动具有无可比拟的优势，能方便处理复杂三维船体曲面，容易保证局部网格质量，动态重叠网格在处理结构物具有大幅运动的绕流问题时不需要网格再生，提高了动态网格处理效率。

本文基于粘性流理论使用CFD软件建立了数值波浪水池并结合重叠网格方法，对航速为1.047 m/s(Fn=0.142 )缩尺比为58的KVLCC2模型在波浪上的运动进行了数值模拟，同时研究了KVLCC2在波浪中的阻力增加，分析了波阻成分。

1 数值波浪水池与重叠网格技术

本文的数值模拟是在一个三维数值波浪水池中进行，水池尺度随入射波长的不同而不同，在出口处和两侧边界附近设有人工阻尼消波区，消波区段约为2倍波长，船侧消波区距船舷侧约1倍船长，水池前端距船艏约1倍波长，尾端消波区距船艉约1倍波长。数值波浪水池及坐标系如图 1所示，固定坐标系O₀x₀y₀z₀原点位于船体中纵剖面、中横剖面以及静水面交点，x₀轴指向船艏，y₀轴指向船舷左侧，z₀轴垂直向上；参考坐标系Oxyz原点位于重心处，各坐标轴方向与固定坐标系对应一致。

1.1 流体控制方程

整个流场以连续性方程和N-S方程为控制方程：

$ \frac{{\partial \rho }}{{\partial t}} + \frac{{\partial \left( {\rho {{\mathit{\boldsymbol{\bar u}}}_i}} \right)}}{{\partial {x_i}}} = 0\;\;\;\left( {i = 1,2,3} \right) $

(1)

$ \begin{array}{*{20}{c}} {\frac{{\partial \left( {\rho {{\mathit{\boldsymbol{\bar u}}}_i}} \right)}}{{\partial t}} + \frac{{\partial \left( {\rho {{\mathit{\boldsymbol{\bar u}}}_i}{{\mathit{\boldsymbol{\bar u}}}_j}} \right)}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left[ {\mu \frac{{\partial {{\mathit{\boldsymbol{\bar u}}}_i}}}{{\partial {x_j}}} - \rho \overline {{{\mathit{\boldsymbol{u'}}}_i}{{\mathit{\boldsymbol{u'}}}_j}} } \right] - }\\ {\frac{{\partial \bar p}}{{\partial {x_i}}} + \rho {f_i}\;\;\left( {i,j = 1,2,3} \right)} \end{array} $

(2)

式中：u_i为流体微团在i方向上的速度，f_i为质量力，为流体压力，流体密度定义为$\rho = \sum\limits_{q = 1}^2 {{a_q}{\rho _q}} $，其中体积分数a_q表示单元内第q相流体体积占总体积的比例，并且有$\sum\limits_{q = 1}^2 {{a_q} = 1} $，μ为相体积分数的平均动力粘性系数，与密度定义的形式一致。

湍流模式为SST k-ω两方程模型，其控制方程为

$ \begin{array}{*{20}{c}} {\left\{ \begin{array}{l} \frac{\partial }{{\partial t}}\left( {\rho k} \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho k{{\mathit{\boldsymbol{\bar u}}}_i}} \right) = \frac{\partial }{{\partial {x_j}}}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_k}\frac{{\partial k}}{{\partial {x_j}}}} \right) + {\mathit{\boldsymbol{G}}_k} - {\mathit{\boldsymbol{Y}}_k}\\ \frac{\partial }{{\partial t}}\left( {\rho \omega } \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho \omega {{\mathit{\boldsymbol{\bar u}}}_i}} \right) = \frac{\partial }{{\partial {x_j}}}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_\omega }\frac{{\partial \omega }}{{\partial {x_j}}}} \right) + {\mathit{\boldsymbol{G}}_\omega } - \\ \;\;\;\;{\mathit{\boldsymbol{Y}}_\omega } + {\mathit{\boldsymbol{S}}_\omega } \end{array} \right.}\\ {\left( {i = 1,2,3} \right)} \end{array} $

(3)

式中：Γ_k、Γ_ω为湍动能k和比耗散率ω的有效扩散系数，Y_k、Y_ω为k和ω湍流耗散，G_k为平均速度梯度引起k的产生项，G_ω为ω的产生项，S_ω为交叉扩散项。

流体体积输运方程为

$ \frac{{\partial {a_q}}}{{\partial t}} + \frac{{\partial \left( {{{\mathit{\boldsymbol{\bar u}}}_i}{a_q}} \right)}}{{\partial {\mathit{\boldsymbol{x}}_i}}} = 0\;\;\;\;\;\left( {q = 1,2;i = 1,2,3} \right) $

(4)

本文在后面的数值计算中采用有限体积法离散上述控制方程，对流项使用二阶迎风插值格式，扩散项的离散采用中心差分格式，应用SIMPLE算法分离式求解；考虑重力影响，使用多重网格方法迭代求解离散代数方程组。

1.2 数值造波与消波

船舶在波浪中航行时的运动模拟，需要提供满足计算要求的波浪环境，即首先要构建数值波浪水池，本文采用速度边界法造波。假定波浪沿x轴负方向传播，根据波浪理论，参考坐标系下入射势、遭遇频率以及波面表达式为

$ \left\{ \begin{array}{l} \mathit{\Phi } = {C_0}{\xi _a}\frac{{\cosh \left[ {k\left( {z + H} \right)} \right]}}{{\sinh \left( {kH} \right)}}\left( {kx + {w_e}t} \right)\\ {w_e} = {w_0} - k{U_0}\\ \eta = {\xi _a}\cos \left( {kx + {w_e}t} \right) \end{array} \right. $

(5)

参考坐标系下流体质点的速度表达式为

$ \left\{ \begin{array}{l} u = {w_0}{\xi _a}{{\rm{e}}^{kz}}\cos \left( {kx + {w_e}t} \right)\\ v = 0\\ w = {w_0}{\xi _a}{{\rm{e}}^{kz}}\sin \left( {kx + {w_e}t} \right) \end{array} \right. $

(6)

在数值波浪水池的尾部设置有人工阻尼消波区，其长度约为入射波长的2倍。在消波区，对流体质点垂向速度做强迫衰减，为保证流动的连续性，水平速度不做衰减处理^[13]。强迫衰减的公式为

$ \left\{ \begin{array}{l} {w_r}\left( {x,y,z;t} \right) = w\left( {x,y,z;t} \right)\mu \left( {x,z} \right)\;\;\;\;\;\;\;\;\;\left( 7 \right)\\ \mu \left( {x,z} \right) = \alpha {\left( {\frac{{x - {x_s}}}{{{x_e} - {x_s}}}} \right)^2} \cdot \frac{{{z_b} - z}}{{{z_b} - {z_f}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( 8 \right) \end{array} \right. $

式中：w_r(x, y, z;t)为衰减后的垂向速度；μ(x, z)为衰减函数；x_s≤x≤x_e，z_b≤z≤z_f；下标s和e分别代表阻尼区沿x方向的起点和终点；b和f分别代表沿z方向的底部和自由面；α为阻尼控制参数。

当流场中有船体存在时，还要在数值波浪水池的侧面设置人工阻尼消波区，同样只对流体质点垂向速度做强迫衰减，强迫衰减公式为

$ \left\{ \begin{array}{l} {w_r}\left( {x,y,z;t} \right) = w\left( {x,y,z;t} \right)\eta \left( {y,z} \right)\;\;\;\;\;\;\;\;\left( 9 \right)\\ \eta \left( {y,z} \right) = \alpha {\left( {\frac{{y - {y_w}}}{{{y_n} - {y_w}}}} \right)^2} \cdot \frac{{{z_b} - z}}{{{z_b} - {z_f}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {10} \right) \end{array} \right. $

式中：η(y, z)为衰减函数；y_w≤y≤y_n，z_b≤z≤z_f；下标w和n分别代表阻尼区沿y方向的起点和终点；b和f分别代表沿z方向的底部和自由面；α为阻尼控制参数。

1.3 重叠网格技术

为更好地模拟船体在波浪上的运动，本文采用重叠网格技术。重叠网格，也叫嵌合体网格，是一种区域分割与网格组合的策略，涉及到两种区域：背景区域和重叠区域，两种区域以任意方式叠加。如图 2所示，整个流场称之为背景区域，该区域生成的网格为背景网格；包裹着船体随船体做刚性运动的小区域为重叠区域，这个区域产生的网格称为重叠网格，两套网格之间通过交界面来识别。

重叠网格技术将网格分为有效单元和无效单元。只有有效单元参与离散控制方程的求解，无效单元指的是重叠区域内被挖去的背景网格，不参与控制方程的求解，当重叠区域随船体做刚性运动时背景网格中的部分有效单元和无效单元会相互转化。重叠网格与背景网格的信息交换通过受者单元和贡献单元之间插值来完成的，一般采用拉格朗日插值的思想，进行线性插值，具有以下形式^[14]：

$ {\varphi _{{\rm{rec}}}} = \sum {{\xi _i}{\varphi _i}} $

(11)

式中：φ_rec为受者单元的流动物理量；φ_i为其相邻贡献单元的流动物理量；对于二维问题，ξ_i为以相邻贡献单元中心为顶点组成三角形对应的形函数；如果是三维问题，ξ_i是四面体所对应的形函数，该四面体的顶点由相邻贡献单元的中心组成。

1.4 船体运动方程

本文数值模拟了KVLCC2在长波中迎浪航行工况下的运动，纵荡运动较小且对其他运动以及波浪增阻几乎没有影响，此处只考虑垂荡和纵摇两个自由度，其运动方程为

$ \left\{ \begin{array}{l} m\frac{{{\rm{d}}\omega }}{{{\rm{d}}t}} = {F_3}\\ M\frac{{{\rm{d}}\dot \beta }}{{{\rm{d}}t}} = {F_5} \end{array} \right. $

(12)

式中:m为船体的质量，ω为船体重心处的垂向速度，F₃为船体受到的合力的垂向分量, M为船体绕旋转轴(过重心平行于y轴)的转动惯量，${\dot \beta }$为纵摇角速度，F₅为作用于船体绕旋转轴的合力矩。

为避免作用于船体的力或力矩过大造成船体剧烈的运动，在初始时给船体受到的力和力矩乘上一个阻尼函数f_r，定义如下

$ {f_r} = \left\{ \begin{array}{l} \;\;0,\;\;\;\;\;\;\;\;\;t < {t_s}\\ \frac{{t - {t_s}}}{{{t_r}}},\;\;\;\;\;{t_s} \le t < {t_s} + {t_r}\\ \;\;\;1,\;\;\;\;\;\;\;\;t \ge {t_s} + {t_r} \end{array} \right. $

(13)

式中：t_s为船体可以开始运动的时间，t_r为船体开始运动后作用于船体的力和力矩受到人工抑制的时间。

2 波浪中航行船舶的数值模拟

本文数值模拟的是迎浪工况，流场关于船体中纵剖面对称，计算域取流场的一半即可。数值模拟计算主要考虑船体垂荡和纵摇两个自由度，部分短波工况，也采用固定模进行了数值模拟。

2.1 船体几何参数

KVLCC2是II型MOERI(maritime and ocean engineering research institute)油轮，相较于广泛研究的Wigley、S60以及S175船模具有更复杂的船体外形，其主要船型参数如表 1所示，各站横剖面如图 3所示。

2.2 边界条件与网格划分

波浪中航行船舶水动力计算域的边界条件设定如下：前端及上下面均为速度入口；左右为对称边界条件；长波时尾端采用压力出口，短波时尾端可以视为远场，用速度入口作为边界条件；船体设定为壁面边界条件。

计算区域网格的划分主要从如下几个方面考虑：沿波浪传播方向保证足够数量的网格，以避免数值耗散引起的波浪幅值的沿程衰减；垂向沿自由面附近要保证足够数量的网格，以精准捕捉自由液面；为准确捕捉尾迹，在开尔文兴波区范围内也进行了网格加密；在流场变化比较剧烈的地方如球鼻艏以及船艉附近要进行网格加密；为了获得准确的船体周围流场信息，船体附近网格划分较远场网格划分应精细一些；在消波区，保证网格平缓过渡的情况下适当拉大x和y方向网格的尺度，既减少计算量又可以起到数值消波的作用。船体壁面设置5层边界层网格，y⁺取50左右，为减少网格量，甲板不设边界层。

主要工况计算所用网格参数及时间步如表 2所示，计算所用网格划分如图 4所示。

2.3 数值模拟结果的处理

本文采用傅里叶级数展开法，使用傅里叶级数对数值模拟得到的垂荡、纵摇以及阻力时历曲线进行拟合。

以工况为Fn=0.142、L=4.965 m、ξ_a=0.075 m时的波浪力时历曲线为例，如图 5所示(运动时历曲线处理分析方法与此一致)，待流场达到稳定后，选取10个周期的波浪力时历曲线用傅里叶级数进行拟合，波浪力关于时间的函数可以表达为

$ {R_{{\rm{wave}}}} = {R_t} + {R_{1{\rm{st}}}}\sin \left( {{\omega _e}t + {\gamma _1}} \right) + \cdots $

(14)

式中：R_wave为波浪力；ω_e为遭遇频率；R_t为拟合得到的傅里叶级数的定常值，即平均波浪力；R_1st为一阶波浪力幅值；γ₁为一阶波浪力的相位。

规则工况，垂荡、纵摇和阻力随时间变化具有周期性，且其周期与波浪遭遇周期相等。对于阻力，关注的是平均阻力，即拟合阻力时历曲线得到的傅里叶级数表达式中常数项就是船舶在波浪中航行的平均阻力R_t；对于垂荡和纵摇，关注的是运动的幅值和相位，即拟合运动时历曲线得到的傅里叶级数表达式中各阶系数和相应相位。

3 KVLCC2在波浪中的运动与增阻

势流理论使用的是由上海交通大学水动力学实验室基于二维势流理论开发的切片程序计算的结果，实验结果取自大阪大学(OU)、意大利船舶研究所(INSEAN)以及挪威科技大学(NTNU)三家单位发表在相关文献上的实验数据^[15]。

3.1 垂荡和纵摇

采用傅里叶级数拟合运动时历曲线分析可得不同频率遭遇波下船体运动响应的幅值和相位。KVLCC2垂荡和纵摇幅值分别按式(15)进行无因次化，其幅频响应曲线如图 6、7所示。

$ {H_{{\rm{RAO}}}} = \frac{{{X_3}}}{{{\zeta _a}}},{P_{{\rm{RAO}}}}\frac{{{X_5}}}{{K{\zeta _a}}} $

(15)

式中：ξ_a为遭遇波波幅，K为波数，X₃和X₅分别为垂荡和纵摇的幅值。

从图 6、7中可以看出，数值模拟的结果与实验结果整体吻合较好。在0.7＜L_pp/λ＜1.1时，随着波长减小，垂荡幅值和纵摇幅值都迅速减小，在L_pp/λ=1.4时垂荡和纵摇运动已经可以忽略不计了，说明波长较短时，由波浪激励引起的船体辐射运动很小。对于纵摇运动，在峰值附近，CFD计算值更接近INSEAN水池的试验值，切片法的计算值更接近大阪大学的试验值，两个水池的试验结果略有差异。

图 8、9给出了KVLCC2垂荡和纵摇运动与遭遇波的相对相位随L_pp/λ变化情况，数值模拟的结果同实验结果整体符合较好。1.0＜L_pp/λ＜1.1时，垂荡相位存在一个波谷；0.75＜L_pp/λ＜1时，垂荡相位随波长变化比较明显；L_pp/λ＜0.75时，随L_pp/λ减小垂荡相位缓慢增加，最终收敛于0°附近。纵摇相位在0.5＜L_pp/λ＜1和L_pp/λ＞1两个区间内分别随L_pp/λ增大而单调递减；0.5＜L_pp/λ＜1时随着L_pp/λ的减小垂荡相位有收敛于-90°的趋势。对于长波工况，当船中附近波高为0时，波面高度沿船长分布单调(递增或递减)，纵倾角达到最大或最小值，此时船艏与船艉受到的浮力随船舶纵倾发生相反的变化正好互相补偿，所以垂荡运动处于平衡位置；当船中附近为波峰或者波谷时，艏艉浮力相当，船舶表现为正浮，总的浮力达到最大值或者最小时，船体需要通过上浮或者下沉来平衡重力；数值模拟和实验的结果都反映了长波工况时，船体处于随浪运动状态，跟真实的物理情况相符。

3.2 KVLCC2在波浪中的增阻

根据定义，船舶在波浪中的增阻R_aw等于船舶在波浪中航行的平均阻力R_t减去同一航速下在静水中航行的阻力R_calm，即R_aw=R_t-R_calm，R_t通过使用傅里叶级数拟合阻力时历曲线获得。为了便于不同数值模拟方法计算的结果以及实验的结果相互之间的比较，波浪增阻值按式(19)进行了无因次化：

$ {C_{aw}} = \frac{{{R_{aw}}}}{{\rho g\zeta _a^2{B^2}/{L_{{\rm{PP}}}}}} $

(16)

图 10为Fn=0.142 时KVLCC2波浪增阻的CFD模拟结果与切片法计算结果以及实验结果的对比，其中切片法采用辐射能量法计算波浪增阻，可以看到，CFD模拟的结果同实验结果吻合较好。切片法计算的增阻值整体上小于实验值，增阻幅值同实验值差异较大，主要原因是：势流方法没有考虑粘性的影响，流体的粘性也贡献了一部分增阻；本文势流切片方法采用辐射能量法计算波浪增阻，没有考虑绕射对增阻的贡献。

从作用力方向分析，总的波浪增阻由压力增阻、摩擦增阻两部分组成；从物理问题角度出发，波浪增阻又可分解为辐射增阻和绕射增阻。图 11给出了遭遇频率对增阻成分的影响。由图 11(a)可以看到，在全波长范围内波浪增阻基本上由压力增阻贡献，由摩擦力引起的增阻很小。从图 11(b)可以看出，波长对绕射增阻影响较小，辐射增阻在L_pp/λ＞1.3时趋于0，在长波范围，辐射增阻随着波长的增大而减小。

3.3 KVLCC2在短波中的增阻

由上述分析可知，船舶在短波中的增阻主要由绕射增阻贡献，L_pp/λ＞1.3时辐射增阻可以忽略不计，现在对KVLCC2在短波中的增阻成分作进一步验证。

图 12中，CFD方法计算的是Fn=0.142 时KVLCC2固定模在短波中的增阻值，实验测量的是Fn=0.142 时KVLCC2自由模在短波中的增阻值。其中“Faltinsen”是用Faltinsen的渐近公式^[14]计算得到的波浪增阻，顶浪工况下渐进公式为

$ {R_{aw}} = \int\limits_{C\_l} {\frac{1}{2}\rho g\zeta _a^2\left[ {{{\sin }^2}\theta + 2\frac{{{\omega _0}U}}{g}\left( {1 - {{\cos }^2}\theta } \right)} \right]\sin \theta {\rm{d}}l} $

(17)

式中:C_l表示非遮蔽水线，ω₀是波浪的频率，θ是水线的切线与x轴的夹角，ζ_a表示波幅。

可以看到，KVLCC2在短波范围内的增阻值随波长基本变化不大，不计入辐射增阻的CFD方法计算值以及Faltinsen渐进公式的计算值都同实验值吻合的非常好，也说明了本文数值模拟方法预报船舶在短波中的波浪增阻相当有效。

4 结论

1)结合重叠网格技术的CFD方法能够准确处理船舶在波浪上的运动响应，相较于势流理论，CFD方法能够在全波长范围内更准确地预报船舶在波浪上增阻；

2)船体遭遇短波时，船体辐射运动非常小，可以忽略，随着波长增加，垂荡幅值和纵摇幅值都迅速增大，船体在长波中航行时无因次的运动幅值趋于1；遭遇足够长的入射波，船体垂荡与波浪同步，相位差0°，纵摇运动相位则落后波浪四分之一周期，船体处于“随浪逐流”状态；

3)L_pp/λ=0.9左右时波浪增阻达到峰值；在短波中船体在波浪中的运动响应比较小，波浪增阻主要由绕射增阻贡献，辐射增阻可以忽略不计；长波范围，船体辐射运动引起的增阻随波浪增大逐渐减小；船舶在波浪中的增阻主要由压力场的改变引起，摩擦对增阻的贡献非常小。