活塞式造波池中箱式浮体运动的数值研究

廖 斌，陈善群

(安徽工程大学 建筑工程学院，安徽 芜湖 241000)

浮体一般是指漂浮于液体表面的物体，诸如浮标、船只、浮桥以及大型网箱等，在海洋及河流等大型水体中均极为常见。浮体在水体中易受波浪的驱动而产生一定的摇晃，当波浪作用较强时甚至会出现浮体的倾覆以及相互拍击等现象，直接影响浮体的稳定性，进一步危及浮体上设备的正常工作以及人员的安全[1]。同时，波浪中浮体运动的研究涉及到流固耦合、自由面精细追踪等一系列学术问题。因此，研究波浪中的浮体运动具有较为重要的现实意义和学术价值。

国内外对于波浪中浮体运动的研究由来已久，使用的手段目前主要集中在实验和数值研究两方面。Zhao[2]等对三种类型波浪(规则波、聚焦波和规则波与聚焦波的组合)驱动下倒T型箱型浮体运动进行了实验观测，获取了浮体在3种类型波浪作用下的起伏位移以及摇荡角度等相关数据。He[3]等使用电荷耦合器件摄像机捕获浮体的实时运动，对6种波高和6种周期条件下波浪驱动中的箱型浮体运动的非线性行为进行了实验研究。Wang[4]等对规则和不规则波浪驱动下沙漏型浮体运动进行了实验研究，获得了波浪作用下浮体运动的摇荡周期、粘性阻尼和运动响应等实验数据。Hashiura[5]等和Sun[6]等对长宽比较大的桥式浮体在规则波浪中的运动响应分别进行了实验研究，重点对浮体所受的波浪力、弯矩等数据进行了深入分析。Ji[7]等对并列水平圆柱和矩形单浮桥浮式防波堤的水动力性能进行了实验研究，对浮体的系泊索力和三自由度运动(摇摆，升沉和滚动)进行了系统评估。Yang[8]等设计了一种新型的压载水浮式防波堤，并对该浮体在二维规则波驱动下的流体力学效率和运动响应进行了深入分析。显然，实验研究可以获取浮体运动的大量直观数据，但所需的测试场地以及采集设备的经济成本较高，能实现的波浪类型和浮体种类也相对有限。基于此，部分学者将目光聚焦于波浪驱动下浮体运动的数值研究，其中采用的数值方法主要分为有网格方法和无网格方法两大类。采用有网格方法，Wu[9]等基于二维三叉树方法生成运动网格，数值研究了具有漂浮结构的拖曳罐与非线性波浪的相互作用。Hu[10]等研发了一套流体动力学(CFD)程序AMAZON-SC，并数值研究了规则波驱动下波能转换器的运动响应规律。Jeong[11]等发展了一种基于改进标记密度法的数值模型来模拟浮体和自由表面的运动，计算的固定圆柱阵列与波浪的相互作用以及海上平台的响应幅度与研究数据相符。鉴于有网格方法在处理流固耦合问题时网格易在流体自由面与固体交界面处产生畸变，需对网格进行调整或者重构，耗费计算资源且影响计算精度。相较之下，无网格方法将固体和流体均离散成粒子，粒子之间的拓扑关系可不受物性参数的影响。因此，由粒子构成的物质界面随时间更新时可避免网格的调整或者重构，在提高计算效率的同时还可消除物质传递带来的不利影响。光滑粒子流体动力学方法(Smoothed Particle Hydrodynamics，SPH)以及移动粒子半隐式方法(Moving Particle Semi-implicit，MPS)为目前较为流行的两种无网格方法。Omidvar[12]等使用具有嵌入式黎曼求解器的任意拉格朗日-欧拉公式实现可变粒子质量分布，开发了一种可变质量分布的SPH方法，并数值研究了二维楔形和三维圆锥形自由浮体随波浪的运动规律。曹文瑾[13]等在MPS中建立了刚体与流体之间的单向及双向耦合模型，并对活塞式造波池中长板型浮体运动进行了数值分析。

研究拟采用SPH方法建立一种适用于研究活塞式造波驱动下浮体运动的数值模型，并获取箱式浮体六自由度条件下运动数据的变化曲线，发掘波浪驱动下箱式浮体运动的响应规律，为浮体运动的数值研究提供一个较为适用和有效的手段。

1 数值模型

1.1 SPH方法简介

光滑粒子流体动力学(SPH)是一种拉格朗日无网格方法,使用一组粒子离散化连续介质，用于流体动力学模拟时，SPH根据每个粒子周围的物理特性，在每个粒子的位置对离散的Navier-Stokes方程进行局部积分。每个粒子有一个相邻粒子的集合，基于距离的函数决定，相关的特征长度或平滑长度通常用h表示。通过在每个时间步长计算出每个粒子新的物理量，然后每个粒子根据更新的值发生移动。

SPH使用基于插值函数的积分方程将连续流体动力学的守恒定律从其偏微分形式转化为适合于基于粒子模拟的形式。SPH函数核近似的表达式为：

(1)

式中，r为粒子的矢径；下标a表示流体粒子；下标b为该粒子的相邻粒子；mb为相邻粒子的质量;ρb为相邻粒子的密度;W(r,h)为核函数，表达式为：

(2)

1.2 控制方程

粒子移动满足连续介质的动量守恒方程为：

(3)

式中，v为粒子速度矢量；ρ为粒子密度;p为压力梯度;g为重力加速度;Γ为粘性耗散项;这里采用层流粘度和亚粒子尺度(Sub-Particle Scale，SPS)湍流模型[14]。

SPH求解流体运动时整个流场的密度是弱可压缩的，即每个粒子的质量保持不变，但粒子相关的密度会产生波动。这些密度变化是通过求解连续性方程来计算的：

(4)

式中，vab为粒子间的速度矢量差。

在流体被视为弱可压缩性时，需要引入状态方程来确定基于粒子密度的流体压力。压力和密度之间的关系遵循以下表达式：

(5)

时间步长控制方程采用可变步长格式，满足(Courant-Friedrichs-Lewy,CFL)条件，表达式为：

(6)

式中，|fa|为单位质量力;cs为瞬时声速;η为阻尼系数;CCFL=0.2为CFL常数。

浮体运动是通过浮体与流体粒子的相互作用力来驱动，相互作用力的求解主要通过以下步骤实现： ①假设浮体是刚性的；②根据指定的核函数和平滑长度，判定浮体边界粒子所有的相邻流体粒子；③计算作用在每个边界粒子上的净力并累积之和为浮体所受的作用力。每个边界粒子k承受的单位质量力为：

fk=∑fka,

(7)

式中，fka为流体粒子a在边界粒子k上施加的单位质量的力。

使用刚体动力学的基本方程式求解浮体的运动，表达式为：

(8)

(9)

式中，M为浮体的质量；V为浮体的速度矢量；mk为边界粒子k的质量;I为浮体的转动惯量;ω为浮体转动的角速度矢;rk为边界粒子k的矢径;R0为浮体的质心。

在求得每一时间步的浮体运动的速度和角速度矢基础上，还需对浮体的边界粒子k的速度矢量进行更新，计算式为：

vk=V+ω×(rk-R0),

(10)

2 算例验证

2.1 活塞式造波池中的波浪传播

图1 活塞式造波池计算模型

造波板往复运动计算式为：

(11)

式中，H为波高；T为周期；d为静水深度；L为波长;S0为活塞行程;e(t)为造波板t时刻的位移。

生成的二阶斯托克斯波水位计算式为：

(12)

式中，x为水平方向坐标；η(x,t)为t时刻的水位。

波浪水位SPH数值计算结果与二阶斯托克斯波解析解的对比如图2所示。图2给出了活塞行程S0=0.165 832 m，波高H=0.15 m，周期T=2.0 s，波长L=4.523 67 m条件下，x=6.0 m处波浪水位数值解与斯托克斯波解析解的对比结果。从图2中可以看出，波浪水位的数值解与斯托克斯波水位解析解在4 s之后重合度较高。进一步提取速度测点(x=6.0 m，z=0.4 m)x与z方向速度并与斯托克斯波速度解析解进行对比，结果如图3所示。由图3可知，测点x与z方向速度数值解与斯托克斯波速度解析解4 s之后同样高度吻合。需要说明的是在图2和图3中，前4 s数值解与解析解存在明显差异，原因是数值模拟时，波浪从生成到传播至测点有一段距离，在这段时间内测点水位和速度变化幅度很小，等波浪传播过来才会产生明显水位波动，进而产生相应的x与z方向速度。综合以上说明，研究数值模型数值模拟活塞式造波池中的波浪传播过程精度较高。

图2 波浪水位SPH数值计算结果与二阶斯托克斯波解析解的对比

图3 测点x与z方向速度SPH数值计算结果与二阶斯托克斯波解析解的对比

2.2 水平圆柱自由跌落入水

以水平圆柱自由跌落入水作为算例来验证数值模型模拟流体驱动浮体运动的正确性和适用性。水平圆柱自由跌落入水的计算模型如图4所示。计算域x方向长10.0 m，z方向高16.0 m，静水深度为14.0 m，水的密度ρl=1 000 kg/m3。水平圆柱直径D=1.0 m，密度ρs=1 200 kg/m3。z方向重力加速度大小为9.81m/s2。初始时刻圆柱质心与静水面高差为0，释放后圆柱穿过自由液面完成入水。

图4 水平圆柱自由跌落入水计算模型示意图

以自由液面为z方向坐标原点，分别提取水平圆柱质心z方向位移和速度随时间的变化曲线，如图5和图6所示。从图5、图6中可以看出，由于圆柱密度大于水的密度，圆柱入水过程中不断下降，质心z方向位移随时间增加而一直呈现减小趋势。反观圆柱质心z方向速度随时间的变化曲线，呈现先减小后些许增加再些许下降的振荡趋势。将水平圆柱质心z方向位移和速度的SPH数值模型计算结果与Fekken[15]的数值计算结果对比可以看出，位移与速度曲线的吻合度均较高，说明本研究中数值模型对于流体驱动浮体运动的模拟具有较高的精度和适用性。

图5 水平圆柱质心z方向位移随时间的变化曲线 图6 水平圆柱质心z方向速度随时间的变化曲线

3 活塞式造波池中箱式浮体的运动

活塞式造波池中箱式浮体运动计算模型如图7所示。计算模型的几何尺寸为长×宽×高=16.0 m×6.0 m×8.0 m。活塞式造波板尺寸为长×宽×高=0.3 m×6.0 m×7.0 m，离计算模型左边界0.7 m。造波板行程为1.0 m，其水平往复周期运动过程为先向x正方向运动0.5 m后返回原位置再向x负方向运动0.5 m再回到原位置，其运动周期为10/3 s。正方体箱式浮体位于离左边界10.0 m的水池中央，尺寸为长×宽×高=2.0 m×2.0 m×2.0 m，密度为ρs=1 250 kg/m3。造波池中的静水深度为4.0 m，水的密度为ρl=1 000 kg/m3。z方向重力加速度大小为9.81 m/s2。初始时刻计算模型的流场速度为0，箱式浮体运动时长为16.7 s。

图7 活塞式造波池中箱式浮体运动计算模型及几何尺寸

为对波浪驱动下箱式浮体运动做深入剖析，现将箱式浮体六自由度运动描绘出来，如图8所示。其中，沿3个正交坐标轴的平移运动包括垂向、横向和纵向运动；绕3个正交坐标轴的旋转运动包括偏航、俯仰和横摇转动。进一步，将箱式浮体六自由度运动随时间变化的数据进行，位移和转动结果分别如图9和图10所示。从图9中可以看出，在波浪驱动下，箱式浮体纵向和垂向位移曲线均呈现规律的波动状态，说明浮体位移对波浪的响应以纵向和垂向为主，横向运动受波浪驱动的影响较小。由于计算模型右边界没有进行消波处理，活塞持续水平往复运动产生的波浪激励逐步增强，使得箱式浮体纵向位移振幅呈现明显增加趋势。而垂向的浮力作用主要由密度差决定，所以箱式浮体垂向位移振幅的增加趋势并不明显。箱式浮体绕正交坐标轴的旋转运动角度曲线如图10所示。由图10中可以看出，绕y轴俯仰转动以及绕x轴横摇转动的角度曲线均呈现规律的波动状态，说明箱式浮体转动对波浪的响应以俯仰转动和横摇转动为主，偏航转动影响受波浪驱动的影响较小。同样，受波浪激励逐步增强的影响，箱式浮体绕y轴俯仰转动的角度振幅呈现逐渐增强趋势，绕x轴横摇转动角度振幅的增加趋势并不明显。

图8 箱式浮体六自由度运动示意图

图9 箱式浮体沿正交坐标轴平移运动位移 图10 箱式浮体绕正交坐标轴旋转运动角度

4 结论

研究通过SPH方法建立了一个用于研究活塞式造波池中浮体运动的数值模型，并通过活塞式造波池中的波浪传播以及水平圆柱自由跌落入水分别验证了该数值模型在活塞式造波以及流体驱动浮体运动两个方面的正确性和适用性。在此基础上，对活塞式造波池中箱式浮体的运动进行了数值研究，并对箱式浮体六自由度运动的数据和响应规律进行了深入剖析，得到了以下结论：箱式浮体位移对波浪的响应以纵向和垂向为主，横向运动受波浪驱动的影响较小；箱式浮体纵向位移振幅易受波浪激励的影响，而垂向位移振幅受波浪激励的影响较小。箱式浮体转动对波浪的响应以俯仰转动和横摇转动为主，偏航转动受波浪驱动的影响较小；箱式浮体绕y轴俯仰转动的角度振幅易受波浪激励的影响，而绕x轴横摇转动角度振幅受波浪激励的影响较小。