滑行艇二自由度运动数值模拟方法

摘 要：利用商业计算流体力学软件CFX，基于RANS方程结合VOF法处理自由液面，建立了具有纵倾和升沉两个自由度的自由模拖曳的数值模拟方法。并用一条双断级滑行艇模型的阻力计算来验证了该方法的可信度。数值计算的结果与试验值相比，阻力值相差不大并且有相同的趋势，姿态的计算结果也再比较接近。建议在进行滑行艇自由面扰流问题的数值计算时，使用自由模方法模拟。

关键词：双断级滑行艇 两自由度 RANS方程 VOF方法

滑行艇的阻力性能是滑行艇最重要的性能之一，是滑行艇设计师在设计过程中所面临的首要问题，因此滑行艇的阻力预报一直都是滑行艇研究中的热点，并逐渐发展出了一系列的阻力计算方法。下面利用商业计算流体力学软件CFX为求解器来求解滑行艇的粘性流场，并利用其自带的六自由度刚体求解器（Rigid body Solver）结合动网格技术来求解滑行艇的运动姿态，建立滑行艇自由模模拟的数值方法，并用于一条双断级滑行艇模型的阻力计算。

数值计算方法

1、控制方程及运动方程

求解粘性问题就是求解N-S方程，本文对N-S方程的求解，采用RANS方程方法，对引入湍流模型后构成的封闭方程组求解得到湍流要素的时均值。连续性方程和动量方程为：

■ （1）

■（2）

式中，Ui，Uj为速度分量时均值（i，j=1，2，3）；P为压力时均值；ρ为流体密度，μ为动力粘性系数；为雷诺应力项。计入雷诺应力项后，本文计算中选取Shear Stress Transport（SST）湍流模型来封闭RANS方程。

2、湍流模型

SST模型结合了k-ω模型和k-ε模型的优点，能够较好的模拟存在流动分离和强逆压梯度的复杂流动问题。

k方程为：

■（3）

ω方程为：

■（4）

其中Gk和Gω为由于平均速度梯度引起的湍流动能的产生，Yk和Yω为关于k和ω的湍流耗散项。

3、自由液面的处理方式

自由液面的捕捉采用VOF法，设计算区域是V，流体1所在的区域记为V1，而流体2所在的区域记为V2。定义函数

■（5）

对两种不相容的流体组成的流场， 满足：

■（6）

其中=（u，v，w）为流体的速度场，定义VOF函数Cijk为在网格单元上的积分除以单元体积，即：

■（7）

并且Cijk满足：

当C=1时，网格充满流体1；C=0，网格不含流体1；0

4、运动方程

滑行艇模型航行时所受到的力和力矩分别为

■（8）

■（9）

在求解过程中艇模所受到的力和力矩由以上两式得到，其中为剪切应力，p为压力，为艇模外法线方向，m为艇模的质量，为艇模重心处的位置矢量，下标C代表艇模的重心即转动中心，为艇模表面任意一点的位置矢量。

艇模的六自由度控制方程为：

■（10）

■（11）

通过控制方程可得到艇模的运动姿态，其中，为艇模的平动速度矢量，是惯性矩张量，为艇模的角速度矢量。

在计算过程中通过每次迭代所得到的力和力矩来解算艇模的姿态，并实时的更新艇模表面网格的节点位移来实现艇模的运动模拟。

船池模型及网格划分

1、船池模型

本文计算的双断级滑行艇如图1所示，艇底沿纵向设置两个断级，断级将艇底分为三个滑行面，断级线在铅锤方向的投影与船模中心线的夹角为80°。

图1 双断级滑行艇模型

计算网格采用结构化网格，并在艇体附近进行加密，计算网格如图1，图2所示。令船长为L，则所采用的计算域为：沿船长方向向上游延伸1.5L，向下游延伸4L，沿船宽方向向外侧延伸1.5L，垂直方向向上延伸1.5L，向下延伸1.5L。由于断级处结构比较复杂，在断级处进行了网格的加密。

图2 计算域网格 图3 滑行艇表面网格

2、边界条件

入口：指定来流速度，即为船模试验中各点所对应的船速

对称面：采用对称边界条件

出口：出口处指定压力分布为静压

船体：不可滑移壁面

其他壁面：滑移壁面

计算结果及分析

1、典型的收敛时历

在本文的数值模拟中，艇模的运动包含两个自由度，即升沉和纵倾。图4到图6给出了在特定的航速下艇模的阻升力及纵倾的收敛时历图。

图4阻力及升力的收敛时历图 图5 纵倾收敛时历图

图6 升沉收敛时历图

2、计算值与试验值的比较

从图4到图6可以看出，由于初始时刻滑行艇姿态为任意姿态，而非一个稳定的平衡位置，因此在初始升力和阻力都会发生较大的变化，同时姿态也会发生较大的改变；经过一定的时间步之后，阻升力逐渐收敛至比较稳定的状态，即在某一个值附近微幅震荡；相应的纵倾也会相对稳定，此时所对应的艇模的姿态即为自由模前进的平衡姿态。

图7 阻力计算结果 图8纵倾计算值

图9 计算升沉值与试验值

从图7到图9中可以看出，数值计算所得到阻力、升沉、纵倾有一定差别，但总体上符合较好。计算阻力值在与试验值相比整体偏小，随着速度的增加其误差也逐步增大；而纵倾与升沉值则在中速段与试验值契合较好，在低速及高速段则有一定的偏差。

由于低速时艇模处于排水航行状态，此时姿态的变化对阻力性能的影响不如高速时明显，因此此时的计算阻力值仍然能够保持一定的精度；而高速时由于升沉值的计算结果较试验值偏大，因此计算时艇体的抬升将更严重一些，造成高速时的阻力偏差较大。

3、计算结果的后处理

图10到图13给出了航速为8m/s和10m/s时的艇模绕流场的后处理结果，从图上可以看出，后处理结果可以很好的反映流场的细节，从压力分布图上可以看出，随着航速的增大，艇底的动压力也逐渐增大，并且压力中心也逐渐后移；从自由表面的波形图上可以看出，随航速增大，艇体的抬升，滑行艇两侧的兴波也逐渐减弱，艇后的鸡尾流也向后移动，这与试验现象是相符合的。

图10 8m/s时艇底的压力分布 图11 10m/s时艇底的压力分布

图12 8m/s时自由液面波形 图13 10m/s时自由表面波形

结论

根据本文的艇模的计算结果，表明根据采用本文的数值计算方法，能够较好的模拟滑行艇自由模拖曳前进的情况：①各速度点所对应的阻力计算结果，与阻力试验值吻合较好，但由于高速时升沉值的差别较大，造成高速段的阻力值存在一定的误差。②从后处理结果上可以清晰的看出流场随航速的变化情况，这在试验中是难以观察的，对试验来说是一个很好的补充。

因此在处理滑行艇等高速或新船型的运动模拟问题上，上述方法具有较好的适用性，尤其是在新船型的开发中可以代替模型试验，以大量节省时间和人力，大幅缩短设计周期。

参考文献：

[1]董文才，郭日修.滑行艇阻力研究进展[J].船舶力学，2000，4（4）：68-81页.

[2]Azcueta R.，Rousselon N. CFD APPLIED TO SUPER AND MEGA YACHT DESIGN. International Conference on Super &. 20th-21st May 2010.

[3]Rodrigo Azcueta.STEADY AND UNSTEADY RANSE SIMULATIONS FOR PLANING CRAFTS.The 7th International Conference on Fast Sea Transportation . 2003： 41-48P.

[4]Caponnetto.M.Numerical Simul-ation of PlaningHulls.3rd NUTTS Conference，2000.

[5]Caponnetto M.Sea keeping simu-lation of fast hard chine vessels using RANSE.5th NUTTS，2002.

[6]Caponnetto M.，Soding.H.，Azcu-eta.R..Motion Simulations for Planing Boat in waves，Ship Technology Research，2003.

[7]王福军.计算流体动力学分析―CFD软件原理与应用[M].北京：清华大学出版社，2004.66-76.

（第一作者单位：南宁船舶检验局）