0 引　言

自航绞吸船为布置型船舶，根据布置和功能需求，在船首布置钢桩定位系统，或者挖掘系统，前者称为“定位艏型”，后者称为“挖掘艏型”^[1]。对于“定位艏型”自航绞吸船，为了满足钢桩台车布置和行走空间的需求，需要在船体上开槽，从甲板贯通，至与船底海水相连，形成月池结构。航行时，其内部流体流动复杂，会产生水流运动和波浪破碎等强非线性现象，使航行阻力显著增加。

针对月池内部水流运动特征，Fukuda^[2]通过试验，研究了船舶在固定和自由运动状态下，矩形和圆形开口的月池内部水流运动特征，分析了其对于船舶阻力和运动的影响，提出了月池固有频率的计算方法，试验值与理论值比较相符，对于降低月池内部运动提出了有效的解决方案；Molin^[3]提出了船舶月池内部水流的2种运动形式，一是晃荡模式（Sloshing mode），即水流在月池内前后运动；二是活塞模式（Piston mode），即水流在月池内做升沉运动。采用线性势流理论，提出了在这2种模式下计算二维和三维矩形月池的固有频率的理论公式；Newman^[4]在Molin研究工作的基础上，针对具有1个台阶形式的矩形月池，给出了其晃荡模式下各阶数固有频率的解析解；Senthuran等^[5]基于线性势流理论，采用边界元法（BEM）和CFD方法，对具有台阶形式的二维月池在晃荡模式和活塞模式进行了求解并与模型试验结果进行了对比。其认为，在预报2种模式下月池内部自由液面的最大运动响应方面，CFD具有更好的精确度；鲜于晨松等^[6]采用势流理论和CFD相结合的方法，研究了钻井船月池内流体运动，通过模型试验，验证了CFD计算的准确性。

在月池引起的船舶航行阻力增加方面，Sivabalan等^[7]采用势流理论对具有矩形月池和无月池的船舶在低航速下的阻力进行了计算和对比，结果表明，在相同航速下，前者的阻力大于后者，且随着航速增加，二者的差距也更加明显；Yang等^[8]对具有不同形状的月池活塞和晃荡2种模式下的运动响应开展了模型试验，验证了Molin理论公式对于月池固有频率计算有较好的精确度，提出了优化月池形状的有效减阻方案；SON等^[9]对具有不同台阶形式月池的钻井船进行了数值计算和船模试验，结果表明，相对于常规矩形月池，月池形状优化的船舶减阻效果达到10%以上；Guo等^[10]采用数值计算和模型试验的方法，对1艘具有台阶形式矩形月池的钻井船开展研究，结果表明，CFD方法具有较高的精确度，也能够很好地模拟月池内部水流和自由液面的物理特性和运动情况；李文娟等^[11]针对大洋勘探船，采用CFD方法分析了不同月池方案对于船舶总阻力及月池增阻情况，研究了摩擦阻力和剩余阻力分量变化的规律；许建龙等^[12]研究了环形月池对于海洋核能平台拖航阻力的影响，通过CFD计算，认为月池的存在，使得平台的拖曳阻力增加5%～15%；仉永超等^[13]采用CFD方法，针对一艘潜水支持船进行了分析，结果表明，月池位置变化对于其阻力变化影响不大，月池的设计更多地考虑船舱内设备的布置需要。

然而，以上研究主要是以具有规则形月池的钻井船为主要研究对象，对于超大型自航绞吸船的相关研究较少。这2种船型因布置和功能需求不同，因此，在水动力性能的侧重点上存在较大差异；另外，由于钢桩台车的特殊性，即作业工况放下使用，航行工况收起固定的特点，可以考虑在满足其布置需求的前提下，对于开槽进行设计和优化，或者采用合理的封盖，来进一步减小航行阻力。本文基于船模试验和CFD方法，对一艘“定位艏型”超大型自航绞吸船（下简称“目标船”）进行数值计算和对比分析。首先，通过船模试验，验证CFD方法的可靠性；通过对比不同开槽设计优化方案，分析航行阻力成分变化趋势，了解封盖的减阻机理；进一步地，为了更好地理解开槽内部水流运动特征和物理现象，采用时域和频域方法来追踪和分析自由液面的变化，确定目标船开槽内部水流的运动模式，分析了其对于航行阻力的影响。得到了减阻效果良好、实际操作可行的钢桩开槽封盖方案，为提高实船的水动力性能提供技术支撑。

1 基本理论 1.1 开槽内水流运动固有频率

开槽内水流运动固有频率为：

$ f = \frac{1}{T} = \frac{1}{{\dfrac{{2\text{π} }}{\omega }}} = \dfrac{\omega }{{2\text{π} }} 。$

(1)

对于矩形月池，当开槽内水流运动为活塞模式时：

$ {\omega _n} = \sqrt {\frac{g}{{d + d'}}} 。$

(2)

式中：d为月池内吃水，即船舶吃水，m；$ d' $为附连水引起的吃水增加值，m。Fukuda^[2]给出其经验值为：

$ d' = 0.41\sqrt S 。$

(3)

式中：$ S $为月池水线面面积，m²。

而Molin^[3]给出了附连水引起的吃水增加值解析解：

$ \begin{split} d' = &\frac{b}{\text{π} }\left\{ {{\sinh }^{ - 1}}\left( {\frac{l}{b}} \right) + \frac{l}{b}{{\sinh }^{ - 1}}\left( {\frac{b}{l}} \right) + \frac{1}{3}\left( {\frac{b}{l} + \frac{{{l^2}}}{{{b^2}}}} \right) - \right.\\ &\left.\frac{1}{3}\left( {1 + \frac{{{l^2}}}{{{b^2}}}} \right)\sqrt {\frac{{{b^2}}}{{{l^2}}} + 1} \right\} 。\end{split}$

(4)

式中：b为月池宽度，m；l为月池长度，m。

当开槽内水流运动为晃动模式时，Fukuda^[2]给出经验值：

$ {\omega _L} = \sqrt {\frac{{n\text{π} g}}{l}}。$

(5)

式中，n为阶数，特别地，当n = 1时，即为1阶晃荡固有频率。

而Molin计算值为：

$ {\omega _L} \cong \sqrt {\frac{{n\text{π} g}}{l}\left(\frac{{1 + {J_{n0}}\tanh (n\text{π} d/l)}}{{{J_{n0}} + \tanh (n\text{π} d/l)}}\right)}。$

(6)

式中：$ {J_{n0}} $为关于$ n\text{π}b/l $的积分函数，其结果不大于1。Molin已绘制成图谱，可根据$ n\text{π} b/l $的计算值进行查询得到^[3]。

1.2 控制方程和湍流模型

对于不可压缩流体，连续性方程（质量守恒方程）如下：

$ \nabla \cdot u = \frac{{\partial {u_x}}}{{\partial x}} + \frac{{\partial {u_y}}}{{\partial y}} + \frac{{\partial {u_z}}}{{\partial z}} = 0 。$

(7)

式中：$ \nabla $为哈密顿算子；$ {u_x} $、$ {u_y} $、$ {u_z} $分别为x、y、z方向上的流体的平均速度。

流体的动量守恒方程：

$ \frac{{\partial {u_i}}}{{\partial t}} + {u_j}\frac{{\partial {u_i}}}{{\partial {x_i}}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial {x_i}}} + \upsilon \frac{\partial }{{\partial {x_j}}}\left(\frac{{\partial {u_i}}}{{\partial {x_i}}}\right) + {f_i} 。$

(8)

式中：ρ为流体密度；t为时间；p为微元流体压力；$ \upsilon $为流体运动黏度；$ {f_i} $为微元流体单位质量力；下标i、j分别为笛卡尔坐标系横坐标与纵坐标的方向。

由于N-S方程中增加了雷诺应力项，使得控制方程未知量个数大于方程个数。为解决该问题，引入SST k-ω湍流模型对N-S方程予以封闭，即可求解^[14]。

1.3 虚拟水池技术

采取面重构网格、切割单元体网格和边界层网格的网格划分模型。在首尾处、湍流层、船体、首尾开槽以及呆木周围采取网格加密技术，以便更加精确地捕捉流场细节。

采取虚拟水池的方法进行数值模拟计算。由于船体对称，简化成一半进行处理。虚拟水池长度为4倍船长，进口距离船首为1倍船长，出口距离船尾为2倍船长，半宽为2倍船长。整个计算域以水面为界，上部为空气，下部为水，计算域高度方向为水面以上1倍船长，水面以下2倍船长。设置速度入口和压力出口边界，中纵剖面和计算域侧壁面设为对称面边界，船体表面为不可滑移壁面，近壁面处设置边界层条件。边界条件设置见图1。

边界层划分对于准确预报船舶阻力非常重要，尤其是第1层网格的厚度，引入无因次参数$ y^{+ } $来表示。计算公式为：

$ y{\text{ + }} = 0.172\left( {\Delta y/{L_{PP}}} \right)R{e^{0.9}}。$

(9)

式中：$ \Delta y $为第1层边界层厚度；$ {L_{PP}} $为船长；$ Re $为雷诺数。对于船舶船模尺度的数值模拟而言，$ y^+ $的取值要求为30～300^[15]。

2 研究对象

目标船船体上需设置钢桩台车开槽，并为钢桩台车配置包括沿纵向行走的水平轨道和多条侧向限位的竖直轨道。钢桩台车开槽的长度和宽度主要取决于钢桩直径、钢桩台车型式和钢桩台车行程。钢桩台车开槽长度大致等于钢桩台车长度和钢桩台车行程之和加上一定裕量。钢桩台车开槽宽度大致等于钢桩台车主体结构宽度与两侧的水平限位滑块结构宽度之和，船底处的宽度还可合理设计减小。一般而言，钢桩台车的横向开槽宽度为上宽下窄，这样便于安装和移去钢桩台车。

针对目标船，在满足钢桩台车布置和行走的需求的前提下，开展开槽设计优化。为了对比，以完整船体作为基准，即船首无开槽方案，记为F₀；按照一般设计，采用矩形开槽，称“常规开槽”方案，记为F₁；钢桩台车在船体开槽中移动时，船底处的最小开槽尺寸需满足钢桩移动所需空间需求，此时在船底增加部分固定浮箱结构，称“钢桩开槽”方案，记为F₂；对优化的目标船钢桩台车船底方案，拟配置移动式的封盖结构，目标船调遣航行前安装，把钢桩台车航海固定时钢桩下抱箍前的船底开槽封盖，见图2，称“钢桩开槽封盖”方案，记为F₃。船体及开槽模型见图3。

目标船主要尺度要素见表1，开槽尺寸见表2。模型缩尺比为25。

在开槽吃水高度处设置3个波高测量点，见图4。

3 阻力计算验证

针对方案F₁，开展阻力试验。在上海交通大学船海工程试验中心船模拖曳水池进行。水池全长为300 m，宽为16 m，深为7.5 m，并配有最大速度为10 m/s的高速拖车。船体模型由松木制作而成，并对其表面进行了打磨及喷涂油漆处理，保证船体的光滑度。模型艏部加装了激流丝以模拟实际湍流状态。实验前，在水平台上调节压铁在船体不同位置处的配重，以保证船体首尾吃水、排水量、重心高度、转动惯量等参数与试验要求保持一致。在吃水为0.28 m工况下开展试验。试验为淡水，水温为10.1℃，水密度ρ为999.59 kg/m³，运动粘性系数ν为1.3027$ \times $10⁻⁶ m²/s。

船体及虚拟水池网格划分及边界条件设置见1.3节相关内容。图5为船体网格划分，图6为船体$ y^{ + } $分布。数值计算采用隐式非稳定求解，应用SST k-ω湍流模型，采用流体体积函数法（Volume of Fluid Method，VOF）来捕捉自由液面，时间步长满足ITTC推荐值^[16]，取为0.02 s。

采用ITTC推荐的方法，对方案F₁进行数值模拟误差和不确定度分析^[17]。基于网格收敛性研究，以网格细化率$ {r_G} $为$ \sqrt 2 $生成3套不同计算网格，3套网格从细到疏，分别为Grid-1、Grid-2、Grid-3。网格数如表3所示。

对总阻力数值计算开展不确定度分析。计算结果如表4所示。总阻力检验如表5所示。表中，$ {R_G} $为网格收敛率，$ {P_G} $为准确度阶数精度，$ {C_G} $为修正因子，$ {U_G} $为计算数值不确定度，$ \delta _G^* $为网格误差，$ {U_{GC}} $为修正后数值不确定度，$ {S_C} $为修正解。有效性验证如表6所示。表中，$ |E| $为比较误差绝对值，$ |{E_C}| $为修正比较误差绝对值，$ {U_V} $为比较误差不确定度，$ {U_{VC}} $为修正后比较误差不确定度。

由表5和表6可知，网格收敛率满足单调收敛条件0<$ {R_G} $<1，故3套计算网格为单调收敛，满足网格收敛性；$ \delta _G^* $和$ {U_{GC}} $均小于2%$ {S_C} $，表明总阻力数值模拟满足数值误差和不确定度要求，CFD方法得到检验；由$ |E| $<$ {U_V} $，$ |{E_C}| $<$ {U_{VC}} $可知，表明总阻力数值模拟的有效性也得到了验证。

综合考虑计算精度和计算时长，选取Grid-2作为数值计算网格划分方案，开展后续研究。

4 结果与分析

在吃水为0.28 m工况下，选取3个航速点，即船模航速V_m为0.823、1.029和1.235 m/s，针对不同船体模型开展数值计算。

4.1 阻力特性

方案F₁和F₀总阻力比较见图7。各航速下，相对于完整船体F₀，开槽使得模型F₁总阻力显著增加，其中，设计航速处增加了70.21%。可知，有必要对开槽进行进一步优化设计，以降低航行阻力。

按照傅汝德观点^[18]，即：

$ {R_\text{tm}} = {R_\text{fm}} + {R_\text{rm}}。$

(10)

式中：R_tm为模型总阻力；R_fm为模型摩擦阻力；R_rm为模型剩余阻力，单位均为N。

为了比较不同开槽方案的航行工况下的阻力特性，以阻力试验模型，即方案F₁为基准，对3个方案进行阻力成分分析，见图8。

各航速下，相对于F₁，F₂和F₃的总阻力均有所减小，F₂平均减小13.20%，F₃平均减小25.04%；摩擦阻力均有所增加，F₂平均增加2.24%，F₃平均增加2.46%；剩余阻力均有所减小，F₂平均减小15.44%，F₃平均减小27.51%。

由此可知，开槽封盖的主要作用是降低了船体的剩余阻力。虽然摩擦阻力有所增加，但是剩余阻力的减小量大于摩擦阻力的增加量，因此，总阻力随之减小，且随着航速增加，台车船底开槽封盖的减阻效果越明显。

4.2 开槽内水流运动

在船模航速为1.235 m/s工况下，船舶航速稳定后，对这3个方案钢桩台车开槽内的水流运动进行分析，见图9。对于F₁，在周期为0T时，在开槽底部靠船首形成了第1个小的漩涡，在1/4T时，该漩涡向后移动并逐渐增大，同时，在靠近开槽艉部也逐渐形成了第2个漩涡，在1/2T时，这2个漩涡逐渐增大，随着第1个漩涡向后移动，在3/4T时，这2个漩涡融合成1个大的漩涡，并随水流流向船后。该过程具有周期性；对于F₂，其过程与F₁类似，也具有周期性。差别在于：第1个漩涡和第2个漩涡生成时与水平面存在一定的角度；F₂漩涡融合流向船后时，受到了开槽底部靠船尾的封板的阻挡，有一定的减缓作用；F₁开槽内部存在水气混合区域，分界线较为明显，而F₂开槽内部以水流运动为主。

对3个方案开展波高测量，得到了无因次波高随时间变化的时域曲线，通过傅里叶变换：

$ F\left( \omega \right) = \mathcal{F}[f(t)] = \int\nolimits_{ - \infty }^\infty {f(t)} {e^{ - i\omega t}}{\mathrm{d}}t。$

(11)

式中：f(t)为时域函数；F(ω)为频域函数。

得到振幅随频率变化的频域曲线（见图10）。频率分析结果见表7。

对于F₁，从无因次波高时域图中可以看到，在3个观测点处，波高变化较为复杂，3个位置处的固有频率均为0.750 Hz。数值计算得到固有频率，更接近于Fukuda的计算结果。而采用Molin计算方法，比上述结果均为偏大。主要原因是，Fukuda和目标船均考虑了船舶航速的影响，而Molin假定船舶为零航速。

对于F₂，在W₂观测点处，波高变化相对复杂，说明此处的水流运动存在较强的非线性，这可以从图9中得到验证。而在W₁和W₃观测点处，波高呈正余弦变化；从频域图中可以了解到，W₁和W₃处水流运动的固有频率为0.450 Hz，W₂位置处固有频率接近0.875 Hz。

对于F₃，在3个观测点处，波高几乎均维持不变；从频域图中可以了解到，波高变化不存在固有频率。

5 结　语

1）通过阻力试验，对船模总阻力进行验证，结果表明，CFD方法得到了检验和验证（V&V），总阻力计算误差不超过5.0%，证明数值计算方法的可靠性。

2）相对于完整船体，常规开槽方案的总阻力显著增加，船舶在设计航速航行时，最大增阻达70.21%。封盖的主要作用是降低了船体的剩余阻力。随着航速增加，封盖的减阻效果越明显，设计航速处可达28.55%。

3）在开槽内水流运动模式上，常规开槽方案以晃荡为主，钢桩开槽方案在靠船首的一侧主要为活塞，在靠船尾的一侧主要为晃荡的混合模式；钢桩开槽封盖方案开槽内部水流高度保持不变，几乎不存在水流运动。

本文基于常规开槽船型阻力试验，进行了基于总阻力的数值计算验证，并未对开槽内部的水流运动开展测量试验，这是需要后续进一步考虑和深入研究的。同时，船舶采用了固定状态，需进一步考虑在运动状态下，其阻力和开槽内水流运动特性；另外，对于开槽优化设计的减阻效果，也需要通过船模试验进行验证，予以闭环。从开槽封盖的实用性和便利性上，也需要在目前方案的基础上，进一步优化和完善，如封盖的形式和重量、航海固定方法等。