结合Savitsky 方法和重叠网格技术的滑行艇阻力数值计算与分析

郭军，扈喆，朱子文，陈作钢，崔连正，李贵斌

1 集美大学 轮机工程学院，福建 厦门 361021

2 上海交通大学 海洋工程国家重点实验室，上海 200240

3 上海交通大学 船舶海洋与建筑工程学院，上海 200240

4 中国船舶及海洋工程设计研究院 喷水推进技术重点实验室，上海 200011

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.228 6 m，型深0.143 m，排水量7.26 kg，底部斜升角为10°。滑行艇的艇体型线和几何模型如图1 所示。沿艇长方向，将滑行艇等分为20 站，从第2 站到20 站之间的型线相同，第0 站到第2 站之间的型线如图1（a）所示，图中数值表示型线图的站位。

图1 滑行艇模型Fig. 1 Geometry model of planing craft

滑行艇的无量纲速度为体积弗劳德数Fr∆：

式中，B为滑行艇艇宽。

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

表1 滑行艇的载荷工况Table 1 Load conditions of planing craft

2 数值计算方法

2.1 计算方法

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

连续性方程：

动量方程：

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

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

2.2 计算域和边界条件

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

图2 计算域和边界条件Fig. 2 Computational domain and boundary conditions

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

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

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

图3 计算域网格Fig. 3 Grid of computational domain

2.3 航态估算

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

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

图4 滑行艇的航态估算Fig. 4 Attitude estimation of planing craft

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

式中：Rt为滑行艇总阻力；Awet为滑行艇的动态湿表面积。

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

表2 航态预估对阻力预报的影响Table 2 Influence of trim estimation on resistance prediction

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 升沉计算值与试验值对比Fig. 5 Comparison of heaving between CFD and EFD results

图6 纵倾计算值与试验值对比Fig. 6 Comparison of trim between CFD and EFD results

图8 平均浸湿长度计算值与试验值对比Fig. 8 Comparison of average wet length between CFD and EFD results

3.2 压力分布

滑行艇压力系数CP的定义为：

式中：P为滑行艇表面的绝对压力；P0为参考点压力（1 个大气压）。

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

图9 不同载荷工况下的压力分布Fig. 9 Pressure distribution at different load conditions

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

图10 不同航速下的压力分布Fig. 10 Pressure distribution at different speeds

3.3 自由面波形

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

图11 不同航速下y = 0 剖面的自由面波高对比Fig. 11 Wave pattern of y = 0 at different speeds

图12 不同航速下的自由波面分布Fig. 12 Free wave surface distribution at different speeds

3.4 水气分布

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

图13 不同航速下的滑行艇表面的水气分布Fig. 13 Volume fraction of water of planing craft at different speeds

4 结 论

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

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

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

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

致 谢

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