三体滑行艇纵向运动稳定性的数值模拟

邹 劲 王瑞宇 孙寒冰 蒋 一

（哈尔滨工程大学 船舶工程学院 哈尔滨150001）

引 言

高速三体滑行艇是常规滑行艇、高速多体船和气膜减阻船的组合船型，综合以上几种船型的优点，在高速滑行时，片体底部与水面接触，艇体产生的兴波与喷溅迅速被吸入槽道内，槽道内会形成气相区域、液相区域、气液二相混合区以及喷溅水层，在一定程度上减小了滑行阻力而且还提供了减震、缓冲的作用［1］，但也使其运动机理变得较为复杂，尤其是在高速航态下的“海豚运动”。目前，对于常规滑行艇纵向运动稳定性的理论判定方法是通过研究运动微分方程解的稳定性来判断系统的稳定性［2］，但是三体滑行艇底部复杂的结构使求解其在时域内任意时刻艇底的水流体动升力及力矩变得极为困难。而在CFD 领域，由于主要解决两相流问题，要求在自由液面附近的网格具有很高的分辨率，这一特殊性使商用软件在船舶六自由度运动具有较大航态变化的高性能船舶，如快艇、高速复合船方面的模拟具有局限性。然而重叠网格的出现使结构网格的自动化成为可能，也使得以上难题获得解决［3］，尤其是在非结构重叠网格技术推广到包含任意单元类型子网格重叠的情况之后，非结构重叠网格方法可以既具有非结构网格的灵活性，又具有结构重叠网格的计算精度［4］。

借助商业CFD软件STAR-CCM+，运用结构化背景网格与非结构子域网格相结合的重叠网格方法对某三体滑行艇在高航速下的“海豚运动”现象进行模拟，随后比照船模试验结果制定不同重心位置的计算工况，借助二分法对纵向运动稳定性界限曲线进行逼近，计算结果与试验结果相符，为类似船型的运动预报以及设计优化提供了一定参考。

1 船型简介

计算船型为某超高速三体滑行艇，主要参数如表1所示。图1为模型剖线示意图，缩尺比为1∶5.625，运用CATIA软件进行建模。

表1 模型参数

2 数值方法

借助商业CFD软件STAR-CCM+进行仿真计算，使用改进型的Realizablek-ε湍流模型来封闭RANS方程，该方法可较好地模拟存在流动分离和逆压梯度的复杂流动问题。自由液面处用VOF法进行追踪。

图1 模型横剖线示意图

图2 CATIA模型

运动模拟中采用自由模方法，释放船舶的升沉和纵摇两自由度，完全模拟船模在静水中的直航状态。船模在拖曳前进时，会对其周围流场产生影响，使压力场和剪应力场发生变化，也使船模所受的力和力矩发生变化：

船模刚体的六自由度控制方程为：

3 网格划分

采用重叠网格技术，将控制域划分为主域与子域两部分，在计算中每迭代一次都需重新记录交叉在一起的不同两种网格的信息并让数值在其中传递，虽耗时较长，却保证了计算对象网格形态的稳定，有利于实现姿态大幅变化的运动模拟［5-6］。

控制域尺寸为：宽为3L，入口距离船尾为3L，出口距离船尾为4L，气相、液相区各为1.5L。背景网格采用正交结构化网格，子域网格选择切割体网格（Trimmer），在自由液面处、船体运动范围内进行网格加密，且加密处网格尺寸与子域网格大致相当。背景网格基准尺寸为0.12倍船长，子域以及加密区域网格基准尺寸为0.012倍船长。

由于三体滑行艇航行速度极快，当越过阻力峰之后雷诺数Re达到107，因此壁面第一层网格厚度会在一定程度上影响计算机结果，边界层网格的最小尺度可由以下公式［7］求出：

式中：y+代表壁面无量纲距离；Δy表示壁面第一层网格厚度；L为船长；Re表示雷诺数。

以=4，Re=2.79×107为例，边界层层数设为三层，增长率取1.1，选取不同的y+并且生成四套网格，将平衡后得到纵倾角相对于试验值的误差作为依据，以此来确定网格方案。在表2中可以看到方案3误差较小，因此选择第三套方案。

表2 网格方案

图3 控制域

4 数值计算与结果分析

4.1 “海豚运动”现象数值模拟

常规滑行艇在静水中航行时，随着速度不断增加，会出现一种周期性、有界的垂直平面内的运动，这是一种垂荡和纵摇的耦合运动，称为“海豚运动”［8］。其产生原因主要是随着速度的增大，滑行艇尾部单位长度上的负荷逐渐加大。尽管垂荡及纵摇的小幅度变化只引起浸湿长度的微小改变，但是由于速度较大，所以整个水动压力的变化很大，使得滑行艇难以维持原有平衡。

下页图4为当xg= 0.561 5 m，速度v= 10.5 m/s时，不同时刻的艇底压力分布图。由图4（a）可以看到，t= 1.35 s时，水动升力集中在艇尾的一小部分区域，使浸湿长度减小，整体动升力下降，艇体开始下沉并伴随着埋首。当动浮力点不断前移并达到平衡时，由于惯性作用，运动并不停止，直到动浮力集中在船艏附近并形成小范围高压区，如图4（b）所示，此时艇体又开始抬首。这一过程循环往复，就形成了周期性的垂荡与纵摇组合——“海豚运动”现象。

图4 艇底压力分布图（xg = 0.561 5 m，v = 10.5 m/s）

4.2 纵向运动稳定性界限的求解

三体滑行艇的这种纵向不稳定性与 和重心位置xg有关。王庆旭等［9］通过修改重心纵向位置xg、改变航速v对该三体滑行艇的船模进行静水直航试验，并将得到的纵向失稳点进行曲线拟合，得到关于的曲线函数，即纵向稳定性界限曲线，其表达式为：

式中：CB=Δ/（ρV2B2/2），表示动负荷系数。

由式（6）可见，若单纯增加重心距离尾板的距离xg，失稳的临界会变大，也就是失稳速度会加大。为了验证数值计算结果，将数值计算工况与船模试验对应（文中所有实验数据均参照文献［9］），运用“二分法”对四种不同重心纵向位置处的失稳点进行逼近：速度相隔1 m/s，在试验值附近选取三个工况分别计算，若在所选工况中不存在未失稳点则将工况向下调整0.5 m/s；若包含，则对未失稳点和下一个与其相邻的工况进行二分继续计算，精度为0.5 m/s。初始工况如表3所示。

表3 数值计算工况

图5为垂荡时历曲线，下页图6为纵摇时历曲线。

图5 垂荡时历曲线

在图5（a）、6（a）中可以看到，xg= 0.561 5 m、速度在8.5 m/s时，纵向运动虽有规律，但幅度微小，升沉幅度不足0.01 m，纵摇幅度不到0.5°；当速度从8.5 m/s提升到9.5 m/s时，升沉以及纵摇幅度明显加剧，升沉幅度达到0.03 m，纵摇2.5°左右。不难看出，xg= 0.561 5 m时，纵向运动失稳点在8.5 m/s之前且就在附近。事实上，将速度降到v= 8.0 m/s时，“海豚运动”消失。同样，xg= 0.611 5 m时，v= 10 m/s依然失稳；当速度降到9.5 m/s时，海豚运动消失。因此，xg= 0.561 5 m、xg= 0.611 5 m时的临界失稳速度分别在8.0～8.5 m/s、9.5～10 m/s之间。而另外两个重心位置，在预先选取的最小速度点处均未出现“海豚运动”，对最小速度点与相邻速度点进行二分，继续计算，分别得到未失稳点v= 9.5 m/s、v= 12.5 m/s，将数值结果与实验结果比较，得到表5。

图6 纵摇时历曲线

表5 模拟结果与实验结果

在表5中，v1表示稳定速度上限、v2表示失稳速度下限。从计算结果中看到，随着重心位置的前移，纵向失稳速度变大，趋势与实验结果相吻合。与实验结果相比，xg= 0.561 5 m时的v2误差最大，为19%；xg= 0.731 5 m时的v2误差最小，为3.7%；在其他结果中，除xg= 0.561 5 m时的v1之外，误差均小于10%。图7为纵向稳定性界限。

表4 数值计算结果

图7 纵向稳定性界限

如图7所示，实曲线以上区域为试验得到的纵向失稳区域（“海豚区”），以下为稳定区，而通过模拟得到的稳定性上限曲线与失稳下限曲线的中间地带即为纵向运动稳定性界限曲线所在位置。不难看出与实验值相比较，模拟值的稳定区域较小，失稳区域较大，随着体积傅汝德数的增加误差减小，模拟值曲线更加贴近实验值。

5 结 论

本文运用重叠网格技术对三体滑行艇在高速航行时产生的“海豚运动”进行CFD仿真计算，并结合二分法对纵向稳定性界限进行逼近计算。结果显示：

（1）随着重心纵向位置的前移，失稳临界速度随之变大，趋势与试验结果相吻合。

（2）二分法逼近得到的纵向稳定性界限位于实验值下方，即失稳区大于实验结果。

（3）在高傅汝德数条件下，计算结果更加贴近实验值，误差极小。

综上所述，结果显示文中运用的数值计算方法可以较好实现对三体滑行艇纵向失稳运动的模拟，且失稳速度较为准确，相比于船模试验具有高效能、低成本等优势，可为类似船型的设计与优化等工作提供可信的参考数据。

［1］ 孙华伟.三体滑行艇船型与阻力性能研究［D］.哈尔滨：哈尔滨工程大学，2010：1-51.

［2］ Azcueta R，Caponnetto M，Soding H.Planing boats in waves［C］// Proceedings of IWSH’ 2003.Third International Workshop on Ship Hydrodyamics，Wuhan ：Wuhan University of Technology Press，2005，2.

［3］ 赵发明，高成君，夏琼.重叠网格在船舶CFD中的应用研究［J］.船舶力学，2011（4）：332-341.

［4］ 田书玲.基于非结构网格方法的重叠网格算法研究［D］.南京：南京航空航天大学，2008：56-80.

［5］ Thomas P D，Middlecoff J.F.，Direct Control of the Grid Point Distribution in Meshes Generated by Ellipitic Equations［J］.AIAA Journal.1980，18 ：652-656.

［6］ Royk，Bharats，Rajkeshar S A.Comprehensive Generalized Mesh System for CFD Application［J］.Mathe and Computers in Simulation，2008，78：605-617.

［7］ 邓锐.阻流板对双体船水动力性能影响的数值研究［D］.哈尔滨：哈尔滨工程大学，2010：31-38.

［8］ 赵连恩.高性能船舶水动力原理与设计［M］.哈尔滨：哈尔滨工程大学出版社，2007：195-199.

［9］ 王庆旭，邹劲，史圣哲，等.高速型三体滑行艇简介及纵向稳定性初步研究［A］.第十五届中国海洋（岸）工程学术讨论会论文集，2011：226-230.