0 引　言

水翼艇是近年来快速发展的一类高性能船舶，具有阻力低、经济航速高、耐波性好等优点^[1]，在军用与民用方面都有较为广阔的市场前景和发展空间。

水翼艇在航行时水翼产生升力将艇体托出水面, 有效地减小了阻力。Suastika等^[2]研究发现在高航速下，与滑行艇相比，安装水翼减阻率可达10%。随着水翼艇的发展，国内外学者针对水翼艇的水动力性能开展了大量研究。周俊麟等^[3]通过模型试验研究了水翼翼型、安装高度与攻角等参数对水翼艇阻力性能的影响，为水翼艇研究设计提供参考。Li等^[4]基于CFD方法对水翼双体船在迎浪航行中的阻力特性与运动响应进行研究，并进行实船试验，研究结果表明，合理设置水翼布局方案的水翼方案可使船舶总阻力降低26%，同海况条件下航速提升达21%。Bi等^[5]研究了首部水翼对水翼艇耐波性能的影响，并与无水翼的裸艇进行对比，结果表明水翼可以有效减小波浪中运动响应的幅度。Haekal等^[6]讨论了尾部水翼攻角对水翼艇阻力性能的影响，研究发现，Fr在0.5～0.75范围内，加装合适的尾翼可以减阻9%～26%。Liao等^[7]对水翼船水翼支柱的几何形状进行优化，优化后的水翼有效抑制了空化的产生。

海上航行的船舶通常会遭遇波浪，波浪的扰动会使船舶受到的阻力较静水中有所增加，这种现象被称为波浪增阻^[8]，波浪增阻的研究方法主要有试验方法、势流方法以及计算流体力学方法。Naito等^[9]对波浪中船首水翼与波浪增阻之间的关系进行试验研究，试验结果表明安装水翼船舶在波浪中的性能优于没有安装水翼船舶。基于势流理论的方法主要有切片法^[10]与三维面元法^[11]。Havelock^[12]最早使用纵向积分的方法对波浪中的无航速船形水平切片圆柱体的阻力增值进行研究。随着计算机技术的发展，计算流体力学方法在预报船舶波浪增阻方面得到了广泛应用。Steen等^[13]基于CFD方法对比加装水翼与裸艇体在波浪中的性能，并与试验结果进行对比，结果表明在各个波浪周期下，数值模拟的阻力值低于试验值，主要是由于试验过程中水翼安装方式导致的误差。Cai等^[14]基于CFD方法对高速双体船的波浪增阻及纵向运动进行数值模拟，并与势流方法进行对比，结果表明CFD方法的计算精度更高，而势流方法在波浪增阻的预报上略有误差。

目前，国内外对水翼艇的试验与数值研究大多针对单水翼船（翼滑艇）的阻力性能，对双水翼艇在波浪中的运动与波浪增阻研究较少。为此本文以某水翼艇为研究对象，基于CFD方法研究该艇在迎浪规则波中不同波浪工况下的波浪增阻与运动响应的变化，分析波浪因素对水翼艇波浪增阻的影响，最后研究水翼安装位置对水翼艇波浪增阻的影响。

1 数值方法 1.1 控制方程

雷诺平均纳维斯托克斯方程（RANS方程）是目前较为广泛使用的湍流数值模拟方法。本文采用雷诺平均方法对水翼艇的水动力性能进行求解。在笛卡尔坐标系下，RANS方程的表达形式为：

$ \frac{{\partial {u_i}}}{{\partial {x_i}}} = 0,\quad i = 1,2,3，$

(1)

$ \rho \frac{{\partial {u_i}}}{{\partial {x_i}}} + \rho {u_j}\frac{{\partial {u_i}}}{{\partial {x_j}}} = - \frac{{\partial p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left( {\mu \frac{{\partial {u_i}}}{{\partial {x_j}}} - \rho \overline {u_i^\prime u_j^\prime } } \right) + {S_i}。$

(2)

式中：$ \rho $为密度；$ p $为压力；$ {u_i} $和$ {u_j} $为速度分量时均值；$ {x_i} $和$ {x_j} $均为空间坐标分量；$ \mu $为动力粘度系数；$ t $为时间；$ {u_i} $与$ {u_j} $为速度分量时均值；$ - \rho \overline {u_i^\prime u_j^\prime } $为雷诺应力项；$ {S_i} $为动量方程广义源项。

1.2 湍流模型

在式(2)中出现雷诺应力项，导致方程组不封闭，需对雷诺应力项进行求解或模拟。本文采用Realizable k-ε模型^[15]对雷诺应力进行模拟。该模型与其他模型相比可节省计算时间，并有研究表明Realizable k-ε模型在处理航速较高雷诺数较大的船舶水动力问题时，其结果与试验值更为接近^[16]。

Realizable k-ε模型控制方程为：

$ \frac{{\partial (\rho k)}}{{\partial t}} + \frac{{\partial (\rho k{u_i})}}{{\partial {x_j}}} = \frac{\partial }{{\partial {x_j}}}\left[(\mu + \frac{{{\mu _t}}}{{{\sigma _k}}})\frac{{\partial k}}{{\partial {x_j}}}\right] + {G_k} - \rho \varepsilon，$

(3)

$ \begin{split} \frac{{\partial (\rho \varepsilon )}}{{\partial t}} + \frac{{\partial (\rho k{u_i})}}{{\partial {x_j}}} =& \frac{\partial }{{\partial {x_j}}}\left[(\mu + \frac{{{\mu _t}}}{{{\sigma _\varepsilon }}})\frac{{\partial \varepsilon }}{{\partial {x_j}}}\right] +\\ &\rho {C_1}E\varepsilon - \rho {C_2}\frac{{{\varepsilon ^2}}}{{k + \sqrt {v\varepsilon } }} 。\end{split}$

(4)

式中：系数$ {C_1} = \max [0.43,\eta /(\eta + 5)] $；$ \eta = Sk/\varepsilon $；S为平均应变率张量的模数；$ {C_2} = 1.9 $；$ {\sigma _k} = 1.0 $；$ {\sigma _\varepsilon } = 1.2 $；$ {G_k} $为k的产生项。

1.3 自由液面处理方法

在数值计算中，使用了VOF方法和高分辨率界面捕捉来捕捉自由液面：

$ \frac{{\partial {\alpha _1}}}{{\partial t}} + \nabla \left( {{\alpha _1}\vec U} \right) = 0。$

(5)

式中：$ {\alpha _1} $为气体体积分数；$ \vec U $为流场速度矢量场。

1.4 数值造波与消波

常见的造波方法分为两类：一是源造波方法，包括质量以及动量源造波；二是边界模拟法，主要包括模拟物理模型的推板造波、摇板造波、明渠造波、以及边界流体流动速度函数的直接输入法^[17]。本文采用速度边界法造波。设定波浪沿x轴负方向传播，根据波浪理论，规则波的波面方程与速度场分别表达式为：

$ \eta = {\xi _a}\cos (kx - {\omega _0}t), $

(6)

$ \left\{ \begin{gathered} u = {\omega _0}{\xi _a}{e^{kz}}\cos (kx - {\omega _0}t)，\\ v = 0 ，\\ w = {\omega _0}{\xi _a}{e^{kz}}\sin (kx - {\omega _0}t)。\\ \end{gathered} \right. $

(7)

式中：$ {\xi _a} $为波幅；k为波数；$ {\omega _0} $为波浪圆频率；x轴方向设置为波浪传播方向；通过在入口边界将式(5)作为边界条件输入来实现造波。

为了避免波浪到达计算域出口时产生反射波，对数值模拟产生干扰，因此在压力出口处设置了消波区域，用于减小出口所产生的反射波，消波区的长度取波长的1～2倍。

2 数值模型及网格划分 2.1 几何模型

本文以一艘高速水翼艇为研究对象，计算模型的主尺度及几何模型如表1与图1所示，仿真对象采用缩尺比为2的模型。前后水翼剖面均采用NACA6409翼型，首翼为V型水翼，上反角为45°。尾翼为梯型水翼，裸船体的三维模型如图1(a)所示，V型水翼的三维模型如图1(b）所示，梯型水翼的三维模型如图1(c)所示，水翼攻角为4.7°。

2.2 计算域及网格划分

采用 DFBI方法模拟水翼艇的纵摇和垂荡运动，背景域与重叠域均采用长方体。为了降低计算成本，提升计算效率，采用对称边界，对半船体进行计算。计算域的边界条件设置如图2所示，其中入口边界距水翼艇艇体为2.5L，出口边界距艇体为4L，侧部边界距艇体为2L，艇体距离上边界为1L，距离下边界为2.5L；中纵剖面所在平面设置为对称平面。入口边界、顶部以及底部边界条件设置为速度入口，出口边界设置为压力出口。

网格离散采用切割体与棱柱层网格相结合的形式，并采用重叠网格技术，在艇体附近设置重叠网格区域。为更好地捕捉自由表面处与流场，在相关位置采用六面体密度盒进行局部加密。在距离艇体较远的区域，网格尺寸逐渐放大，以此减小网格数量，达到节省计算时间的效果。计算域整体以及艇体表面的网格划分网格情况如图3所示。

2.3 数据处理

水翼艇在规则波中航行时，其纵摇与垂荡具有周期变化的特点，本文采用傅里叶级数展开法对周期性运动进行处理。

根据傅里叶级数展开原理，随时间变化的函数按傅里叶系数展开得到级数，其形式为：

$ {\varphi _t} = \frac{{{\varphi _0}}}{2} + \sum\limits_{n = 1}^N {{\varphi _n}} \cos (n{\omega _e}t + {\gamma _n})，n = 1,2\cdots。$

(8)

式中：$ {\varphi _t} $为分析对象；$ {\varphi _0} $、$ {\varphi _1} $…$ {\varphi _n} $为各阶傅里叶展开级数；$ {\gamma _n} $为相位。其相关计算式为：

$ a_n=\frac{2}{{T}}\int_0^{{T}}\varphi(t)\cos(\omega_ent)\mathrm{d}t，$

(9)

$ b_n=-\frac{2}{{T}}\int_0^{{T}}\varphi(t)\sin(\omega_ent)\mathrm{d}t，$

(10)

$ {\varphi _n} = \sqrt {a_n^2 + b_n^2}，$

(11)

$ {\gamma _n} = {\tan ^{ - 1}}\left(\frac{{{b_n}}}{{{a_n}}}\right)。$

(12)

根据数值模拟得到的垂荡和纵摇幅值计算垂荡和纵摇运动响应：

$ C_{\mathrm{heave}RAO}=\dfrac{x_{\mathrm{heave}}}{\xi_a}，$

(13)

$ C_{\mathrm{pitch}RAO}=\dfrac{x_{\mathrm{pitch}}}{\xi_ak}。$

(14)

式中：$ C_{\mathrm{heave}RAO} $、$ C_{\mathrm{pitch}RAO} $分别为垂荡和纵摇的传递函数；$ x\mathrm{_{heave}} $、$ x\mathrm{_{pitch}} $指的是垂荡、纵摇和横摇运动的傅里叶展开式中的一阶幅值；$ {\xi _a} $为波浪的幅值；$ {\xi _a}k $为波陡；$ k $为波数，$ k = 2 {\text{π}} /\lambda $。

水翼艇在波浪中的增阻$ {R_{aw}} $为水翼艇在迎浪规则波中的平均阻力减去水翼艇在静水中的阻力$ {R_{calm}} $，得到波浪增阻系数$ {C_{aw}} $的无因次化公式:

$ {C_{aw}} = \frac{{{R_{aw}} - {R_{calm}}}}{{\rho g{\xi _a}^2{B^2}/L}}。$

(15)

2.4 数值水池波浪验证

为验证数值造波的精度，对规则波进行数值模拟，并将计算结果与理论值进行对比。以波高H=0.06 m、$\lambda /L = 3$处验证对象，在距离入口边界2.875、7.875 m处设置监测点以监测波高的变化，分别为测点1与测点2。

图4为自由面波高云图，图5分别给出了2个波高监测点的时历曲线，图中虚线为对应的理论解，实线为数值模拟的结果，不同监测点监测到的波高与理论波高对比的相对误差均小于 6%。

2.5 网格收敛性分析

为验证数值模拟方法的有效性，对数值结果进行网格收敛性分析。采用波长船长比为 1，波高为0.02 m，$F{r_\nabla }$=1.659的计算工况。

通过改变自由液面附近水平方向上和垂直方向上的网格尺度，分为高、中、低密度三重系列网格模型，其中中等密度网格的基础尺寸设置为0.079 m。3套网格的网格数目分别为6.70×10⁶、2.91×10⁶和 1.27×10⁶，3套网格水平方向及垂直方向的尺度参数见表2，${G_1}$、${G_2}$、${G_3}$分别对应高密度、中密度和低密度网格划分方案，不同网格下的波浪增阻、垂荡和纵摇的模拟结果见表3。

为了评估结果的收敛性，引入收敛性$ {R_G} $，公式为：

$ {R_G} = \displaystyle\dfrac{{{\varepsilon _{k21}}}}{{{\varepsilon _{k32}}}} = \dfrac{{{\varphi _{k2}} - {\varphi _{k1}}}}{{{\varphi _{k3}} - {\varphi _{k2}}}}。$

(16)

式中：$ {\varphi _{k1}} $、$ {\varphi _{k2}} $和$ {\varphi _{k3}} $分别为高、中和低密度网格尺寸下的设置结果。根据$ {R_G} $的大小判断数值结果的敛散性，以$ {C_{aw}} $为参考值，求得$ {C_{aw}} $=0.592，满足0<$ {R_G} $<1,表明数值结果单调收敛。

3 水翼艇迎浪耐波性数值预报 3.1 波长对水翼艇波浪增阻的影响

为研究迎浪规则波中水翼艇的耐波性能，选取波高H=0.06 m，开展不同波浪工况下的数值模拟研究，具体的计算工况如表4所示。

图6为不同速度下水翼艇的运动与波浪增阻曲线。由图6(a)可知，当$\lambda /L < 1$时，不同速度下的纵摇响应幅值均较小；$F{r_\nabla }$=1.659时，纵摇幅值随着波长的增加而逐渐增加，当$\lambda /L > 2.5$纵摇响应趋于稳定，$F{r_\nabla }$=3.456时，纵摇幅值在$\lambda /L \approx 2.5$达到峰值。由图6(b)可知，不同航速下，垂荡运动均呈现先增大而后趋于稳定的趋势。纵摇与垂荡幅值均随着航速的增加而减小。

由图6(c)可知，不同速度的波浪增阻系数均随着波长的增加先增加后减小，当$0.5 < \lambda /L < 2.5$时，波长与艇长相近，艇体的响应较为剧烈，因而波浪增阻较大。当$\lambda /L > 2.5$时，波长远大于艇长，对于艇体相当于处于静水状态，波浪增阻随之减小。当$F{r_\nabla }$=1.659时，在$\lambda /L = 2$的波浪增阻系数最大，随着航速的提高，波浪增阻系数最大值出现在$\lambda /L = 2.5$处，${C_{aw}}$=9.485。

图6(d)为不同速度下艇体重心加速度幅值响应曲线，通过读取重心垂向加速度幅值${A_a}$并对其进行无量纲化，得到${A_a}L/(g{\xi _a})$。可知，不同航速下重心加速度响应趋势相似，均为随着波长的增加逐渐增大，经过峰值点后逐渐减小；当$F{r_\nabla }$=1.659时峰值点在$\lambda /L = 2$附近，当$F{r_\nabla }$=3.456时在$\lambda /L = 2.5$附近取得峰值，且高于$F{r_\nabla }$=1.659时的峰值。

3.2 水翼安装位置对水翼艇波浪增阻的影响

水翼的纵向安装位置容易影响水翼艇载荷分布，因此为研究水翼安装位置对水翼艇波浪增阻与运动的影响，本节对不同波浪工况下加装不同安装位置水翼的水翼艇进行数值模拟。

通过调整首部水翼位置，设计出3种不同水翼安装方案，首尾水翼的距离分别设置为2.36、2.46、2.56 m，如图7所示。

图8所示为$F{r_\nabla }$=1.659时不同首部水翼位置的水翼艇运动与波浪增阻曲线，并与未安装水翼的艇体进行对比。由图8(a)、图8(b)可知，安装水翼后，艇体的纵摇与垂荡在不同波浪工况下相较于未安装水翼的艇体均有所减小，说明安装水翼可以有效改善其耐波性能；在不同安装位置下纵摇与垂荡产生的影响规律基本一致，均为随着波长船长比的增加先增加后趋于稳定；当$\lambda /L < 2$时，水翼安装位置对纵摇与垂荡的影响较小，当$\lambda /L > 2$时，水翼安装位置对纵摇与垂荡的影响较为明显；首部水翼离尾翼越远纵摇幅值越大，主要由于当首水翼与尾水翼距离越大时，能产生更大的纵倾力矩，从而产生更大的纵摇。

由图8(c)可知，安装水翼在不同波浪工况下均可有效地减小艇体的波浪增阻；不同安装位置的波浪增阻的整体变化趋势基本一致，均为随着波长的增加先增大后减小；首尾翼之间距离越大波浪增阻越大，主要由于首翼的安装位置前移，导致首翼对艇体重心的力矩增大，艇体与水面的接触面积减小，导致船体的波浪增阻减小。

由图8(d)可知，当$\lambda /L = 2$时，各位置间的重心加速度幅值差距最大，当首尾水翼间距在2.56 m时，重心垂向加速度最小，而其他安装位置的重心加速度与未安装水翼的艇体相差不大。

4 结　语

1）水翼艇垂荡的幅值响应在不同速度下呈现先增大而后趋于稳定的趋势；纵摇与垂荡的幅值响应均随着航速的增加而减小；不同航速下，波浪增阻系数均随着波长的增加先增大后减小；在$F{r_\nabla }$=3.456，$\lambda /L = 2.5$，波浪增阻系数最大，${C_{aw}}$=9.485。

2）安装水翼在不同波浪工况下均可有效提高船舶的耐波性能。

3）不同水翼安装位置的波浪增阻系数均随着波长的增加先增大后减小，在$\lambda /L = 2$附近增阻系数达到最大；当首尾水翼间距为2.56 m时波浪增阻系数最小，纵摇幅值与垂荡幅值最大；在$\lambda /L = 2$附近不同水翼位置间的重心加速度幅值差距最大，当首尾水翼间距为2.56 m时，重心加速度幅值最小。在确保满足耐波性的设计基础上，适当首尾水翼的距离可行。