蔡 威，张金凤

(天津大学 水利工程仿真与安全国家重点实验室，天津 300072)

在近些年来海洋资源的开发利用中，浮式结构由于其建造成本低、可预制、可拆卸、施工安装方便等诸多优点，越来越受到工程应用领域的重视。浮式结构在海洋环境中的水动力特性的相关研究方法主要有实验研究和数值模拟研究，其中实验研究具有成本高、时间长、操作复杂的特点，而数值模拟具有效率高的特点。浮式结构和流体相互作用的数值模拟研究问题的本质是一个具有自由表面的流固耦合问题，而且需要模拟物体的大幅度运动。这一系列特点与经典的物块入水抨击数值模拟相同，因此本文通过对物块入水抨击的数值模拟来验证数值模型的准确性。

对于两相流的自由表面追踪问题，研究者目前使用的主要有Level set方法[1]和VOF方法[2]。Gu等[3]使用Level set方法追踪自由表面，模拟固体的垂向和斜向入水过程；王平等[4]使用VOF方法对圆柱体的入水和下潜过程进行了模拟。对于流固耦合问题，数值模拟方面有很多的研究成果。胥飞等[5]使用SPH(Smoothed Particle Hydrodynamics)方法对二维船体入水砰击过程进行了模拟；Khayyer等[6]利用MPS(Moving Particle Semi-implicit)方法模拟了二维波浪对海岸结构物的冲击力；Johnson等[7]基于IBM(Immersed Boundary Method)方法模拟流动问题中的移动固壁边界；任安禄等[8]利用ALE(Arbitrary Lagrangian-Eulerian)方法来模拟圆柱绕流涡致振动；Rahman等[9]使用多孔体模型研究了半潜式浮式防波堤的运动；Shen Z[10]等使用重叠网格(Overset grids)方法研究了船只的自推进和机动；Quon等[11]使用重叠网格方法分析了波浪能发电装置的性能；Cheng等[12]使用重叠网格方法研究了浮式风机的气动性能。无网格的SPH方法和MPS方法都存在计算量较大的问题，而目前IBM方法不能较好地解决物体大幅度运动的问题，重叠网格方法可以模拟物体的大幅度运动，计算效率较高且对原有求解器的继承性好。

入水抨击是海洋工程和海军工程中的一个重要课题，通常做法是将物体假设为二维楔形体。Von Karman[13]最早采用附加质量近似代替流体作用来分析入水抨击问题，提出附加质量法计算水上飞机入水抨击荷载； Dobrovol[14]假设流体具有自相似性，忽略重力影响，使用自相似法，推导出了定速度入水时抨击压强分布的理论解。Mei[15]在Dobrovol的基础上使用边界有限元(Boundary Element Method，BEM)离散和求解方程，得到了更为一般情况下的解析解。

本文求解Navier-Stokes方程，在OpenFOAM开源程序的基础上进行二次开发， 利用VOF方法追踪自由表面，使用重叠网格技术，建立了流体(波浪)和浮体结构相互作用的二维数值模型，通过模拟楔形物块垂向和斜向入水过程，验证数值模型的准确性，并简单分析楔形物块入水规律。

1 数学模型

1.1 控制方程

本文通过求解有限体积法离散的RANS方程，采用VOF方法追踪自由表面，来模拟气液两相流的运动。流体基本控制方程RANS为

(1)

(2)

式(1)和式(2)分别为连续方程和动量方程，其中：u为速度；ρ为流体密度；v为运动粘滞系数；S为源项，其表达式为

S=Fσ+ρg=Cκα+ρg

(3)

(4)

式中：Fσ为作用在曲面微元上的表面张力之和；C为表面张力系数；κ为界面处的曲率；α为VOF方法中定义的体积分数，并且满足相方程

(5)

使用Menter[16]提出的SSTk-ω两方程紊流模型来闭合RANS方程，该模型综合了标准k-ω模型在近壁面区的计算优点和标准k-ε在远场的计算优点，其紊动能方程和紊流耗散率方程的表达如下

(6)

(7)

式中：k为紊动能；ω为紊动能耗散率。其他各项参数的表达式如下

(8)

(9)

(10)

(11)

(12)

模型中的系数α、β、σk和σω的计算方法为(用φ代表)

φ=F1φ1+F2φ2

(13)

式中：α1=5/9，α2=0.44；β1=3/40，β2=0.082 8；σk1=0.85，σk2=1；σω1=0.5，σω2=0.856；β*=0.09；a1=0.31。

1.2 重叠网格

图1 重叠网格示意图Fig.1 Schematic diagram of overset grids

楔形体的运动模拟采用重叠网格法。在这种方法中，计算区域网格根据所覆盖区域的特点，划分为多块具有重叠或者嵌套部分的子网格，如图1所示流场划分为背景网格(虚线)、固体划分为贴体网格(实线)。当物体运动时，贴体网格随之做无变形的刚体运动。数值计算在各个子网格中分别进行，同时在重叠的区域通过网格节点插值来实现计算信息的传递。

重叠网格方法的关键在于建立域连接信息(DCI， domain connectivity information)，用于子网格间的计算信息传递，其处理流程为：

(1)搜索洞单元，找出不参与计算的节点单元；

(2)为边缘(边界)节点从重叠区域的另外一套网格中寻找合适的插值贡献节点；

图2 重叠网格方法求解器的计算流程Fig.2 Flow chart of the overset grids solver

(3)根据边缘(边界)节点和插值贡献节点的位置关系，求解插值系数；

(4)优化重叠区域，减少计算量。

(14)

式中：φL为边缘(边界)节点上的某一流场计算信息，如压强、相分数等；αk是其对应的第k个插值贡献节点的插值系数；φk是其对应的第k个插值贡献节点的流场计算信息值。

如图1所示为重叠网格示意图，其中方形空心点为洞点，方形实心点为背景网格的边缘(边界)节点，圆形空心点为物体贴体网格的边缘(边界)节点。

对于定常边界问题，在求解物理问题之前，只需对DCI的建立进行一次处理，并将其存储在计算机内存中。对于移动边界问题，需要在流场计算期间重复这些步骤。图2给出了所开发的重叠网格方法求解器的整个计算流程。它与常规的CFD程序非常相似，区别在于包含了挖洞、DCI建立和插值。

2 数值模拟参数设置

表1 计算分组Tab.1 The groups of calculation

2.1 计算分组设置

图3-a所示为二维对称楔形体定速度入水的示意图，其中α和β分别为物块左侧和右侧与水平面所成夹角，θ为物块对称轴与竖直方向所成夹角，l为物块顶部长度，v和u分别为物块在Z方向和X方向的分速度。算例设置如表1所示。

2.2 网格划分

二维模型的背景网格为均分网格，尺寸为X=16 m、Z=16 m，X方向和Z方向的网格数均设为160，单个网格尺寸为0.1 m×0.1 m。楔形体部分的贴体网格设置为渐变网格，尺寸为X=6.2 m、Z=5.4 m，楔形体的顶部长度取2.4 m。为了较为准确地模拟楔形体底部的压力分布，楔形体底部使用snappyHexMesh进行五级渐变加密，最小网格尺寸为3.125 mm×3.125 mm。同时使用extrudeMesh只对Y方向的网格进行处理，将网格变为单层网格，可以显著减少网格数量，3个算例的计算总网格数分别为194 608、131 800和186 168。图4给出了算例一网格划分示意图。

3-a 二维对称楔形体3-b 模型示意图图3 二维对称楔形物块入水示意图Fig.3 Schematic diagram of two-dimensional symmetrical wedge block entering water

图4 算例一楔形体周围嵌套网格示意图Fig.4 Sketch of nested grid around wedge of case 1

2.3 边界条件

如图3-b所示，二维模型边界设置为：左右壁面采用固壁边界条件，法向梯度为零；前后面采用侧壁边界条件，法向始终不参与求解；楔形体采用物块边界条件，边界上的通量恒定为零；大气边界为自由表面边界条件；底边采用底床边界条件，边界上的速度恒定为零。

3 结果分析

二维水槽的水深设置为8 m，物块初始时刻下端点位于水面处且为水槽正中央的位置。

3.1 垂向入水

为了便于分析结果，避免液体密度和速度的影响，引入无量纲化的附加压强参数Cp，其表达式为

(15)

式中：ρ为水的密度，取1 000 kg/m3；v为物块入水的速度，取2 m/s；p为物块入水产生的附加压强，其表达式为

p=P总-p0

(16)

式中：P总为总压强；p0为静水时此处的压强。

图5给出了算例1和算例2的垂向入水时楔形体底部壁面上的附加压强分布。分别给出了0.05 s、0.06 s、0.07 s时的无量纲化压强分布曲线，结果表明不同时刻的数值模拟结果几乎一致，满足自相似性。作为对比，给出了Mei[15]的解析解和Dobrovol[14]的理论解，结果与Mei[15]的解析解吻合良好，与Dobrovol[14]的理论解存在明显的偏差。Mei[15]的解析解和Dobrovol[14]的理论解均基于势流理论和自相似性求解，但是Mei[15]的解析解进一步考虑了射流对壁面压强分布影响，更符合实际。

图6给出了算例1和算例2的垂向入水时楔形体底部壁面上的附加压强场分布。结合图5和图6，可以看出：(1)最大附加压强并不一定出现在楔形体下端点处；(2)高于静水面位置处的附加压强依然较大。算例1中的楔形体，最大附加压强位置在下端点处，而且附近的压强变化较小，压强大的区域分布较为集中。算例2中的楔形体，最大附加压强位置分布在距离下端点一段距离的底部壁面两侧，位置高于静水面，附近的压强变化较大。

5-a 算例15-b 算例2图5 垂向入水楔形体底部壁面的附加压强参数分布Fig.5 Distribution of additional pressure parameter on the bottom wall of the vertical water entry wedge

6-a 算例16-b 算例2图6 垂向入水楔形体在0.05 s时的附加压强场分布Fig.6 Distribution of additional pressure field around the vertical water entry wedge at 0.05 s

7-a u=0.6 m/s,v=-2 m/s7-b u=-0.6 m/s,v=-2 m/s图7 斜向入水楔形体底部壁面附加压强参数分布Fig.7 Distribution of additional pressure parameter on the bottom wall of the oblique water entry wedge

3.2 斜向入水

实际情况中，物块入水会有一定的斜向速度，斜向入水速度会对底部壁面上的压强分布产生影响。Xu[17]在Mei[15]的基础上，进一步考虑了入水时的射流对压强分布的影响，使用边界元分析了楔形体斜向入水问题。本文模拟了水平向速度分别为0.6 m/s 和 -0.6 m/s，垂直速度为-2 m/s，即水平速度和垂直速度绝对值的比值u/|v|为0.3和-0.3时，算例3的底部壁面在0.05 s时刻的附加压强分布，结果与Mei[15]的BEM结果吻合良好(如图7所示)。同时给出了水平与垂向速度绝对值之比u/|v|为0.1、0.5、-0.1和-0.5的情况，其结果如图8所示。

对比分析不同水平与垂向速度之比的情况，可以看出：(1)附加压强的极值出现的相对位置几乎不变；(2)附加压强的最大值随着水平与垂向速度之比的变化而出现明显的变化。对于算例3的楔形体，附加压强的最大值对应的相对位置X/vt≈-4.14。由于算例3中锲形体入水角度α=20°和β=40°，楔形体入水时左侧更贴近水面，因此入水时，左侧水面受到挤压更大，压强极大值会出现在左侧。当水平速度为0时，无量纲化压强参数极大值为16.87，位于楔形体底部壁面左侧。当楔形体斜向左下方入水时，由于左侧的挤压效应增强，随着水平向的速度增大，附加压强的最大值增大。当水平与垂向速度绝对值之比为-0.5时，无量纲化附加压强参数的极大值为21.46。当楔形体斜向右下方入水时，由于左侧的挤压效应减弱，随着水平向的速度增大，压强的最大值减小。当水平与垂向速度绝对值之比为0.5时，无量纲化压强参数极大值为13.10。

8-a 楔形体斜向左下方入水8-b 楔形体斜向右下方入水图8 不同u/|v|时斜向入水楔形体底部壁面压强分布对比Fig.8 Distribution of additional pressure parameter on the bottom wall of the oblique water entry wedge in different u/|v|

4 结论

本文采用开源计算流体力学软件OpenFOAM的两相流求解器和重叠网格方法，建立了可以模拟物体大幅度位移运动的流固耦合模型。通过对经典的物块入水问题的模拟，分析和验证了楔形体物块垂向和斜向入水状况下，物块底部壁面上的压强及其极值位置的分布。数值模拟结果与理论分析或已有数值模拟结果吻合良好，同时阐明了楔形体垂向和斜向入水的压强分布规律，认为水平速度的增大会改变附加压强的极值大小，但是不会改变其出现的相对位置。