0 引　言

穿浪双体船最早于20世纪80年代初提出，其航行时船首频繁出水入水，由于船体与波浪之间剧烈的相对运动，船首处容易发生砰击现象。在入水瞬间，结构的运动引起周围流体介质运动，流体介质同时对结构物施加反作用力，此时结构物将遭受巨大的流体砰击载荷，这会对结构体的结构强度造成严重威胁。本文将穿浪双体船船首剖面简化为2个完全相同的二维楔形体。

Zhao和Ma等^[1-2]研究了无重力影响时二维楔形体压力分布的峰值位置，表明底升角大约大于40°时峰值位置由射流根部转移到楔形体尖端。Chen^[3]建立的二维直升机V形浮标的数值模拟结果也验证了Zhao和Ma等^[2]的结论。陈震^[4]表明二维楔形体入水压力峰值的传播速度随着入水角度的增大而减小。张健^[5]、王玲等^[6]、Chen^[7]和陈光茂等^[8]分别对二维楔形体的入水砰击问题进行仿真分析，研究压力峰值与底升角、入水角度和入水速度的关系。Ma^[9]模拟了二维楔体垂直和斜入的入水和水动力冲击过程。Wu等^[10]表明有限高度非对称二维楔形体斜入水速度增大时，一段较长时间内的重力效应可以忽略。Sun等^[11]首次分析了二维楔形体在波浪中的入水过程，并研究了重力对楔形体入水过程和物面压力分布的影响。陈翔等^[12]表明了楔形体入水底升角越小流体沿物面爬升速度越快，射流越明显。Zhang等^[13]表明对称楔形体入水过程中两侧受到的压力也是对称的。

本文采用OpenFOAM开源软件，建立具有不同底升角单楔形体和双楔形体在静水中入水的CFD数值模拟模型，采用有限体积法（FVM）离散平均雷诺数方程和连续性方程组成的控制方程组，采用流体体积法（VOF）追踪自由表面。将基于粘流理论的CFD数值结果与前人基于势流理论的数值结果进行对比分析，研究粘流方法与势流方法在处理入水问题过程中的差别，并重点分析重力和粘性在此过程中所起的作用。

1 数值过程

对于不可压缩粘性流体，可采用连续性方程和雷诺平均方程来建立流场的控制方程组。连续性方程形式为：

$ \nabla \cdot u = 0 \text{，} $

(1)

式中，u为速度矢量。

雷诺平均方程可以表达为：

$ \dfrac{{\partial {\boldsymbol{u}}}}{{\partial t}} + \left( {{\boldsymbol{u}} \cdot \nabla } \right){\boldsymbol{u}} = - \frac{{\nabla p}}{\rho } + \nu {\nabla ^2}{\boldsymbol{u}} + {\boldsymbol{F}} \text{。}$

(2)

式中：ν为流体运动粘性系数；F为质量力；ρ是流体密度；p为压力。

采用VOF方法逼近空气与水的交界面，在模型中不同的流体状态共用着一套控制方程组，通过引进相体积分数这一变量，实现对每一个计算单元相界面的追踪。对于波物相互作用流场，存在空气和水2种流体流动，用 $ {\rho _{{\text{air}}}} $ 表示空气密度，用 $ {\rho _{{\text{water}}}} $ 表示水密度，流场内流体密度可表达为：

$ \rho = \alpha {\rho _{{\text{water}}}} + \left( {1 - \alpha } \right){\rho _{{\text{air}}}} \text{，}$

(3)

$ \left\{\begin{array}{l}\alpha =\text{0}\text{，}空气 \text{，} \\ \text{0} \lt \alpha \lt 1\text{，}自由表面 \text{，} \\ \alpha =\text{1}\text{，}水 \text{。} \end{array} \right.$

(4)

使用SST k-ω两方程湍流模型来闭合RANS方程，SST k-ω湍流扩散方程为：

$ \frac{\partial }{{\partial t}}\left( {\rho k} \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho k\overline {{u_i}} } \right) = \frac{\partial }{{\partial {x_j}}}\left( {{\Gamma _k}\frac{{\partial k}}{{\partial {x_j}}}} \right) + {G_k} - {Y_k} + {S_k} \text{，} $

(5)

$ \frac{\partial }{{\partial t}}\left( {\rho \omega } \right) + \frac{\partial }{{\partial {x_j}}}\left( {\rho \omega \overline {{u_j}} } \right) = \frac{\partial }{{\partial {x_j}}}\left( {{\Gamma _\omega }\frac{{\partial \omega }}{{\partial {x_j}}}} \right) + {G_\omega } - {Y_\omega } + {D_\omega } + {S_\omega } \text{。}$

(6)

式中： $ k $ 为湍流动能； $ \overline {{u_i}} $ 和 $ \overline {{u_j}} $ 为速度的水平分量和垂直分量； $ \omega $ 为湍流特征频率； $ {\Gamma _k} $ 和 $ {\Gamma _\omega } $ 分别为 $ k $ 和 $ \omega $ 的有效扩散项； $ {G_k} $ 和 $ {G_\omega } $ 分别为 $ k $ 和 $ \omega $ 的生成项； $ {Y_k} $ 和 $ {Y_\omega } $ 分别为 $ k $ 和 $ \omega $ 的耗散项； $ {D_\omega } $ 为横向扩散项； $ {S_k} $ 和 $ {S_\omega } $ 为自定义源项。

2 数值结果分析 2.1 物理模型和网格划分

如图1所示，笛卡尔坐标系 $ x $ 轴方向水平向右， $ y $ 轴方向竖直向上，初始时刻楔形体尖端位于原点且静置于静水面上。流体以速度 $ W $ 自下而上流动，等价于楔形体以速度 $ W $ 沿 $ y $ 轴负方向运动。楔形体固定高度 $ h $ 为5 m，底升角 $ \gamma $ 分别采用 $ {\text{45}}^\circ $ 和 $ {\text{30}}^\circ $ ，边界条件见图1。

使用blockMesh工具绘制背景网格，计算域尺度为宽度6 m，厚度0.5 m，高度5 m，其中空气域为2 m，水深为3 m，如图2所示。为了较为准确地模拟楔形体的压力分布，使用snappyHexMesh网格生成器分裂加密楔形体周边网格，最小网格尺度为2.6 mm×2 mm，如图3所示。

2.2 数值结果分析

为了便于分析问题，在下面数值结果中，如无特殊说明，本文将采用垂直速度 $ W $ 、重力加速度 $ g $ 和水密度 $ \rho $ 对数值结果进行无量纲化，将长度、时间 $ t $ 、速度、压力 $ P $ 的量纲分别用W²/g，W/g， $ W $ 和 $ \rho {W^2} $ 表示。

2.2.1 45°底升角单楔形体入水

$ \gamma = {\text{4}}5^\circ $ 单楔形体不同工况下的结果如图4～图6所示。其中图4（d）和图5（d）给出了无量纲自由液面形状，图6的相似解引自Sun等^[11]的基于势流理论的45°底升角无限高度楔形体匀速入水时的压力分布曲线。

由图4和图5可知，当入水时间较短时，自由液面沿物面快速爬升。随着入水时间的增加，射流尖端到达楔形体顶部，射流沿楔形体顶部飞出，形成了射流分离现象。入水时间相同，重力的作用效果也基本相同，而图5比图4射流更高的直接原因是其对应的楔形体入水速度更高，被楔形体扰动的射流获得了更高的分离速度。不考虑重力和射流分离影响时，楔形体的入水问题是自相似的，而图5中 $ t = 0.1s $ 时刻的自由面无量纲曲线与其余2个时刻产生偏差，图5对应的楔形体入水距离更大，重力影响变明显，射流无法无限期维持。

图6（a）中 t= 0.01 s 的数值结果与其余3条曲线有明显差别。这是由于采用VOF方法确定的气液交界面是一条密度逐步过渡的细带状区域，这将导致采用CFD方法计算的压力在入水瞬间计算的不够精准。速度提升了1倍后，相同时刻楔形体入水距离已经达到一定距离，楔形体表面压力不再受自由液面影响，见图6（b）。但同时在0.1 s，物面压力又急剧降低，这是由于相似解没有考虑重力，而基于粘流理论的CFD结果考虑了重力，随时间增加重力的作用越来越明显，导致物面压力急剧下降，并且远离相似解。此外，所有CFD压力数值结果都是小于相似解的，前面描述的重力只是其中的一个次要原因，另外一个重要原因是砰击发生时物面压力梯度较大，CFD方法由于网格尺度限制在处理大压力梯度问题时存在一定困难。

2.2.2 30°底升角单楔形体入水

$ \gamma = {\text{30}}^\circ $ 单楔形体不同工况下的结构如图7～图9所示。其中图7（d）和图8（d）给出了无量纲自由液面形状，图9的相似解引自许国冬^[14]的基于势流理论的30°底升角无限高度楔形体匀速入水时的压力分布曲线。

相同入水速度相同时刻， $ \gamma = 30^\circ $ 楔形体入水瞬间引起的射流均比 $ \gamma = {\text{45}}^\circ $ 楔形体的更加细长，符合许国冬^[14]和Sun^[11]等的研究成果。当楔形体底升角越小，流体质点流动速度更快，射流更长。图8（d）在 t = 0.1 s 时刻分离射流达到较大的高度和体积，由于自由射流获得的惯性向上速度来源于楔形体，当它脱离楔形体以后所受的作用力只有重力，使得自由射流的惯性速度变小。越接近自由射流尖端，脱体时间越长，受到的重力影响就越明显，因此射流尖端附近的流体质点速度首先降低，出现向右下方弯曲的趋势。

除上一个算例对压力曲线的分析原因，影响压力另外一个重要原因是射流分离和重力，见图9（b）的0.1 s压力结果，附体射流层的0压力线在曲线中消失，且压力值降低。

2.2.3 双楔形体以给定速度入水的数值结果分析

如图10所示，双楔形体A和B两个顶点的距离设置为0.5 m，其余设置与图1一致。

$ \gamma = {\text{4}}5^\circ $ 双楔形体不同工况下的结果如图11～图14所示。其中图13和图14分别引入图6与之对应的单楔形体数值结果。

由图11可知，自由面沿物面快速爬升，形成粘附于物面的附体射流，最终脱离物体顶部，形成一小段自由射流，在双楔形体中间，自由液面连同自由射流形成一个开口气腔。图12中双楔形体入水速度提高一倍，附体射流较快地演变为自由射流， $ t = 0.1\;{\rm{s }}$ 时双楔形体之间的液面形式为一个凸起的水柱。这是由于入水速度足够大时，楔形体扰动的自由液面范围扩大，使得双楔形体内的两股液面融合，自由液面以水柱形式抬升。

由图13和图14可知，压力曲线整体趋势和峰值出现的相对位置几乎不变，双楔形体压力峰值都大于单楔形体压力峰值。图13（d） t = 0.1 s 的数值结果与其余不同，双楔形体压力显著高于单楔形体压力。这是因为对于单楔形体入水，一部分流体受到扰动以后以射流形式沿着物面向上运动，一部分流体沿着自由面向远方流动，而对于双楔形体，向远方流动的流体被前方楔形体阻挡，被阻断的流体只能以更高的速度向上运动来满足动量守恒定律，因此物面压力显著上升。可以预见，在砰击发生的短暂时间内，入水速度越大，入水时间越长，双楔形体之间的相互压力影响便越显著。当时间足够大时，重力的影响变得明显，这一结论不一定适用，但是长时间的流体与结构相互作用问题并不是本课题的研究内容。14（d）的差别更大是由于速度越高，流体扰动速度越大，流体更快达到相邻楔形体，相邻楔形体压力也升高的更快。

$ \gamma = 30^\circ $ 双楔形体不同工况下的结果如图15～图18所示。其中图17和图18分别引入图9与之对应的单楔形体数值结果进行对比。

不同于图11的气腔开口，图15的气腔基本封闭了，原因是小底升角导致射流速度更高，射流长度变得更长，使得气腔更容易封闭，所以图16（c）中的双楔形体中间的拱形射流也具有更高的高度。与上一个算例类似，在图17和图18中，由于双楔形体之间的相互影响需要时间，时间越长，单楔形体和双楔形体压力差别越大。与上一个算例不同的是，本算例楔形体底升角更小，受到扰动的流体运动得更快，因此本算例单楔形体和双楔形体压力之间的距离较前一个算例要大。由此可以得出另外一个结论，在楔形体间距不变的条件下，小底升角会使楔形体之间的影响时间提前，影响幅度扩大。

3 结　语

本文采用基于粘流理论的CFD方法，模拟重力作用下不同底升角单楔形体和双楔形体以给定速度在静水中的入水过程，得到如下结论：

1）当入水时间较短时，不同时刻的无量纲自由液面曲线基本重合，随着时间的增加，射流从物体顶部分离，自由射流在重力作用下发生弯曲向下的现象。

2）当入水时间较短时，排除初始阶段的不准确压力，不同时刻的物面无量纲压力曲线是基本重合的，并且与基于势流理论的自相似解非常接近，随着时间的增加，重力的作用变得显著，压力在重力作用下降低。

3）在相同入水速度和相同底升角条件下，双楔形体物面压力大于单楔形体物面压力，并且在砰击发生的短暂时间范围内，刚开始影响不大，从某一时刻开始，受一个楔形体扰动的流体到达相邻楔形体以后，相邻楔形体压力差值瞬间变得显著。此外，若提升速度，被楔形体扰动的流体速度加快，相邻楔形体受到的压力影响更大。

4）当降低楔形体底升角时，相当于增加楔形体周围被扰动的流体速度，因此相邻楔形体受到的压力也提升的更快。