刚性气封局部气垫双体船波浪运动数值分析

杨静雷 孙寒冰 李晓文 扈喆

（1.集美大学轮机工程学院，福建厦门 361000；2.哈尔滨工程大学船舶工程学院，黑龙江哈尔滨 150001）

局部气垫双体船（PACSCAT）是一种以细长型双体船为基础，辅以局部气垫作为支撑的新型高性能船舶［1］。它不仅具有较好的快速性，而且具有良好的航行稳定性和操纵性，能够以滑行艇的姿态高速航行。由于局部气垫双体船拥有宽大的侧船体，因此可以运用螺旋桨或者喷泵作为推进装置，使该船型可以拥有比全垫升气垫船更好的操控性能。

局部气垫双体船的两侧片体的底部宽大扁平，这种船型特点有利于船体的冲滩、坐滩和退滩。气垫位于两侧片体和艏艉柔性气封装置之间，由于侧片体较为宽大，其气垫排水体积占全船排水体积的30%左右［2］。因此，影响运动性能的主要因素除了船体本身的特性外，气垫作用也会直接影响到内部的自由液面，进而会影响到总的运动性能。周伟麟等［3-4］对具有简单剖面形状的气垫船侧体采用试验数据拟合出不同遭遇频率下的二维切片水动力系数，但受制于船型因素，仅当气垫船侧体的剖面形状与试验船型相似时才可以根据该方法得到相对准确的水动力系数。三维势流理论由于可以考虑流场的三维效应，目前在船舶水动力学领域得到了广泛应用，其中Lee等［5-7］均使用频域零航速格林函数法或基于该方法的计算软件针对零航速气垫平台侧体的水动力性能进行研究。

在气垫船的运动仿真计算中，除了风机特性方程由于无法完全模拟风机打气外，其他的力学模型都是非线性，运动特性也是强非线性［8-9］，这正好与局部气垫双体船由于多重因素引起的强非线性现象有着高度的吻合。而计算流体力学（CFD）方法由于考虑粘性对流场的影响，可以更为精确地模拟非流线形物体的三维绕流场，理论上更适合模拟有航速船舶在波浪中的运动。其中尤忠强等［10-13］均采用数值计算方法开展过船体总阻力等方向的研究工作，并基于内部流场影响因素的分析，研究了内部流场特性对气垫船船体总阻力的影响。本文作者在对局部气垫双体船的研究工作中，基于试验方法和CFD 方法开展了大量细致的研究［14］，尤其是通过CFD方法对气垫系统各要素的影响因素对该船型在波浪中的运动特性进行了分析。所不同的是，受制于当时研究水平的限制，对研究模型和气动参数进行了较多简化，造成研究结果具有一定的误差。本文正是在前面研究的基础上，改进简化的计算模型，建立与试验模型一致的气封结构，基于不同的网格形式，在刚性气封数值计算方法的基础上，对局部气垫双体船的静水阻力计算精度进行研究分析，进而在优化的网格基础上，开展其在波浪中的运动特性研究。

1 数值计算方法

1.1 控制方程及理论

对于传统的计算流体力学，数值计算对于不可压缩的粘性流动和多相流理论，都会遵循工程应用中最广泛的N-S 方程［15-16］，忽略密度脉动的影响，考虑平均密度变化，表达如下：

式中：ui、uj为速度分量的时均值；ρ为流体密度；μ为流体动力粘性系数；ui'、uj'为速度分量的脉动值；p为压力的时均值；Si为动量方程的广义源项。

在运动界面追踪问题的数值模拟方法中，本文采用的是VOF自由液面法［17-18］，其特点是将运动界面在空间网格内定义为一种流体体积函数，并建立这种流体体积函数的方程。

本文中采用边界流体速度函数直接输入法进行造波［19］，根据波浪水质点的速度和波面方程定义入口处的速度来实现。该方法对于波浪的波长、波幅和频率等均能有较好的控制，数值模拟收敛性更好。对于目标波浪为有限水深平面线性正弦规则波，有关波浪理论方程表达如下。其中水平方向速度方程为

垂直方向速度方程为

式中，ζa为波幅，k为波数，ω为圆频率，dw为水深，c为速度，t为时间。

波浪周期定义为

波长λ定义为

有限水深条件下波浪周期和波长的关系为

式中，λ为示波长，g为重力加速度，T为波浪周期。

1.2 数值计算模型

局部气垫双体船实验模型艏艉气封采用传统的柔性装置来封闭空气，从而达到垫升减阻的目的。其中艏部采用囊指式围裙，艉部采用双囊式围裙，这种构型除了能保证气垫良好的密封性之外，还能提高局部气垫双体船的越障性能和让浪性能。两舷设计为宽大片体，艏部结构类似于滑行艇［2］的滑行壁面，连接桥为浮箱式结构用作增压室，局部气垫双体船的船型示意图如图1所示，计算模型主要尺度如表1所示。

表1 模型主要参数Table 1 Main dimension of ship model

图1 局部气垫双体船的船型特征Fig.1 Hull characteristics testing of PACSCAT

当气垫系统工作时，气垫系统的进气形式采用增压室式进气的方式，高压气体从风机打压后直接进入到增压室内，一部分经整流后通过增压室内一排气孔直接进入气垫舱室内，另一部分高压气体经气囊注入气垫舱室内。其垫升系统气体输送原理如图2所示。在对气垫系统进行数值模拟时，气体进气口设置在增压室顶部，由此设定恒定的气体流量向增压室内注入高压气体，气体通过与气垫舱室联通的气孔注入，数值计算的进气形式与实验保持一致。

图2 数值计算中垫升系统各个气孔的设置Fig.2 Hole’s setting of the air cushion system in the numeri⁃cal calculation

本文的艏艉气封均设置为刚性结构，刚性艏艉气封形式除了结构是固体的，其与船体的连接也是一体的，即作为船体刚性体的一部分随船体一起做六自由度运动，刚性艏艉气封形式如图3 和图4所示。

图3 刚性艏气封结构Fig.3 Rigid bow seal structure

图4 刚性艉气封结构Fig.4 Rigid stern seal structure

鉴于局部气垫双体船的船型特点，计算域采用整体型式以保证气垫系统的完整性。对于数值计算流场的模拟，计算域的设置尽量保证计算效率及计算精度的有效平衡，并保证消波区至少为波长的一倍以上，这里对不同波长采用了不同的计算域，其中短波2、5、7 m计算域的设置为：船前1倍船长，船尾3 倍船长，船侧1 倍船长，自由液面以上为1倍船长，自由液面以下1.5 倍船长。长波9、12 m计算域的设置对船尾消波区进行了延长，尺度为5倍船长，其他参数设置不变。

计算域入口边界条件为速度进口，出口边界条件为压力出口，顶部、底部以及两侧设置为滑移壁面。计算模型作为六自由度刚体运动位于自由液面附近。如图5计算域和边界条件情况。数值波浪采用VOF方法来捕捉自由液面的运动情况，波浪型式为一阶Stokes波。计算中流场为欧拉多相流，采用隐式非稳态进行连续体迭代。数值计算采用与试验完全一致的波浪参数，即规则波波幅为25 mm，计算波长分别为2.0、5.0、7.0、9.0、12.0 m，计算航速为3.6、5 m/s。

图5 数值计算域Fig.5 Numerical computational domain

1.3 数值计算网格

本文中分别采用滑移网格和重叠网格，来对局部气垫双体船在不同航速下的静水阻力进行数值计算，判断这两种网格形式对于研究此类船舶的适用性。对于滑移网格划分除了对气垫系统进行额外的网格加密处理外，其他的网格加密设置与传统数值计算方法类似，即在计算模型周围以及自由液面处进行网格加密设置，加密区域分别为BLOCK1、BLOCK2、BLOCK3.其中BLOCK1 和BLOCK2 加密网格基数设置分别为船体周围网格设置的最小值和最大值，而BLOCK3为总体计算域网格设置的最小值。重叠网格对于计算模型的加密区域和运动网格类似，但是需要在船体周围设置额外的流场域，该流场域边界条件与大的流场进行重叠网格设置。图6为两种网格形式的划分示意图。

图6 数值计算域网格加密形式Fig.6 Numerical calculation domain grid encryption form

需要说明的是，BLOCK1 加密区域为局部气垫双体船增压室和气垫舱室，包含艏艉刚性气封结构。计算中发现，由于气流从增压室进行打压和整流后进入气垫舱室内，该区域网格的加密是为了能够更好地捕捉该高压气流场的注气过程，防止因紊流导致部分气体能量缺失进而引起气垫压力和增压室压力的计算失真。除此之外，虽然艏部囊指进行了刚性结构处理，但不同裙指间存在极为细小的间隙，需要在此处加密网格数量，保证两个裙指间至少2个网格数量。

BLOCK3 为自由液面网格加密区域，设置的波浪波高为25 mm，短波和长波的波长船长比具有较大差异。而为了更好地捕捉波浪的特性，本文对不同波长下的网格采用了不同的划分形式，以保证在整个计算域内，单位波长下的网格数不少于70个网格，单位波高下不少于20个网格。如图7所示为波长船长比（λ/L）为1.67时的自由液面网格划分示意图。

图7 λ/L=1.67时自由面网格划分示意图Fig.7 Mesh scene for the free surface with λ/L=1.67

一般来说，数值计算的界面捕捉和边界层模拟问题对于精度和计算效率有较大的关联。从历来研究人员大量的试验和数学分析中可知：紧贴近固壁区域的流体可以分为粘流底层和对数律层。粘性底层内分子的粘性效应比较明显，而湍流粘性效应相对较小；在粘性底层之外的对数律层内流体的无量纲速度及温度分布服从对数分布律；粘性底层和对数律层之间存在着一个过渡区，该过渡区内分子粘性效应与湍流粘性效应相当。如果距离固壁表面的第1层网格节点距离这个区域较远，那么计算所得到的结果将是不准确的，也就是说，固壁表面以外第1层计算网格节点所在的位置对物体受力的计算影响较大。因此，对于船体附近的网格，船体表面以外第一层网格的分布非常重要，尤其是第一层网格节点的高度，是影响计算结果的重要参数。通常定义船体表面以外第一层网格节点到船体表面的无因次距离为y+，则无因次参数y+可以给出如下计算公式：

式中，y为第1层网格节点距离船体表面的高度，L为几何体长度，Re为相对几何体长度所定义的雷诺数。邓锐［19］在对双体船的阻力进行数值模拟时得到结果，证明y+值大约处于112.2 ≤y+≤281.2范围内时，第1 层网格节点位于过渡层附近，这部分节点的分布有利于捕捉船体表面附近流场中的物理量及物理量的梯度。特别提到的是，结合文献，对于3 m 左右的双体船船模而言，在实际计算时，第1 层网格节点距固壁表面的y+值取100 ≤y+≤200 就能够满足计算精度的要求。本文所研究的局部气垫双体本身拥有两个宽大的侧片体，严格意义上来说也拥有高速双体船的共同特性，尤其是研究的模型在3 m 左右，因此，对于边界层网格节点的划分和计算，严格借鉴了邓锐［19］的有关结论。同样采用了3 种不同尺寸的边界层参数进行对比计算，边界层设置为5 层，相对增长率为1.5，第1 层网格的基本尺寸分别为0.000 94、0.001 17、0.001 40 m，满足y+值范围在（100，200）之间的要求，边界层的计算工况设置如表2、表3所示。

表2 v=3.6 m/s静水性能计算工况Table 2 Hydrostatic performance calculation conditions at v=3.6 m/s

表3 v=5.0 m/s静水性能计算工况Table 3 Hydrostatic performance calculation conditions at v=5.0 m/s

2 计算结果及分析

2.1 静水中船体阻力计算结果及分析

对局部气垫双体船在波浪中的运动结果进行分析前，文中就静水中的阻力特性进行了计算分析，以判断该数值计算方法的精度以及存在的问题，数值计算的工况与波浪水池试验的参数一致，船重为145 kg，风机流量为150 m3/h，计算的速度点分别为3.6和5.0 m/s。数值计算域的设置与短波计算域的设置一致，即船前1 倍船长，船尾3 倍船长，船侧为1 倍船长，自由液面以上为1 倍船长，自由液面以下为1.5倍船长，VOF 计算模型的设置为平静水面。由于局部气垫双体船气垫系统的存在，高压气体和船体高速航行状态双重作用会导致气液两相流的复杂扰动，液面的破碎会带来较为困难的界面捕捉。

表4 和表5 分别为使用重叠网格和滑移网格得到的不同速度节点下数值计算值和实验值的结果，分别用Rcal和Rexp表示，其中计算误差用e表示，从中可以看出，重叠网格相对于运动域网格在计算具有刚性气封结构的局部气垫双体船的静水阻力时误差值更小，尤其是在高速v=5.0 m/s，重叠网格相对于运动域网格在两个速度节点下均有更小的误差。

表4 v=3.6 m/s静水性能计算结果Table 4 Clam water performance calculation result at v=3.6 m/s

表5 v=5.0 m/s静水性能计算结果Table 5 Clam water performance calculation result at v=5.0 m/s

2.2 波浪中船体阻力计算结果及分析

本节首先对基于数值造波理论的波浪情况进行了计算，以判断数值波浪衰减特性能否满足计算的需要。计算网格采用重叠网格，计算速度节点为3.6 m/s，波长为5 m、波幅为25 mm，选择的计算网格状态为表4 中的C3 形式。文中在计算域波面y=0 m 的纵剖面上平行设置波浪监测点，波浪检测点的位置沿x方向分别距离速度入口x=0.5λ、x=λ和x=1.5λ。图8为数值计算波高历时曲线。本次算例为计算时长达到35 s后的波形，从中可以看出的是，波浪传播一定时间后基本达到稳定状态，稳定后历时曲线波高产生一定的衰减，而谷值相对于峰值的衰减更大。主要原因是波浪峰值处为密度较小的气体，气体黏性小，耗散小，而谷值处为液体，密度大、黏性大从而导致耗散更大。因此，从波浪检测点看波浪虽有一定的耗散，但波浪基本形态比较稳定，总体上起伏不大，这里认为该造波方法可以应用于局部气垫双体船波浪运动的研究。

图8 数值计算波高历时曲线Fig.8 Numerical calculation of wave height duration curve

对于波浪中的阻力特性，图9为速度v=3.6 m/s时，波长分别为5.0、7.0、9.0和12.0 m状态下阻力计算值的历时曲线和试验值的对比，其中exp 表示为试验数据，cal1表示为刚性气封状态下的计算数据。从中可以看出的是，阻力的计算值无论是幅值还是均值都大于试验值，尤其是在波长为7.0 m处，两者的误差值达到最大，此后随着波长的增加误差值相对减小。从历时曲线的特征规律看，数值计算的结果非线性特性更加突出，各个阻力历时曲线没有表现出有效的规律特征，这也与试验结果形成较大的反差。造成这种较大误差的因素除了上文中静水阻力计算中存在的一定数值计算误差外，另一个主要原因是刚性气封结构在波浪中由于无法产生变形而产生较大的阻力增量，从不同波长下的情况可以看出，这种阻力增量在波长λ在7 m 时达到最大；然后随波长的增加刚性气封产生阻力的值相对减小，尤其是在波长达到12 m 时，由于波长船长比达到了4，波浪相对振荡频率大幅降低，刚性气封的推水阻力大幅减小，从而使得总阻力在波浪中的特性更接近试验值。这种情况与在高速v=5.0 m/s 时各个波长下的数值计算结果较为相似，如图10 所示。因此，对于刚性气封形式下的阻力数值计算，各个波浪中的历时曲线所产生的结果差异与波长有关，主要原因是刚性气封产生瞬时推水阻力。

图9 v=3.6 m/s时不同波长阻力计算值与试验值的对比Fig.9 Comparison between calculated and experimental values of the resistance with different wavelengths at v=3.6 m/s

图10 v=5.0 m/s时不同波长阻力计算值与试验值的对比Fig.10 Comparison between calculated and experimental values of the resistance with different wavelengths at v=5.0 m/s

2.3 各运动参数结果分析

图11～图13分别为速度3.6 m/s时，垂荡H、纵摇θ以及船舯加速度am在不同波浪情况下的计算值和试验值的对比。从图11 看出，垂荡运动计算的历时曲线和试验值相比具有一定的相似性，尤其是在波长7.0 m 时，垂荡运动在幅值和均值方面和试验值虽有误差，但没有如阻力一样展现出复杂的非线性特征，同样的规律在纵摇运动的对比曲线中可以看出，如图12 所示。图13 中加速度的计算值和试验值相比有较大的不同，计算值相对来说更有规律性，与试验值更为突出的非线性特征有所不同。总的来看，垂荡运动和纵摇运动的计算历时曲线与试验结果有一定的相似性，较为接近的波长段是在7.0 m附近。

图11 v=3.6 m/s 时不同波长垂荡运动计算值与试验值的历时曲线对比Fig.11 Time curve comparison between calculated and experimental values of the heave motion with different wavelengths at v=3.6 m/s

图12 v=3.6 m/s 时不同波长纵摇运动计算值与试验值的历时曲线对比Fig.12 Time curve comparison between calculated and experimental values of the pitch motion with different wavelengths at v=3.6 m/s

图13 v=3.6 m/s 时不同波长船舯加速度计算值与试验值的历时曲线对比Fig.13 Time curve comparison between calculated and ex⁃perimental values of the midship acceleration with dif⁃ferent wavelengths at v=3.6 m/s

图14～图16为各个运动参数的幅值随波长和船长比（λ/L）的变化曲线对比情况，将各运动参数作无因次化处理。从图14 可以看出，垂荡运动在短波λ/L小于2.0 时计算值和实验值误差较小，长波时误差明显增大，低速工况下的长波阶段误差要大于高速工况。而纵摇运动的对比曲线则不同，短波时误差较大，长波误差较小；在v=3.6 m/s 工况时的误差要大于v=5.0 m/s 时。船舯加速度的计算对比和垂荡运动相似。纵摇运动的不同主要原因来自于刚性气封的作用，在短波中航行时，相当于船体在遭遇高频阶段运动，刚性气封在此阶段内触水频率最大，受到水的反作用力也更加频繁，导致船体姿态没有完成调整后又向下运动，从而引起计算误差偏大，这种现象在高速时更明显。

图14 垂荡运动计算值与试验值的幅值响应曲线Fig.14 Heave motion amplitude response curves of calculated values and test values

图15 纵摇运动计算值与试验值的幅值响应曲线Fig.15 Pitch motion amplitude response curves of calculated values and test values

图16 船舯加速度计算值与试验值的幅值响应曲线Fig.16 Acceleration amplitude response curves of calculated values and test values

3 结论

本文采用计算流体力学对局部气垫双体船的艏艉气封进行了刚性化处理，对数值波浪水池、网格形式等影响计算精度的要素进行了研究分析，在此基础上，基于刚性气封初步对局部气垫双体船在规则波中的阻力、垂荡、纵摇和船舯加速度等运动情况进行仿真计算，并与试验数据进行对比，得到结论如下：

1）重叠网格相对于滑移网格有着更高的计算精度，在航速3.6 m/s工况时精度最高达到6.42%，在航速5.0 m/s工况时精度达到5.2%；

2）数值计算中波浪的阻力幅值和均值均大于试验值，在波长为7.0 m 时尤其明显，且从历时曲线的特征规律看，刚性气封模型的数值计算结果非线性特性更加突出；

3）从各运动参数的历时曲线中看出，垂荡运动和纵摇运动的计算历时曲线和试验值相比具有一定的相似性，在幅值和均值方面和试验值相比虽有误差，计算值和试验值没有和阻力一样展现出复杂的非线性特征；

4）对比各运动参数的幅值变化情况看出，各运动参数的计算值与试验值随波长的变化规律相近，且计算值小于试验值，而在λ/L为2.33处恰好相反。