高速滑行艇在规则波中的纵向运动数值研究

王硕，苏玉民，庞永杰，刘焕兴

(哈尔滨工程大学水下智能机器人国防重点实验室，黑龙江哈尔滨150001)

随着海洋权益的争夺和海洋滨海资源的开发日益被瞩目，滑行艇因其良好的机动性能和高速性能得到广泛的应用。不同于传统的排水型船舶，滑行艇在高速航行时依靠水动升力抬升船体，并伴有砰击、上浪和飞溅等现象发生。传统的船舶原理和理论已不适用于滑行艇的性能估算。需要新的数值方法对其水动力性能进行准确的分析和计算。

关于滑行表面的水动力计算研究起源于20世纪40年代NACA在兰利水池的拖曳实验［1］，基于这些实验后续的许多相关研究才得以展开。其中应用最广的是Savitsky于1964年提出的基于实验数据的半经验公式计算方法［2］。1997 年 Zhao等［3］提出了一种基于势流理论的二维半理论方法用于求解滑行艇水动升力问题。随后Sun等［4-6］应用该方法对滑行艇水动力问题进行了深入研究，但始终基于细长体假设。董文才等对深V型滑行艇进行了规则波迎浪纵向运动实验研究［7］，并且分析了波浪因素对滑行艇纵向运动的影响;又根据一系列的试验结果，提出了一种可预报滑行艇高速航行时的纵向运动的新方法，并且还通过模型实验验证了该方法的准确性［8］。近年来，应用RANS方程求解滑行表面水动力问题也已经不仅仅停留在研究阶段，CFD方法已经成为许多快艇设计者的参考工具。Azcueta［9］应用商用软件COMET对赛艇的水动力特性进行了计算，并应用于快艇的设计参考数据。苏玉民等［10］应用FLUENT软件对滑行艇在静水和波浪中的运动规律进行研究，定性的给出了滑行艇在稳定直航和迎浪航行中的运动规律。

滑行艇的耐波性能一直被认为是制约滑行艇发展的最主要问题之一。由于一般的滑行艇尺度相对波浪尺度较小，所以在波浪中的纵向运动是滑行艇耐波性研究中首要关注的重点。在高雷诺数下讨论滑行艇在波浪中的运动性能无论对于数值计算还是模型试验而言都是困难的。本文基于商业软件STAR CCM+中的RANSE VOF求解器，对高速滑行艇在规则波中的迎浪纵向运动进行了数值研究。首先选取Katayama［11］的小尺度滑行艇模型试验，对本文的数值方法计算精度进行了验证;然后给出了滑行艇在纵浪中的3种运动形式(无跳跃、规则跳跃和不规则跳跃)的运动特征，并讨论了波浪参数和航速对3种运动形式发生的影响。进而研究了跳跃的发生对滑行艇纵向运动函数的影响。

1 数值方法

1.1 控制方程

对于不可压缩的粘性流动，在笛卡尔坐标系下，用张量的形式表示的时均连续性方程和RANS方程可以写为

式中:i，j=1，2，3;ui为速度分量的时均值;为速度分量的脉动值为速度分量乘积的时间平均值;μ为流体的动力粘性系数。

数值模拟中采用隐式非定常形式，湍流模型为k-ωSST两方程模型［12］。计算区域应设置足够大以保证滑行艇周围流场不受边界影响。通常采用长方体计算域，船首到前方速度入口距离为1倍船长，舷侧到计算域边界距离为1.5倍船长，船尾流场长度为3倍船长。VOF(volume of fluid)法［13］是由 Hirt和Nichols在1981年提出来的，它是通过追踪流体域来跟踪自由面的，其基本原理是通过研究网格单元中流体和网格体积比函数F(0＜F＜1)来确定自由面，追踪流体的变化。

1.2 计算网格

自适应直角切割网格技术是近十几年来被广泛应用的一种网格划分技术，由于其生成简单，且兼具结构化网格的较高网格质量和非结构网格的复杂表面适应性，得到了广大CFD学者关注。切割网格技术起源很早，但由于其对边界的处理需要占用大量计算机时间，一直未得以广泛应用。直到1986年Keith Clarke等人利用直角切割网格技术和二维欧拉方程求解多段翼型，切割网格技术才重新发展起来。随后十几年间，切割网格技术得到比较快的发展。

切割网格生成的原理可简单概括为先在整个计算域周围空间生成体网格，体网格覆盖计算域内外，包含边界。然后用计算域边界面切割初步生成的体网格，删除计算域外的网格，进一步切割以提高计算域内网格质量。此种网格可以随意控制计算域内的网格疏密程度，大大减少了不必要的计算位置的网格数量，且网格质量较高，缺点是难以变形适应动网格技术［14］。

STAR-CCM+提供的 DFBI(dynamic fluid body interaction)模块是通过移动整个计算域来实现六自由度数值模拟的，计算中网格不需要变形，通过更新边界条件来实现滑行艇的运动。因此可以采用切割网格来实现滑行艇的运动数值模拟。划分网格时，需要对船体周围进行加密，绘制4层边界层网格，厚度设定为0.5 mm以满足各个航速下的壁面函数的需求。船体表面网格及计算域体网格见图1，为节省网格采用对称边界，对半船体进行计算［15］。

图1 滑行艇切割表面和切割体网格Fig.1 The trimming surface mesh and volume mesh of the planing craft

图2中给出了滑行艇表面边界层网格和自由液面加密层网格的绘制情况。在STAR-CCM+中采用棱柱层网格模型绘制边界层网格，可按规定网格层数和网格厚度绘制滑行艇表面网格。对于滑行艇在波浪中的运动数值模拟，波面周围的网格需要加密，尤其沿垂直方向在波高范围内需要10层以上网格，因此VOF模型才能准确捕捉到波面形状。

图2 滑行艇表面棱柱型边界层网格Fig.2 The prism layer m esh of the planing craft

1.3 运动模型

在建立滑行艇运动方程时，需采用2个坐标系，即滑行艇随体坐标系(动系)和固连于大地的固定坐标系(定系)。随体直角坐标系的原点位于滑行艇重心G处，坐标轴Gx、Gy、Gz分别经过重心G的水线面、横剖面和纵中剖面的交线。滑行艇的六自由度方程即为对于动系的质心运动定理和质心动量矩定理表达式为

式中:B为艇动量，Ω为角速度，F为外力合力，K为相对于质心的动量矩，V为航行速度，M为合力矩。

图3 滑行艇在波浪中运动自由液面Fig.3 The free surface of the planing craft sailing in waves

DFBI模块中内嵌的六自由度(6-DOF)求解器采用整个计算域随艇运动的动网格处理方式。即整个计算域内的网格没有变形和相对运动，而是将运动模型与流动求解器连接，整个计算域(包含船体表面和计算边界)按照船体运动方程的计算结果运动，在每一时间步更新控制网格的空间位置和VOF值以及边界条件函数，以保证自由液面位置在整体坐标系中保持不变。为及时更新边界，要求计算域来流方向和周边均为速度入口，仅船后方向为压力出口。每一时间步长根据船体的运动位移，并通过速度入口的两相流分布来更新自由液面位置(VOF体积分数)，来实现基于整体动网格的多自由度运动求解。图3给出了滑行艇在波浪上运动的数值模拟情况，波形被VOF准确捕捉并表达。

1.4 数值方法验证

为验证本文数值方法对滑行艇在波浪中的运动计算准确性，选取了Katayma开展的滑行艇迎浪运动模型试验进行验证性计算。首先建立该滑行艇模型的数值模型，船模型线及数值模型如图4。船模总长LOA=0.625 m，宽 B=0.25 m，型深 D=0.106 m，干舷d=0.059 m，船模重 W=4.28 kg，重心纵向位置LCG=0.285 m。实验中，船模拖点位于重心处，加速拖曳到实验要求速度，测量其在不同航速、不同波浪中的纵摇与升沉运动规律。

选取无跳跃与规则跳跃2个工况进行对比计算。

工况1:长度弗汝德数Fn=1.21，波高与干舷比Hw/d=0.68，波长与船长比λ/LOA=1.6，滑行艇未发生跳跃，有规律运动;

工况2:长度弗汝德数Fn=3.63，波高与干舷比Hw/d=0.51，波长与船长比λ/LOA=3.59，滑行艇发生规则跳跃。

图4 Katayama滑行艇模型型线图和三维模型Fig.4 Katayama’s planning craftmodel lines and 3-D model

图5 工况1、2中重心升沉曲线和纵倾角试验值与计算值比较Fig.5 Comparison of experiment and numerical data of heavemotions and pitch angles in case 1 and case 2

图5给出了2种工况下，纵摇和升沉的试验值与计算值对比情况。计算曲线符合较好，其中低航速下的精度明显高于高速时，但整体计算精度满足工程应用。该部分验证工作为本文前期工作，此处只给出简单计算结果，对于更多详情以及应用STAR-CCM+进行数值造波的验证请参考文献［16］。

2 计算结果分析

本文的计算对象，为自主设计的滑行艇模型，该艇型的型线及三维数值模型见图6，模型主尺度见表1。

图6 滑行艇模型型线及三维模型Fig.6 The body lines and 3Dmodel of the planing craftmodel

表1 船模主尺度Table 1 M ain dimensions of the ship model

2.1 滑行艇规则波中的3种运动情况

Katayama通过滑行艇耐波性模型试验将观察到的滑行艇迎浪运动分为3种情况:无跳跃运动、规则跳跃运动和不规则跳跃运动。本文通过数值方法进行模拟时，也成功观察到了滑行艇规则波中的迎浪3种运动方式。计算发现，这3种运动方式的运动规律不同，但其发生条件之间却没有明显的界限，随着条件改变，一种运动可以发展为另一种。例如，相同波浪参数下，随着航速的增加，无跳跃运动可逐渐发展成规则跳跃，规则跳跃可逐渐发展为不规则跳跃。本节首先选取了某一航速下(滑行状态Fn=1.93)3种运动形式发生的典型工况，对3种运动加以说明分析。

图7 无跳跃运动时滑行艇运动参数曲线Fig.7 The curves of motion parameters in no jump condition

图7给出了无跳跃运动下的滑行艇升沉、纵摇和重心处的垂向加速度(波长λ=4L，L为船长，波高h=L/40)。升沉和纵摇曲线均近似为正弦曲线，运动规律性明显。重心处的垂向加速度曲线也近似正弦分布，加速度峰值约为0.4g。滑行艇在未发生跳跃时，运动曲线近似为正弦分布，与常规船型在波浪中的运动规律无明显区别。

图8给出了规则跳跃运动下的滑行艇运动曲线和重心处的垂向加速度(波长λ=4L，波高 h=L/20)。发生规则跳跃时，滑行艇在波浪上的运动仍然稳定且呈周期性。升沉运动曲线仍近似正弦分布。而纵摇运动对比无跳跃运动时，在波峰与波谷间表现出不对称性。本文定义的纵摇曲线以艏倾为正，艉倾为负。图8中3条曲线的运动时间完全对应，可以结合3幅图，分析滑行艇的运动。例如图8(a)中升沉减小时，表明滑行艇在跳起后从最高点开始下落，此时对应图8(b)中滑行艇处于埋艏的过程，而图8(c)中加速度在经历短暂的自由下落-g状态变为受砰击加速度迅速上升的过程;当重心升沉增大时，船体被抬升向上跳跃，发生抬艏，纵摇角迅速减小，滑行艇逐渐脱离水面，加速度也由最大值逐渐减小为重力加速度。

由图8(c)可以看出，当滑行艇发生规则跳跃时，加速度的峰值接近重力加速度g，表明滑行艇已接近完全越出水面，与水接触面积很小。

图8 规则跳跃运动时滑行艇运动参数曲线Fig.8 The curves of motion parameters in regular jum p condition

图9给出了不规则跳跃运动下的滑行艇运动曲线和重心处的垂向加速度(波长λ=2L，波高h=L/20)。当发生不规则跳跃时，滑行艇的运动很不稳定，虽然从运动曲线上看仍然有周期性，但幅值和频率上均有明显变化。在现实情况下，如果有其他方向上的扰动和干扰，不规则跳动会进而演变成更极端的运动甚至倾覆。不规则跳跃一般发生在波浪遭遇频率较高的情况下，滑行艇经过一个波峰后，在下落过程中遭遇第2个波峰，由图9(c)可知，此时滑行艇受砰击最严重，向上加速度达到4g以上，在真实海况下这是极度危险的，即使在无人滑行艇上，工作设备也无法在这么大的垂向加速度上维持正常运转。在滑行艇实艇的设计和应用中，应避免不规则跳跃的发生。

图9 不规则跳跃运动时滑行艇运动参数曲线Fig.9 The curves of motion parameters in irregular jum p condition

2.2 跳跃发生条件

本文对波浪的数值模拟均基于一阶规则波，波浪参数包含波长与波高，滑行艇则考虑航速。图10给出了滑行艇在3个航速下、不同波高以及不同波长情况下的运动类型。

在较低航速时(Fn=0.96)，跳跃发生的可能性较低。如图10(a)，在4个波长(2L，4L，6L和8L)和5个波高(L/10，L/15，L/20，L/30 和L/40)组成的20 个工况下，仅有3个工况发生规则跳跃，不规则跳跃没有发生。在h=L/15时，仅在λ=6L时发生了跳跃;在h=L/10时，在λ=4L和λ=6L时发生了跳跃。这表明，在较大波高下跳跃容易发生，但跳跃的发生需要一定范围的波长，波长过短运动不能够充分发展而过长则波陡(h/λ)而跳跃也不会发生。

图10 不同航速下跳跃发生条件分布Fig.10 The distribution of jump at different forward speeds

当航速提高时(Fn=1.93)，如图10(b)，跳跃在更小的波高下便发生了。当λ=4L和λ=6L时，h=L/10、L/15、L/20和 L/30均发生规则跳跃，这表明伴随着速度增大，引起跳跃的波高界限大大降低了，高速下滑行艇更容易产生跳跃，但波长仍需要在一定范围内才会产生跳跃。值得注意的是，在λ=2且h=L/10时，发生了不规则跳跃。较大波高和较短波长的情况下发生的不规则跳跃表明，不规则跳跃可能是伴随着较大的波陡发生的。

随着航速进一步加大(Fn=2.89，图10(c))，跳跃运动在分布范围内大面积发生。只有在极小波高和在小波高大波长范围内没有跳跃发生。此时不规则跳跃的发生条件出现了较为明显的规律。在发生跳跃的范围内，不规则跳跃集中出现在波陡较大的区域。经统计可知，当波陡h/λ＞1/60时，均发生了不规则跳跃;在2个h/λ=1/60的工况中，有一个发生不规则跳跃;当h/λ＜1/60时，均未发生不规则跳跃。

综合对比图10中3个航速下的跳跃发生条件可知，高航速和较大的波高是跳跃发生的最基本条件，与此同时波长也需要在一定范围内，如果波长过长或过短，跳跃也不容易发生。在跳跃发生区间中，不规则跳跃的发生主要分布在较大的波陡范围内，即当滑行艇以一定航速遭遇到高且短的规则波时，就可能诱发不规则跳跃的发生。

2.3 跳跃对纵向运动幅度的影响

根据耐波性原理将滑行艇的运动进行无量纲化，定义ζ为滑行艇升沉运动幅度，ζW为波幅，定义ζ/ζW为滑行艇无量纲升沉运动幅值。图11和图12分别给出了λ=2L和λ=4L时的升沉运动幅度。θ为纵摇角幅度值，K为波数，纵摇运动无量纲幅值定义为θ/(KζW)，图13、14分别给出了λ =2L和4L时的纵摇运动幅度。图11～14中同时给出了曲线上各点对应的运动类型，以便分析运动类型对运动幅值的影响。

图11～14中升沉和纵摇的无量纲运动幅值给出了相同的运动幅值变化规律。由图中可知，当滑行艇未发生跳跃时，升沉和纵摇的无量纲幅值随波高增大均呈线性减小趋势。当规则跳跃发生时，运动幅值仍然呈下降趋势，但下降速度明显增大。由此表明规则跳跃的出现对滑行艇迎浪运动规律产生了影响，规则跳跃使滑行艇在经过波峰后跳起，越过了下坡面与波谷与下一上坡面与波峰接触，因此运动幅值较小。然而，在无规则跳跃运动中，滑行艇与波浪相互作用的相位是随机的，运动呈现不稳定性。当滑行艇发生不规则跳跃时，无量纲运动幅值随波高变化规律被彻底改变，下降趋势直接转变为上升趋势，并在达到峰值后又出现回落。

结合2.1中对3种运动随时间变化曲线的分析可知，当滑行艇从无跳跃进入规则跳跃时，运动规律只发生了小范围内的变化，从安全性角度而言，运动尚在可接受的范围内。但当滑行艇发生不规则跳跃时，运动规律发生较大变化且无规律可循，运动状态不稳定，垂向加速度表明砰击严重，是实际设计中应避免遭遇的海况。

图11 λ=2L时升沉运动幅度Fig.11 The amplitudes of heavemotionswithλ=2L

图12 λ=4L时升沉运动幅度Fig.12 The am p litudes of heavemotions withλ=4L

图13 λ=2L时纵摇运动幅度Fig.13 The amplitudes of pitch motionswithλ=2L

图14 λ=4L时纵摇运动幅度Fig.14 The amplitudes of pitch motionswithλ=4L

3 结论

本文基于RANSE VOF求解器对高速滑行艇规则波中迎浪纵向运动规律进行了数值研究，得到以下结论:

1)滑行艇在规则波中做迎浪运动时，其纵向运动可分为3种运动形式:无跳跃运动，规则跳跃运动和不规则跳跃运动，3种运动发生的条件与波高、波长及滑行艇航速相关;

2)滑行艇迎浪运动时，在一定波长范围内遭遇较大波高会导致滑行艇发生跳跃，随着航速增加，跳跃在更广的波长范围、遭遇更小的波高时发生;

3)滑行艇不规则跳跃的发生区间分布在波陡较大的波浪参数范围，随着航速增大，不规则跳跃运动更易于发生;

4)规则跳跃和不规则跳跃均对滑行艇迎浪运动无量纲幅值规律产生影响，规则跳跃对滑行艇迎浪运动规律影响较小，而不规则跳跃会彻底改变运动规律，不规则跳跃一般发生在较恶劣的海况下，此时滑行艇运动剧烈，稳定性较差，处于危险航行状态。

［1］CHAMBLISSD B，BOYD G M.The planing characteristics of two V-shaped prismatic surfaces having angles of deadrise of 20°and 40°［R］.Washington:Langley Aeronautical Laboratory:NACA.1953:28-76.

［2］SAVITSKY D.Hydrodynamic design of planing hulls［J］.Marine Technology，1964，1(1):71-95.

［3］ZHAO R，FALTINSEN O M，HASLUM H.A simplified non-linear analysis of a high-speed planing craft in calm water［C］//FAST’97.London，UK，1997:431-438.

［4］SUN Hui，FALTINSEN O M.The influence of gravity on the performance of planning vessels in calm water［J］.JEng Math，2007(58):91-107.

［5］SUN Hui，FALTINSEN O M.Hydrodynamic forces on a planing hull in forced heave or pitch motions in calm water［C］//22nd WWWFB.Plitvice，Croatia，2007:185-188.

［6］SUN Hui，FALTINSEN O M.Numerical study of planing vessels in waves［J］.Journal of Hydrodynamics，Ser B，2010，22(5):468-475.

［7］董文才，岳国强.深V型滑行艇纵向运动试验研究［J］.船舶工程，2004，26(2):14-16.DONGWencai，YUEGuoqiang.Experimental study on longitudinalmotion of deep-V-shaped planning craft［J］.Ship Engineering，2004，26(2):14-16.

［8］董文才，刘志华，吴晓光，等.滑行艇波浪中纵向运动理论预报的新方法［J］.船舶力学，2007，11(1):55-61.DONGWencai，LIU Zhihua，WU Xiaoguang，et al.A new theoreticalmethod on longitudinalmotion of planing craft in wave［J］.Journal of Ship Mechanics，2007，11(1):55-61.

［9］AZCUETA R.Computation of turbulent free-surface flows around ships and floating bodies［D］.Harburg:Technical U-niversity of Hamburg，2001:33-40.

［10］SU Yumin，CHEN Qingtong.Numerical simulation of a planing vessel at high speed［J］.Journal of Marine Science and Application，2012，11(2):178-183.

［11］KATAYAMA T，HINAMIT，IKEDA Y.Longitudinalmotion of a super high-speed planing craft in regular head waves［C］//4thOsaka Colloquium on Seakeeping Performance of Ships，Osaka，Japan，2000:214-220.

［12］张志荣，李百齐，赵峰.船舶粘性流动计算中湍流模式应用的比较［J］.水动力学研究与进展，2004，19(5):1-3.ZHANG Zhirong，LIBaiqi，ZHAO Feng.The comparison of turbulence models applied to viscous ship flow computation［J］.Chinese Journal of Hydrodynamics，2004，19(5):1-3.

［13］HIRT C W，NICHOLS B D.Volume of fluid(VOF)method for the dynamics of free boundary ［J］.Journal of Computational Physics，1981(39):201-225.

［14］CHUNG M H.Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape［J］.Computers＆ Fluids，2006，35(6):607-623.

［15］WANG Shuo，SU Yumin，ZHANG Xi，et al.RANSE simulation of high-speed planing craft in regular waves［J］.Journal of Marine Science and Application，2012，11(4):447-452.

［16］王硕，苏玉民，杜欣.滑行艇静水直航及波浪中运动的数值模拟［J］.华南理工大学学报，2013，41(4):119-126.WANG Shuo，SU Yumin，DU Xin.Numerical simulation of planing crafts sailing in calm water and in waves［J］.Journal of South China University of Technology，2013，41(4):119-126.