0.   引　言

滑行艇是依靠流体的动力升力将艇体托出水面而使其处于滑行状态的一种高速艇，具有航速高、机动性好等特点。凭借这些优越的性能，滑行艇已被应用于导弹艇、鱼雷艇、侦察艇、交通艇和游艇等船型，在军用和民用领域具有广阔的应用前景。滑行艇在整个航行过程中航态变化大、阻力变化剧烈，如何快速、准确地预报滑行艇的阻力性能，是研究人员普遍关心的问题。

早期的滑行艇阻力性能预报方法主要有半理论半经验公式方法（如Savitsky法^[1]等）、根据系列滑行艇试验资料整理得出的图谱方法，以及滑行艇阻力模型试验。其中，半经验半理论方法和图谱方法的预报精度受限于计算船型与母型船的相似程度，而滑行艇的模型试验虽然能够直观地反映出艇体的受力情况以及流场特性，但成本高、耗时长。随着计算流体动力学（CFD）计算精度的不断提高，采用CFD技术不仅可以预报滑行艇的受力及运动情况，还可以捕捉艇体周围的流场细节，是研究滑行艇水动力问题的一种经济、有效的方法。近年来，越来越多的研究人员采用CFD技术对滑行艇进行数值模拟，并取得了一定的研究成果。Brizzolara和Serra^[2]对棱柱型滑行艇流场的CFD数值模拟精度进行了研究，结果表明CFD技术可以用于滑行艇水动力问题的分析与计算。曹洪建^[3]基于FLUENT软件预报了滑行艇在固定姿态下的实船阻力，结果显示计算误差约为10%，认为采用FLUENT软件对滑行艇进行阻力预报可行，但对升力的预报则不理想。Ghadimi等^[4]采用RANS方法对滑行艇的姿态进行了数值模拟，结果显示RANS方法可以用于滑行艇的初步设计。马伟佳等^[5]采用混合网格对滑行艇的阻力进行了数值模拟，结果显示数值计算结果与试验结果有一定的误差。蒋一^[6]对某超高速三体滑行艇的快速性进行了数值计算，结果显示在速度较低时阻力的计算值与试验值吻合较好，但随着航速的增加，误差会逐渐增大，在最大航速下误差约为30%。Lotfi等^[7]基于 CFX软件对某断级滑行艇的水动力性能进行了研究，结果显示其阻力预报精度约为5%。Frisk等^[8]采用FLUENT软件和STAR-CCM+软件对滑行艇的阻力与姿态进行了数值计算，结果显示滑行状态下数值预报的阻力偏低，在半滑行状态下倾角的误差达32%。De Marco等^[9]采用重叠网格技术和变形网格技术对滑行艇的水动力进行了分析，发现重叠网格技术表现较好。孙华伟等^[10]采用STAR-CCM+软件探讨了船体表面第1层网格节点高度、船体表面网格尺度以及网格节点分布系数这3个因素对棱柱型滑行艇阻力计算精度、收敛速度及计算稳定性的影响。邵文勃等^[11] 采用STAR-CCM+软件探究了时间步长分析与近壁面网格划分对滑行艇阻力计算精度的影响，结果显示时间步长和近壁面网格划分均会对滑行艇的摩擦阻力与艇底的水气分布造成影响。魏子凡等^[12]采用NUMECA软件对翼滑艇进行了数值模拟，发现加水翼和防飞溅条能改善滑行艇的阻力性能。丁江明等^[13]研究了在RANS计算中网格类型、网格尺寸等因素对滑行艇艇底水气分布、兴波形状、艇体姿态以及阻力的影响，并与船模阻力试验结果进行了比较。易文彬等^[14]采用FLUENT，CFX 以及 STAR-CCM+软件探究了时间步长及网格数对预报结果的影响，发现采用CFX和STAR-CCM+软件模拟的滑行艇艇底表面有部分水体积分数约为0.5~0.8，艇底的水体积分数预报偏低，艇底水气分布异常，并认为较小的时间步长以及较细致的网格分布能够改善艇底水气分布的效果。李屺楠等^[15]采用重新建模法对滑行艇绕流场进行了数值模拟，结果显示重新建模方法相比重叠网格方法其计算精度和计算效率更高。王慧等^[16]采用重叠网格技术对滑行艇的兴波进行了模拟，结果显示采用CFD方法能够准确模拟滑行艇在静水中的航行姿态。

综上所述，虽然现在采用CFD技术可以预报滑行艇的阻力性能，但在半滑行状态和滑行状态下，滑行艇艇体与流体之间存在剧烈的相互作用，在数值模拟中常存在艇底水气分布异常以及喷溅区域较难真实模拟的现象，且当滑行艇处于高速航行状态时其阻力预报与模型试验之间仍存在较大的误差，因此需进一步研究半滑行状态和滑行状态下的滑行艇阻力高精度数值预报方法。为此，本文拟基于CFD方法，首先采用Savitsky方法对滑行艇的航行姿态进行预估，再结合重叠网格技术对滑行艇的阻力性能进行预报，以提供一种准确、有效的滑行艇阻力数值预报方法，并分析滑行艇航态预估对阻力预报的影响；然后，在3种载荷工况及6个航速下对滑行艇的阻力性能进行数值预报，并与试验结果进行对比；最后，分析滑行艇的流场特性。

1.   计算模型

本文将以Fridsma^[17]于1969年进行的滑行艇系列模型试验中的一种模型为研究对象。该艇长1.143 m，宽0.2286 m，型深0.143 m，排水量7.26 kg，底部斜升角为10°。滑行艇的艇体型线和几何模型如图1所示。沿艇长方向，将滑行艇等分为20站，从第2站到20站之间的型线相同，第0站到第2站之间的型线如图1（a）所示，图中数值表示型线图的站位。

图  1  滑行艇模型

Figure  1.  Geometry model of planing craft

下载: 全尺寸图片 幻灯片

滑行艇的无量纲速度为体积弗劳德数$F{r_\nabla }$：

式中：U为滑行艇航速；g为重力加速度；$\nabla$为排水体积。

按照$F{r_\nabla }$的不同，滑行艇的运动可划分3种典型的航态：排水型状态（$F{r_\nabla } \lt 1.0$）、半滑行状态（$1.0 \lt F{r_\nabla } \lt 3.0$）和滑行状态（$F{r_\nabla } \gt 3.0$）。

滑行艇的载荷系数${C_\nabla }$定义为：

式中，B为滑行艇艇宽。

滑行艇的不同载荷工况如表1所示。

表  1  滑行艇的载荷工况

Table  1.  Load conditions of planing craft

工况	载荷系数	吃水/m	重心纵向位置/m	重心垂向位置/m
1	0.304	0.026	0.400	0.067
2	0.608	0.042	0.457	0.067
3	0.912	0.057	0.457	0.067

下载: 导出CSV 

| 显示表格

2.   数值计算方法

2.1   计算方法

对滑行艇水动力性能的数值模拟采用商业软件STAR-CCM+，不可压缩流体的湍流基本方程由连续性方程和动量方程组成^[18]。

连续性方程：

动量方程：

式中：${u_i}$和${u_j}$为速度时均量；${x_i}$和${x_j}$为坐标；t为时间；$\rho$为流体密度； p为压力时均量；$\upsilon$为流体运动黏性系数；$- \overline {u_i'u_j'}$为雷诺应力。

本文的湍流模型采用SST k-ω湍流模型，该模型常被用于船舶水动力性能的数值计算，其表达形式为^[18]：

式中：$\mu$，${\mu _{\text{t}}}$分别为流体的动力黏性系数和湍流黏性系数；${\sigma _k}$，${\sigma _\omega }$分别为湍动能k和耗散率ω的湍流普朗特数；${\tilde G_k}$，${G_\omega }$分别为k和ω的衍生项；${Y_k}$，${Y_\omega }$分别为k和ω的耗散项；${D_\omega }$为ω的交叉扩散项。

在滑行艇阻力数值计算中，采用流体体积（VOF）^[19]方法和高分辨率界面捕捉格式（high resolution interface capturing，HRIC）对自由液面进行捕捉，求解器采用基于分离流的黏性求解器，压力与速度的耦合运用SIMPLE法，对流项离散格式为二阶迎风格式。采用动态流体−固体相互作用模型（dynamic fluid-body interaction，DFBI）模拟滑行艇的六自由度运动，在计算中具有升沉和纵倾的自由度，滑行艇的释放时间为0.5 s，缓冲时间为2.5 s。采用隐式不定常格式，国际拖曳水池ITTC推荐的船舶阻力数值计算时间步长$\Delta t$与垂线间长L_PP及航速U的关系为$\Delta t$ =(0.005~0.01)L_PP/U^[20]，本文针对滑行艇阻力计算采用的时间步长为0.005L_PP/U。

2.2   计算域和边界条件

由于模型几何及其绕流场关于中纵剖面基本对称，故数值模拟对象为滑行艇的半个绕流场。数值计算的计算域和边界条件如图2所示。计算域由静止的背景区域和随艇体运动的重叠区域组成，采用重叠网格法处理滑行艇高速航行时艇体姿态变化较大的问题。计算域的长度、宽度和高度分别为13L_pp，5L_pp和5L_pp，速度入口位于艇前2L_pp处，计算域顶部距离滑行艇1.0L_pp。将计算域入口、侧面、顶部和底部设为速度入口边界，计算域的出口设为压力出口边界，艇体表面设为无滑移壁面边界。在计算域的出口前设置消波阻尼以消除传至出口处的船行波，为避免波的反射，阻尼长度为3L_pp^{[21, 22]}。

图  2  计算域和边界条件

Figure  2.  Computational domain and boundary conditions

下载: 全尺寸图片 幻灯片

网格类型为切割体非结构化网格，在滑行艇周围和自由面等区域进行了网格细化。在滑行艇周围设置3个开尔文波系形状的网格加密区以进行局部加密，在设计水线处设置3个长方体的网格加密区以对自由面网格进行局部加密，在背景区域与重叠区域之间设置1个网格加密区以进行网格过渡，以使背景区域和重叠区域的网格尺寸处于相同量级。滑行艇的表面网格尺寸为1.75%L_PP ，滑行艇边界层网格使用棱柱形网格进行划分，共15层，增长率为1.2，边界层的厚度$\delta$为：

式中： L_W为滑行艇的平均浸湿长度；Re为雷诺数。

计算域的自由面网格和滑行艇表面网格的分布如图3所示。背景区域和重叠区域的网格量分别为90万和141万，最终的计算域网格总数为231万。

图  3  计算域网格

Figure  3.  Grid of computational domain

下载: 全尺寸图片 幻灯片

2.3   航态估算

滑行艇在排水型状态下的艇体姿态与静浮时差别不大，而半滑行状态和滑行状态下的艇体姿态则与静浮时差别较大。由于重叠区域内的网格是随滑行艇姿态的变化而做整体运动，此时艇体周围及自由面等区域需要更大的网格加密区，但这样会极大地增加网格数量，降低计算效率。因此在开始计算时，调整滑行艇的航态以使最终计算的航态接近于初始航态，这样不仅能降低网格数量，提高计算效率，还能改善计算的收敛性。此外，由于网格的方向与来流方向较为接近，故还能改善艇底的水气分布异常现象^[23]，从而提高计算精度。本文将采用半理论半经验公式Savitsky方法^[24-25]对滑行艇的航态进行估算。

式中：${C_{{\rm{L}}\beta }}$为升力系数；${F_{{\rm{L}}\beta }}$为升力（滑行状态时约等于重力）；$\beta$为斜升角；${C_{{\rm{L}}0}}$为零斜升角时的升力系数，

式中：$\theta$为滑行艇纵倾角；$\lambda _{\text{W}}^{}$为平均浸湿长宽比；$Fr_B^{}$为艇宽弗劳德数，$Fr_B^{}{\text{ = }}U/{(gB)^{0.5}}$。

式中，${x_g}$为重心距方尾的纵向距离。

开始计算时，在不同载荷系数下设置的滑行艇初始航态如图4所示。

图  4  滑行艇的航态估算

Figure  4.  Attitude estimation of planing craft

下载: 全尺寸图片 幻灯片

滑行艇的总阻力系数C_t定义为：

式中：${R_{\text{t}}}$为滑行艇总阻力；${A_{{\text{wet}}}}$为滑行艇的动态湿表面积。

在载荷系数${C_\nabla }$= 0.608及航速U = 4 m/s（$F{r_\nabla }$= 2.904）的工况下，对滑行艇进行数值计算，分析航态预估对滑行艇阻力预报的影响，计算结果与试验结果的比较如表2（表中，σ为滑行艇的深沉，L_wet为浸湿长度）所示。由表可见，与单纯采用重叠网格相比，结合Savitsky经验公式进行航态预估后，滑行艇航态、阻力和浸湿长度的预报精度均有所提高。

表  2  航态预估对阻力预报的影响

Table  2.  Influence of trim estimation on resistance prediction

方法	σ/B	θ/(°)	Ct	Lwet /B
试验结果 (EFD)	0.064	4.15	0.00682	3.41
重叠网格	0.058	3.94	0.00752	3.54
Savitsky+重叠网格	0.060	4.07	0.00725	3.47

下载: 导出CSV 

| 显示表格

3.   计算结果

3.1   阻力计算结果

在3种载荷工况及6个航速（U = 1，2，3，4，5，6 m/s）下，对滑行艇的阻力性能进行数值计算，并与阻力试验结果进行比较，以验证数值方法的准确性及有效性。

滑行艇升沉计算结果与试验结果的对比如图5所示。由图5可见，升沉的计算结果稍微偏小，3种载荷工况下升沉的平均误差分别为0.12%，0.31%和−5.06%。滑行艇纵倾角计算结果与试验结果的对比如图6所示。由图6可见，纵倾角的计算结果稍微偏小，3种载荷工况下纵倾角的平均误差分别为−0.62%，−2.11%和−5.86%。滑行艇总阻力计算结果与试验结果的对比如图7所示。由图7可见，阻力的数值计算结果偏大，3种载荷工况下阻力的平均误差分别为9.03%，4.09%和6.93%。滑行艇平均浸湿长度计算结果与试验结果的对比如图8所示。由图8可见，3种载荷工况下浸湿长度的平均误差分别为−0.63%，1.23%和1.61%。整体上，滑行艇的升沉、纵倾角、总阻力和平均浸湿长度的计算结果均与试验值吻合良好，说明本文所提的数值计算方法有效，数值计算结果准确可靠，可用于半滑行状态和滑行状态下滑行艇的阻力预报。

图  5  升沉计算值与试验值对比

Figure  5.  Comparison of heaving between CFD and EFD results

下载: 全尺寸图片 幻灯片

图  6  纵倾计算值与试验值对比

Figure  6.  Comparison of trim between CFD and EFD results

下载: 全尺寸图片 幻灯片

图  7  阻力计算值与试验值对比

Figure  7.  Comparison of resistance between CFD and EFD results

下载: 全尺寸图片 幻灯片

图  8  平均浸湿长度计算值与试验值对比

Figure  8.  Comparison of average wet length between CFD and EFD results

下载: 全尺寸图片 幻灯片

3.2   压力分布

滑行艇压力系数${C_{\text{P}}}$的定义为：

式中：$P$为滑行艇表面的绝对压力；${P_{\text{0}}}$为参考点压力（1个大气压）。

当航速U = 4 m/s时，在不同的载荷工况下，滑行艇艇体表面的压力分布规律如图9所示，其中图9（a）为滑行艇龙骨线的压力分布图，图9（b）为滑行艇艇体表面压力分布云图。图中，x为沿船长方向的位置。从中可以看出，龙骨线的压力系数在与水接触的位置显著增加，在驻点处时达到峰值，然后减小并达到稳定，在尾部逐渐减小至负值，至艉封板时达到最小；当${C_\nabla }$ = 0.304，0.608时，压力系数的峰值比较接近，导致滑行艇的纵倾角也比较接近；当${C_\nabla }$ = 0.912时，压力系数的峰值增加，压力中心位置前移，导致滑行艇的纵倾角增加。

图  9  不同载荷工况下的压力分布

Figure  9.  Pressure distribution at different load conditions

下载: 全尺寸图片 幻灯片

当载荷系数${C_\nabla }$ = 0.912时，不同航速下滑行艇艇体表面的压力分布规律如图10所示。从中可以看出，当$F{r_\nabla }$=1.357时，滑行艇处于半滑行状态的初始阶段，在喷溅区域其压力系数是先快速增加至最大，随后迅速减小至常压，然后缓慢增加，在艉部又迅速减小至负值；当$F{r_\nabla }$ = 2.714时，滑行艇处于半滑行状态的结束阶段，此时龙骨线压力分布的变化规律与滑行状态的变化规律基本相同；随着航速的增加，驻线压力系数的峰值减小，驻点线逐渐向艉部方向移动，驻点线与中纵剖面的夹角逐渐减小。

图  10  不同航速下的压力分布

Figure  10.  Pressure distribution at different speeds

下载: 全尺寸图片 幻灯片

3.3   自由面波形

当载荷系数${C_\nabla }$ = 0.912时，不同航速下滑行艇中纵剖面的自由面波高对比如图11（图中，z为自由面的波高位置）所示。滑行艇的自由面波形如图12所示，其中左图为整体波形，右图为喷溅区域波形。由图可见，随着航速的增加，船行波的角度逐渐减小，艏部兴波的高度逐渐降低，艇后鸡尾状水丘的高度逐渐增大，且逐渐向后移动，艇后“空穴”的深度也逐渐减小，而“空穴”的长度则逐渐增大；当$F{r_\nabla }$ = 1.357时，在滑行艇艏部发生了明显的喷溅现象，水流冲击滑行艇艏部后向前翻卷，喷溅的浸湿面积较大；当$F{r_\nabla }$ = 2.714时，在艇体底部发生了喷溅现象，喷溅区域的浸湿形状接近于三角形；当$F{r_\nabla }$ = 4.071时，喷溅区域的浸湿形状已完全为三角形。

图  11  不同航速下y = 0剖面的自由面波高对比

Figure  11.  Wave pattern of y = 0 at different speeds

下载: 全尺寸图片 幻灯片

图  12  不同航速下的自由波面分布

Figure  12.  Free wave surface distribution at different speeds

下载: 全尺寸图片 幻灯片

3.4   水气分布

当载荷系数${C_\nabla }$ = 0.912时，不同航速下滑行艇表面的水气分布如图13所示。由图可见，滑行艇艇底与水接触部分水相的体积分数接近于1.0，说明本文艇底水气分布预报正常，未出现艇底水气分布模拟异常的现象，这从另一个角度也说明了本文滑行艇阻力数值模拟的准确性。

下载: 全尺寸图片 幻灯片

  13  不同航速下的滑行艇表面的水气分布

  13.  Volume fraction of water of planing craft at different speeds

下载: 全尺寸图片 幻灯片

4.   结　论

本文在求解RANS方程的基础上，首先采用Savitsky方法对滑行艇的航行姿态进行了预估，然后结合重叠网格技术对滑行艇的阻力性能进行预报，并提供了一种滑行艇阻力高精度数值预报方法，接着，在此基础上对3种载荷工况下滑行艇的阻力性能进行数值预报，并与试验结果进行了对比，最后，分析了滑行艇的流场特性。通过研究，主要得出以下结论：

1） 随着载荷系数的增加，龙骨线压力系数的峰值增加，压力中心位置前移，滑行艇的纵倾角增加。

2） 随着航速的增加，龙骨线压力系数的峰值减小，驻点线逐渐后移，驻点线与中纵剖面的夹角减小；随着航速的增加，滑行艇船行波的角度逐渐减小，艏部兴波逐渐减弱，艉部鸡尾流的高度逐渐增大，艇后“空穴”的深度逐渐减小，但长度逐渐增大。

3） 本文模拟的艇底水气分布正常，能够模拟出滑行艇的喷溅现象，且滑行艇阻力数值计算结果与试验结果吻合良好，这些均说明本文数值方法准确可靠，可为滑行艇水动力数值研究提供有力的技术支撑。

致　谢

本文得到了中国船舶及海洋工程设计研究院喷水推进技术重点实验室的支持，在此深表谢意！