浮体吃水比对窄缝流体运动特性的影响

曹学珅，陈林烽，孙士艳

（江苏科技大学船舶与海洋工程学院，江苏 镇江 212003）

0 引 言

FPSO 和FLNG 等海洋平台发展迅速，将石油，天然气卸载到穿梭船的过程中，FPSO 和FLNG 平台将与穿梭船组成多浮体系统，平台和穿梭船间形成窄缝。超大型浮体结构[1]由许多独立模块组成，是多浮体系统，模块间存在窄缝。多浮体系统窄缝内的流体在特定入射波条件下会发生振荡甚至共振，对窄缝两侧结构施加作用力，影响操作过程中结构稳定性。这种现象涉及的问题被称为窄缝流体共振问题。双体船两船体间也存在窄缝，设计双体船时，也需考虑窄缝流体共振。

窄缝内流体发生振荡，多浮体系统和波浪间存在复杂的水动力干扰现象[2]。特定入射波下，窄缝内流体的波幅比入射波波幅高，甚至达到入射波波幅的5 倍[3]。窄缝流体共振会影响系统稳定性，该问题得到了广泛关注。大量的理论分析和试验研究了窄缝流体共振的机理和水动力特性[4]。对此问题的研究方法有三种：理论分析、物理实验和数值模拟。

2001 年，Molin 用线性势流理论模型得到窄缝流体共振问题共振频率的求解公式。Kristiansen等[5]发现与试验结果比，线性势流模型会高估窄缝流体共振中流体的振荡响应幅度。为纠正线性势流模型缺陷，Huijsmans等[6]提出在浮体间窄缝加刚性盖，Newman[7]提出在窄缝加柔性盖，Chen[8]提出把人工阻尼引入线性势流模型。

Tan 等[9]为研究浮体舱底形状对窄缝共振的影响进行了多次试验，发现尖锐形状舱底的能量耗散比圆滑舱底大。Zhao 等[10]用瞬态波浪研究三维条件下窄缝流体共振，分析其线性阻尼项。后来Wang[11]重复Zhao的试验，发现瞬态波浪组比典型入射波浪更适合窄缝流体共振的研究。

计算流体动力学数值模拟已成为一种高效研究窄缝流体共振问题的新方法。Moradi 等[12]在考虑粘性情况下模拟二维旁靠驳船的窄缝流体共振，发现船舱底部半径增加,波浪响应幅度下降。Lu 等[1]用修正势流模型和粘性流体方法研究二维窄缝流体共振。Tan 等[13]得出修正势流模型中阻尼系数和波浪耗散率之间的关系。Lin 等[14]基于非定常雷诺平均Navier-Stokes 方程直接求解FLNG-LNG 旁靠浮体模型的水动力学载荷。

本文的结构如下：第一章介绍粘性流动模型的控制方程，VOF 模型定义，数值波浪水池的边界条件，介绍5种算例的具体设置，展示算例草图和算例网格的加密方式。第二章中使用已有试验和以往研究的数据验证粘性流动模型，利用吃水比为1 的多浮体系统讨论水池深度对模拟的影响。讨论浮体吃水对窄缝流体共振的波浪响应幅度和共振频率的影响；提出基于控制体的理论力学模型，分析窄缝流体共振的能量转换；展示不同吃水比的多浮体系统在不同入射波下的涡量图，讨论浮体吃水对窄缝流体共振中浮体周围和窄缝入口处的流体流动和能量转换的影响。第三章对本文工作进行总结。

1 数学模型

数值水池中水气两相流的粘性流体模型[15]为不可压缩流体的Navier-Stokes方程。本文所采用的数值工具为基于Open FOAM开发的Waves2Foam[16]开源程序，Waves2Foam结合水波理论为数值水池生成入射波、吸收远场反射波。

1.1 粘性流体模型

不可压缩流体的粘性流体模型的控制方程为

式中，ui和xi分别代表沿i方向的速度分量和笛卡尔坐标系分量，prgh表示超过静水压力的压力值，prgh=p-ρg(x2-xr)，其中g是重力加速度，xr表示定义在海平面的参考坐标，ρ，υ分别表示流体密度和运动学粘性系数。

1.2 流体体积函数

本文采用Hirt 和Nichols 提出的VOF 方法[17]（流体体积法）捕捉粘性流体流动模型的自由液面波。对于任意一个计算单元，用α表示流体的体积分数，则

VOF方法中流体体积的输运方程遵循以下对流方程：

式中，vr=vwater-vair表示水和空气之间的相对速度即压缩速度，式（4）中第3项为数值求解过程中引入的一个虚拟压缩项，用于改善自由液面的数值耗散。根据流体体积分数，流体密度ρ和有效动力粘性系数μ定义为

其中下标中的w和a分别表示水和空气。

1.3 算例设置

本文所使用的数值水池是一个18.09 m×1.5 m 的长方形水箱，如图1 所示。数值水池中固定的两个矩形浮体被定义为浮体A1和A2。两浮体有相同的宽度B，B=0.5 m,两浮体间有宽度为Bg的窄缝，Bg=0.05 m。两浮体吃水深度分别用h1和h2表示，浮体吃水深度的最小值是0.25 m。算例浮体吃水比值（h2/h1）的范围是2.0，1.5，1.0，0.67，0.5。数值水池的水深H是1 m。浮体下和窄缝内流体所在的空间被定义为控制体，用CV 表示。一阶Stokes 规则波在数值波浪水池的入口处生成，波浪升高幅度ζ、入射波势能φ和速度分量u，w的公式为

图1 数值模拟水池Fig.1 Schemathic of the simulation tank

式中，A表示波浪幅值（在本文中A被设置为常量0.012），k是波数，ω表示固有圆频率。在模拟中，方程（8）～（9）被用来定义速度的边界条件，入口处的压力边界条件[18]设定为法向零梯度；在出口，速度被设定为法向零梯度条件，压力设定为0；数值水池的底部和浮体的固体壁面被设定为无滑移（法向无穿透，切向无滑移）条件。在数值水池的顶部边界，总压力（动压力和静压力之和）恒定为零，允许空气从该区域流进流出。

为消除反射波的非物理方面的影响，数值水池两端设置了两个数值波浪吸收区域来吸收反射波，反射波的吸收通过使用显示松弛区技术来实现，公式为

式中，wR是速度场的计算解和具有目标解的指示场之间的加权函数。这个方法可以在计算压力速度耦合之前根据公式（12）矫正α场和（u，w）。两个消波区的长度都设置为2 个波长，可以有效地消除虚假反射。由四边形单元和三角形单元组成的非结构化网格用OpenFOAM 的blockMesh工具构建，对于每一个算例，网格的设置都是在波长方向上分布有1120 个单元，在波高方向上有220 个单元，这足以精确捕捉水和空气的分界面。由于要考虑浮体周围流体粘性的影响，对浮体周围的单元进行了局部加密，加密网格的放大草图如图2所示。

图2 数值波浪水池的加密网格Fig.2 Refined mesh of numerical water tank

在以前同吃水双浮体窄缝间流体运动研究基础上，本文改变前后浮体吃水比，探索前后浮体吃水比对窄缝流体运动特性的影响。本文讨论了5类不同吃水比的算例，算例的具体细节如表1所示。每一个吃水比情况下考虑7种不同频率的入射波，入射波波高均为0.024 m，周期分别是1.28 s、1.2197 s、1.2114 s、1.187 s、1.1762 s、1.1505 s和1.1108 s。

表1 5组不同算例的参数Tab.1 Parameters for five different sets of cases

2 结果和讨论

2.1数值模拟验证

数值水池中有四个波高计被用来记录波浪的升高幅度，这四个波高计的位置分别是x=2.705 m，3.205 m，4.73 m 和6.255 m（入口的位置在x=-4.02 m 处），第三个波高计位于窄缝的中间。本文通过5组算例来调查浮体吃水对多浮体系统的窄缝共振现象的影响，就窄缝内液面的升高幅度、窄缝流体共振频率（窄缝流体振动幅度达到最大值时的流体振动频率）、浮体周围及窄缝中的能量耗散和能量转换进行讨论。

在进行数值模拟之前，通过数值造波与线性波理论获得的波浪高度曲线进行对比，以验证数值造波的可靠性，图3 显示的是使用数值模拟和线性波理论获得的波浪高度随时间变化的曲线，可以看到，在数值模拟达到稳定状态后使用数值模拟造波获得的波浪高度曲线与理论曲线相吻合。

图3 由数值模拟和线性波理论获得的波浪高度随时间变化曲线Fig.3 Time histories of the wave elevations obtained by the numerical simulation and the linear wave theory

本文用三套网格测试网格收敛性，三套网格的单元数分别是186 450（网格A），246 400（网格B）和302 300（网格C）。图4 显示的是浮体吃水比为1 的算例在入射波周期为1.1108 s，三种网格下的窄缝内波浪高度随时间变化曲线，结果表明单元数为246 400（网格B）的网格收敛性好，该级别的网格足以用来模拟进一步的研究。

图4 吃水比为1，入射波周期为1.1108 s三种网格得到的窄缝间波浪幅度随时间变化曲线Fig.4 Time histories of the wave elevation in the gap obtained using three different meshes for h2/h1=1 and T=1.1108 s

所有算例的计算时间为50 s，包含40多个波浪周期。图5是吃水比为1的窄缝内波浪幅度随时间变化曲线，表明窄缝内的流体的运动达到动态稳定需要花费一段时间，在20 s后能达到动态稳定。曲线表明窄缝内流体运动达到动态稳定后，其运动非常接近周期性。

图5 吃水比为1的窄缝内流体波浪幅度随时间变化曲线Fig.5 Time history of the wave elevation of the fluid in the gap for h2/h1=1

算例1 中，浮体吃水比值是1。将算例1 数值结果与Saitho 的试验及Zhang 等[20]研究结果比较，结果如图6 所示。结果表明，算例1 模拟结果与试验结果及他人的数据吻合良好。需注意的是，本文采用的入射波的周期与试验和Zhang的不同，所以，在图6中曲线无法完全吻合。将本文数据与Zhang的数据对比可以发现，在一些横坐标相近的点处，纵坐标很相近，数据间差值不超过0.9%，且数据变化趋势也一致，因而可以说明本文数据与Zhang的数据吻合，计算精度与Zhang的相近，也可以进一步说明本文得到的数据与试验吻合且可靠。因Saitho 的试验水池水深为0.5 m，所以算例1 中数值水池的水深设为0.5 m。算例1 模拟结果与试验及他人结果的良好吻合表明，基于粘性流体理论通过Open-FOAM 建立的数值水池能对窄缝流体共振精确模拟，结果可靠。为更好比对数据，后处理入射波频率ω和窄缝内波浪响应幅度ηA，处理方式遵循式（11）～（12）。

图6 粘性流体模拟结果与他人研究结果对比Fig.6 Comparison of predictions of viscous model and other studies

2.2 流体振荡幅度

因后续算例浮体吃水深度加深，后续算例中数值水池水深设为1 m。对比算例1在水深为1 m和0.5 m两个条件下得到的结果，结果如图7所示。由图表明，数值水池水深影响窄缝流体共振的波浪幅度和共振频率，随着数值水池水深从0.5 m增加到1 m，窄缝流体共振的波浪幅度明显下降，且水深为1 m的算例窄缝流体共振频率比水深为0.5 m的大。图8 显示5 类不同吃水比的算例在水深为1 m、7 种不同频率入射波下，窄缝内的流体振荡幅度随入射波频率变化的曲线。

图7 0.5 m水深和1 m水深的窄缝流体振动幅度对比Fig.7 Comparison of wave amplitudes of gap resonance with water depths of 0.5 m and 1 m

图8 不同吃水比下的窄缝流体振荡波幅曲线Fig.8 Fluid oscillation amplitude profiles against angular frequency of incident wave for different draft ratios

图8中吃水比（h2/h1）为2.0、1.5的曲线是背浪侧浮体吃水变大的窄缝流体振荡曲线，与吃水比为1的曲线对比发现，背浪侧浮体吃水加深，窄缝流体振荡波幅变大，尤其是在低频率阶段，共振频率减小。吃水比（h2/h1）为0.67、0.5 的曲线是迎浪侧浮体吃水变大的窄缝流体振荡曲线，与吃水比为1的曲线对比发现，迎浪侧浮体吃水加深，窄缝流体振荡波幅显著下降，共振频率也下降。这与试验观察结果[21]一致。由于吃水比为2.0 和1.5 算例的控制体与吃水比为0.5和0.667的近似，所以用来计算固有频率的附加质量（由等效活塞引起的）可能也是近似的。因而可以发现，共振频率从吃水比为1.0到2.0的变化趋势近似于从吃水比为1.0 到0.5 的变化。波浪振荡幅度的变化可以用透射系数和反射系数的变化解释。

不同吃水比的多浮体系统的窄缝流体在不同频率入射波下的振荡频率ωf通过快速傅里叶变换（FFT）得到，如图9 所示。图9 表明，随着入射波频率ω*的增加，窄缝内流体振荡的频率增加，但窄缝流体振荡的幅度没有随振荡频率的增加而增加。窄缝流体的振荡幅度是先增加，在某个特定的入射波频率下达到最大值后，窄缝流体振荡的幅度开始下降。这表明，不同吃水比的多浮体系统在特定的入射波频率下也会发生窄缝流体共振现象。

图9 窄缝流体振荡的频率随入射波频率变化的曲线Fig.9 Variation of the frequency of fluid oscillation in the narrow gap with angular frequency of incident wave

2.3 透射波和反射波系数

图10是三种不同吃水比算例的透射波（Kt）和反射波（Kr）系数随入射波频率变化曲线。透射波系数是透射波波高（Ht=2At）与入射波高（Hi=2Ai）之比，Kt=Ht/Hi，反射波波高（Hr=2Ar）和入射波波高之比是反射波系数，Kr=Hr/Hi。透射波波高使用浮体后的波高计（波高计4）测量得到，反射波系数使用由Sun等[22]提出的方法计算得出：

图10 不同吃水比下透射波和反射波系数随入射波频率变化曲线Fig.10 Transmission and reflection wave coefficient profiles against the incident wave frequency for different draft ratios

式中，Δx表示波高计1 和2 之间的距离，ξ（x1，t）和ξ（x2，t）是分别由波高计1 和2 测得的复合波面升高高度，可以使用希尔伯特（Hilbert）变换来计算,具体细节请查阅Sun等[22]的文献。

图10 表明随着入射波频率的增加，透射波系数变小，这与Jiang 等[23]的研究结果一致，也可以看到，迎浪侧浮体和背浪侧浮体吃水的增加都会使透射波系数减少。因为吃水比为0.5 和2.0 算例的两个浮体的总吃水深度都是一样的，所以这两个算例的的透射波系数很相似。在反射波系数方面，可以看到当迎浪侧浮体的吃水大（吃水比为0.5）时，反射波系数很大，在此算例中，从波浪转移到窄缝流体的能量减少，因此窄缝内流体响应相对较弱。当吃水比为1.0 时，反射波系数下降。当背浪侧浮体吃水加深（吃水比为2.0）时，在低入射波频率段，反射波系数大幅减小，这表明更多的能量从入射波转移到窄缝内流体。由于有更低的透射波系数，故吃水比为2.0 算例的窄缝流体振荡幅度比吃水比为1.0的更大。透射波和反射波系数的变化曲线既反映了从入射波向窄缝间流体的能量输运，也与窄缝间流体的振荡幅度变化趋势一致。而且可以看到，在每个吃水比下反射波系数最小的频率对应于窄缝间流体出现最大振荡幅度的地方。

2.4 粘性耗散

粘性耗散也对窄缝内流体振荡有影响。当流体振荡发生在窄缝间隙时，由于流体粘性，浮体周围剪切力会导致能量耗散，使从波浪转移到窄缝内流体振荡的能量减少。对于不同吃水比的浮体系统，窄缝流体的响应可能会导致不同大小的剪切力，导致粘性耗散也会不同。为说明浮体吃水对窄缝间流体振荡的影响，得到了共振频率下的涡量图，如图11所示。

图11显示的是h1=0.5 m时在共振频率下两浮体周围的涡量图。由于此算例大部分的入射波能量都被反射了，窄缝流体振荡的幅度很小，较小的响应可能导致浮体周围产生较小的剪切力，因此可以看到，在此算例中，只在背浪侧浮体拐角处附近有少量的涡，这说明在h1=0.5 m 时能量耗散少。当浮体的吃水h1=h2=0.25 m（图11）时，窄缝内流体振荡的幅度增加，更强的流体响应导致浮体附近产生的剪切应力变大，在两个浮体下方形成了更多的涡，揭示了h1=h2=0.25 m 的算例有更多的粘性耗散；当背浪侧浮体吃水增加到h2=0.5 m（图11）时,窄缝流体振荡的幅度达到最大值，流体振荡的最大响应导致浮体周围产生最大的剪切应力。因此可以清楚看到在间隙下方产生了很多具有大涡量的涡，也可以说在这种情况下产生了最大的粘性耗散，这个结果也通过图10最小的透射波和反射波系数得到了验证。总之，当窄缝内流体最大振幅随吃水比h2/h1增加时，浮体附近的流体相对运动增强，粘性耗散也同步增强。

图11 不同吃水比情况下浮体外部涡量分布云图Fig.11 Cloud chart of vorticity distribution outside the floating body with different draft ratios

3 结 论

本文使用粘性流体理论和VOF 模型，通过耦合于OpenFoam 的Wave2Foam 建立算例研究了浮体吃水比对窄缝流体振荡的影响。为了验证粘性流体力学理论和VOF 模型的可靠性，将数值结果与试验及他人结果进行了比较，发现吃水比为1 的算例结果与试验和理论结果吻合，验证了Waves2Foam工具包的粘性理论和VOF模型的准确性。

数值水池的水深对数值结果有明显影响。数值水池的水深加深，窄缝流体共振的波幅变小，共振频率变大。同时发现浮体吃水比对多浮体系统窄缝流体共振的波浪幅度和共振频率有非线性影响。浮体吃水增加，共振频率下降；迎浪侧浮体吃水增加，导致反射波系数增加，流体振荡波幅显著下降；背浪侧浮体吃水增加，导致反射波系数下降，流体振荡波幅上升。同时，由流体粘性引起的粘性耗散随窄缝流体振荡幅度的增加而同步增加。最后，对涡量演化的分析说明了窄缝流体振荡幅度的变化。