0.   引　言

两船并靠作业是常见的海上作业方式，包括大型军舰补给、失去动力船舶的海上救援和浮式生产储油船（FPSO）及圆筒型FPSO（CFPSO）过驳作业等。当两船靠帮作业时，由于两船间距较近，波浪作用下会产生水动力干扰，容易引发间隙内大幅波浪振荡和两船大幅运动响应，威胁作业安全。因此，研究波浪中两船并靠响应具有重要的工程意义。

早期对间隙共振问题的研究主要基于线性势流理论。例如，Molin^[1]基于线性势流模型分析了二维并列布置和三维矩形月池的共振频率和自然模态；丛龙飞等^[2]采用特征函数展开法对窄缝内流体的共振频率开展了解析；朱传仁^[3]、Li等^[4]、Li等^[5]和Sun等^[6]均基于经典势流模型进行了相关研究。为了更深入地研究间隙共振现象并验证理论模型的准确性，研究者们进行了大量的物理实验。例如，Saitoh等^[7]针对二维两固定浮体间隙内流体的振荡进行了物理实验，揭示了模型几何参数对间隙内流体共振的影响；Zhao等^[8]对聚焦波作用下双船间的间隙共振进行了模型实验；赵文华等^[9]采用三维物理模型实验方法，研究了浮式液化天然气装置（FLNG船）与液化天然气运输船（LNG船）系统在波浪作用下的运动幅值，给出了特定船型在某一波浪条件下的水动力响应。刘春阳^[10]和Ning^[11]进行了二维实验研究，探讨了双浮体的宽度比和吃水深度比对共振频率和共振幅值的影响。

研究表明，无论是二维^[12-13]还是三维^[14-15]线性势流理论，往往都会高估间隙间共振响应幅值的峰值。因此，许多改进的势流模型被提出以提高计算值和实验值的吻合度，例如增加刚性盖板^[16]或者柔性盖板^[17]以及在自由液面增加阻尼^[18-22]等方法，但是阻尼系数往往需要通过实验测试或者计算流体力学（CFD）数值模拟方法获得^[14,23]。

CFD方法能够充分考虑流体的黏性效应，因此，采用CFD方法预测的双船间隙共振响应与模型试验数据对比良好。同时，CFD模拟的一大优势在于其能够直观地显示流体运动，有助于理解间隙共振问题中黏性的作用机理。在二维条件下，Lu等^[25]基于N-S方程，使用开源软件OpenFOAM，通过流体体积（VOF）追踪自由液面的方法，对规则波作用下固定双浮体以及三浮体间的间隙共振问题进行了数值模拟。Moradi等^[26]系统研究了入口参数对间隙共振的影响，Gao等^[28]采用OpenFOAM研究了规则波和聚焦波群作用下两个固定浮体之间的间隙共振。Meringolo等^[29]着眼于分析共振过程中的能量耗散，采用光滑粒子流动力学（SPH）模型对固定浮体与直墙间间隙内流体共振过程中的能量耗散进行了深入数值分析，旨在更全面地揭示该过程的特性。这些研究表明，CFD模拟预测的数值结果与现有实验数据吻合较好，但相对于势流模型，CFD模拟计算耗时较长。因此，为了提高计算效率，势−黏流耦合方法被广泛应用于计算中。

同时，间隙内流体共振与浮体的位置和运动状态有关。张婧文^[30]对同主尺度固连双浮体窄缝间水体共振进行了研究，揭示了不同吃水与不同间隙宽度对固连双浮体间隙共振的影响。Gao等^[31]研究了一浮体固定而另一浮体自由垂荡的情况，探讨了上游浮体运动对间隙共振的影响。He等^[32]研究了一浮体运动而另一浮体固定工况中的高次谐振。

考虑到浮体运动对间隙共振现象影响的研究较少，且此类研究鲜有对浮体尺度的分析，因此，本文将使用双向耦合方法，对比分析不同尺度与不同迎浪工况下双浮体在不同入射波频率下的垂荡运动响应情况。在研究之前，进行同尺度双浮体数值模拟，并与已有实验结果进行对比，验证数值方法的准确性。

1.   数值水池的建立

本文采用Meng等^[33]的出的场−域分解势−黏流耦合方法，计算域具体划分如图1所示。图中流体计算域分为黏流计算域${\Omega _{{\text{CFD}}}}$和势流计算域${\Omega _{\rm{PF}}}$，两者之间之间存在重叠区域，即耦合区，耦合区域长度由${L_{{\text{couple}}}}$表示，黏流区域长度由${L_{{\text{CFD}}}}$表示，势流域${\Omega _{\rm{PF}}}$由自由表面${S_{\rm{FS}}}$、耦合界面${S_{\rm{CI}}}$、远场边界${S_\infty }$和池底${S_{\rm{BM}}}$包围。在黏流计算域中，求解纳维−斯托克斯（N-S）方程，采用有限体积法（finite volume method，FVM）进行控制方程离散，并使用几何流体体积VOF方法进行自由液面捕捉。在势流计算域中，引入线性自由液面假设，采用Rankine源法进行求解。势流与黏流计算域通过匹配边界和耦合区实现双向数据交换。

图  1  计算域示意图

Figure  1.  Schematic diagram of the computational domain

下载: 全尺寸图片 幻灯片

1.1   黏流计算域流动模型建立

考虑网格运动的非定常不可压缩流体的连续性方程和动量方程可表示为：

式中：$\nabla$为向量微分算子；u为速度场；${{\boldsymbol{u}}_{\rm{g}} }$为网格移动速度；$\rho$为流体密度；$\mu$为流体的动力黏度；p为压力。

动量方程采用有限体积法离散，求解器采用PIMPLE算法求解耦合的速度和压力。在求解控制方程之前，采用算术VOF方法（MULSE）求解流体的气体和液体的相分数：

式中：U为前一时间步求得的流场速度场；${\alpha _0}$为液体相分数；${{\boldsymbol{U}}_{\text{r}}}$为压缩速度。

流体密度$\rho$和动力黏度$\mu$的求解公式为：

式中：${\alpha _{\rm{m}} }$为相分数；${\rho _{\rm{m}}}$为密度；${\mu _{\rm{m}} }$为动力粘度。

本文求解器采用OpenFOAM中的动态重叠网格方法模拟双浮体的运动，并使用cellVolumeWeight方法进行网格间物理量的插值及传递。

1.2   势流计算域流动模型建立

远场的势流计算域假设流体无黏、不可压缩且流动无旋。如图1所示，势流域${\Omega _{\rm{PF}}}$由自由表面${S_{\rm{FS}}}$、耦合界面${S_{\rm{CI}}}$、远场边界${S_\infty }$和池底${S_{\rm{BM}}}$包围。因此，可以引入速度势${\phi _{\rm{T}}}(x,y,{\textit{z}},t)$，其的控制方程为

入射波势通过使用波浪模型（例如斯托克斯波理论）显式求解，而绕射势通过求解边界积分方程${\phi _{_{}D}} = {\phi _{}}_{\rm{T}} - {\phi _{}}_{\rm{I}}$计算，其中${\phi _{_{}D}}$为扰动势；${\phi _{}}_{\rm{I}}$为入射波势。

平均自由表面上的线性运动学和动力学自由表面边界条件${S_{\rm{FS}}}$为：

式中：g为重力加速度。

在耦合界面${S_{\rm{CI}}}$上，应用Neumann边界条件：

式中：${{\boldsymbol{v}}_{\rm{CI}}}$为通过插值从CFD域获得的速度；n为指向流体域外部边界的单位法向量。CFD域求解的总流场和${{\boldsymbol{v}}_{\rm{CI}}}$仅包含绕射部分，因此在传输至耦合界面边界${{S} _{\rm{CI}}}$之前，需将入射波部分从总CFD解中分离。

池底采用不可穿透边界条件：

在边值问题中，远场边界${{S} _\infty }$处，由于辐射条件，$\phi$和$\dfrac{{\partial \phi }}{{\partial {\boldsymbol{n}}}}$为零。在截断的自由表面边界上实现数值海滩，以避免反射波，从而满足辐射条件。自由表面边界条件修改如下：

其中，${\mu _{\rm{b}}}$为阻尼强度，沿着数值消波区域上的径向距离逐渐增大。

式中：${\mu _0}$为常数；${r_{\rm{D}} }$为位于数值海滩中某点的径向距离；${r_0}$为数值海滩的起始半径；L为数值海滩的长度。采用边界元法（BEM）求解上述边界值问题。

势流计算域采用常值面元法进行离散，自由液面采用Adam-Bashforth方法进行时间步进。

1.3   势流−黏流双向耦合模型

首先在重构的自由液面网格中进行搜索，获得耦合边界上各面元中心点水平位置$(x,y)$处的自由液面升高${\eta _{\rm{CFD}}}(x,y)$，其次通过Wheeler拉伸方法，将耦合边界网格进行拉伸，使其与黏流计算域中的自由液面重合。

再次，通过在黏流计算域内插值，得到黏流场内耦合边界自由液面下一系列速度值。最后将插值得到的流场速度进行入射波速度分离，得到耦合边界处的边界条件。

由于势流计算域和黏流计算域的空间离散方法和离散尺度均不同，因此在势流计算域建立虚拟插值网格，其范围略大于耦合区的尺寸。

虚拟插值网格流场信息准备完毕后，采用松弛区方法，通过权重函数将势流场的数据传输至黏流场。

式中：$\varphi$为流场参数，包括流场速度和液相体积分数；${\varphi _{\rm{PF}}}$为势流场参数；${\varphi _{{\text{CFD}}}}$为黏流场参数；$\omega$ 为权重函数，取值范围为0～1。

式中：幂指数p的取值为3.5（默认设置）；$\sigma$为松弛区内的局部坐标值，根据松弛区形状而定。

2.   数值计算结果

2.1   数值方法验证

应用双向耦合数值模型模拟两个固定并排浮体之间的间隙共振问题，并通过与Saitoh等^[34]的实验结果进行比较以验证本数值模型的准确性。模型设置水深h为0.5 m，两个尺寸相同的浮体位于水池中，浮体宽度B为0.5 m，吃水T为0.252 m，两体间隙d为0.05 m。在两浮体间隙设置波高仪测量波高。入射波频率$\omega$为5.28 rad/s，波高H为0.024 m。网格划分如图2所示，在自由面附近一个波高范围内划分约12个网格，当水池中存在浮体时，对浮体周围进行网格加密，并使用重叠网格进行网格间物理量的传递。

图  2  网格划分

Figure  2.  Meshing diagram

下载: 全尺寸图片 幻灯片

图3给出了两浮体间隙波面升高随不同入射波频率的变化规律，图中两条曲线分别为实验结果与本模型数值模拟结果。其中x轴为利用水深h无因次化得到的入射波频率$\omega /\sqrt {g/h}$。y轴为以入射波幅A₀进行无因次化得到的间隙处的平均波高${A_{\rm{g}}}/{A_0}$。

图  3  间隙中心处无因次波高

Figure  3.  Non-dimensional wave height at the center of the gap

下载: 全尺寸图片 幻灯片

由图3可知，本模型数值模拟结果与实验结果吻合良好，因此本模型可以较为准确地模拟双浮体间隙共振问题。

2.2   考虑浮体运动的双浮体数值模拟

进一步开展规则波中不同主尺度双浮体垂荡运动数值模拟。数值模型水深h设置为2.5 m，CFD计算区域长度为4 m，耦合区域长度设置为2倍波长。势流计算域设置为7倍波长，其中5倍波长的势流计算域用于消波模拟。两个尺寸不同的浮体位于水池中，根据浮体大小分别命名为Body1（大）和Body2（小），两浮体的主尺度见表1。两体间隙为0.05 m。选取入射波频率为3.01～7.17 $rad/s$，波高H为0.024 m，分别计算Body1迎浪（工况1）与Body2迎浪（工况2）两种工况下两浮体的垂荡运动响应，两种工况设置如图4所示，采用24核并行计算，单个算例模拟时长40 s，计算耗时28 h。

表  1  浮体主尺度

Table  1.  Dimensions of floating bodies

浮体	宽度B/m	吃水T/m
Body1	0.5	0.252
Body2	0.25	0.126

下载: 导出CSV 

| 显示表格

图  4  工况设置

Figure  4.  Operating conditions

下载: 全尺寸图片 幻灯片

首先进行收敛性分析，设置网格大小分别为0.0054，0.0039与0.002 0 m，时间步长$\Delta t = 0.1$，取入射波频率为4.39$rad/s$，计算结果如图5所示，其中x轴$\Delta s$为网格尺寸，y轴${\zeta _{\rm{a}}}/{A_0}$为以入射波波幅${A_0}$进行无因次化得到的浮体平均垂荡运动响应幅值。结果表明，网格大小为0.003 9 m时，计算结果收敛，后续模拟均采用此网格尺寸。根据时间步长的收敛性分析，时间步长$\Delta t = 0.1$满足要求。

图  5  网格尺寸的收敛性分析曲线

Figure  5.  Convergence analysis curve for mesh size

下载: 全尺寸图片 幻灯片

图6为Body1在两种工况下的垂荡运动响应幅值曲线。由图可见，随着入射波频的额增加，Body1垂荡运动响应幅值呈现先增大后减小的趋势。在工况1下，入射波频率为4.39 rad/s时，其垂荡运动响应幅值达到最大值；在工况2下，入射波频率为4.5 rad/s时，其垂荡运动响应幅值达到最大值。由于上游浮体的遮蔽效应，在大部分入射波频率下，Body1在工况2下的垂荡运动响应幅值大于工况1下的。在相同入射波频率下，两种工况下的垂荡运动响应幅值相差不大。

图  6  无因次垂荡运动响应幅值（Body1）

Figure  6.  Non-dimensional heave motion (Body1)

下载: 全尺寸图片 幻灯片

图7为Body2在两种工况下的垂荡运动响应幅值曲线。由图可见，随着入射波频率的增加，Body2的垂荡运动响应幅值也呈现出先增大后减小的趋势。在工况1下，当入射波频率为5.19 rad/s时，其垂荡运动响应幅值达到最大值；在工况2下，入射波频率为4 rad/s时，其垂荡运动响应幅值达到最大值。不同频率下，Body2在工况1下垂荡运动响应幅值均大于工况2下的，这是由于上游浮体对下游浮体的遮蔽效应。与Body1相比，在相同入射波频率下，Body2在两种工况下的垂向运动响应幅值相差较大。

图  7  无因次垂荡运动响应幅值（Body2）

Figure  7.  Non-dimensional heave motion (Body2)

下载: 全尺寸图片 幻灯片

图8与图9分别为在工况1下，作用在Body2上的垂向力F_Z与间隙波面升高H_A随频率变化的曲线。由图可见，Body2垂荡运动响应幅值与其所受垂向力和间隙波面升高趋势相同，都在频率为5.39 rad/s时达到最大值。这是由于入射波传到Body1时产生反射波，且两浮体间存在水动力干扰，导致间隙波面升高，作用在Body1上的垂向力增加，从而增强其垂向运动。

图  8  垂向力（Body2）

Figure  8.  Vertical force (Body2)

下载: 全尺寸图片 幻灯片

图  9  间隙波面升高

Figure  9.  Wave surface elevation at the gap

下载: 全尺寸图片 幻灯片

图10和图11分别给出了入射波频率为3.01和5.28 rad/s时工况2下两浮体垂荡运动响应的时历曲线。结果表明，当入射波频率较小时，两浮体倾向于同相位运动；而当入射波频率较大时，两浮体倾向于反相位运动。

图  10  入射波频率为3.01 rad/s时两浮体垂荡运动时历曲线

Figure  10.  The timing curve of the vertical motion of the two floats at the incident wave frequency of 3.01 rad/s

下载: 全尺寸图片 幻灯片

图  11  入射波频率为5.28 rad/s两浮体垂荡运动时历曲线

Figure  11.  The timing curve of the vertical motion of the two floats at the incident wave frequency of 5.28 rad/s

下载: 全尺寸图片 幻灯片

图12展示了双浮体相对运动响应幅值随入射波频率变化的趋势。由图可见，当入射波频率从4.39 rad/s增加到6.28 rad/s时，工况1下两浮体相对运动响应幅值远大于工况2下的。这是因为迎浪与背浪条件对Body2的影响更大，而对Body1的影响较小。

图  12  无因次相对垂荡运动响应幅值

Figure  12.  Non-dimensional relative heave motion

下载: 全尺寸图片 幻灯片

通过以上结果发现，无论是Body1还是Body2，作为上游浮体，在不同频率下，其垂荡运动响应幅值均大于其作为下游浮体时的响应幅值，这是由于上游浮体对下游浮体的遮蔽效应。与Body1相比，在相同入射波频率下，Body2在两种不同工况下的垂荡运动响应幅值差异较大。此外，当入射波频率较大时，两浮体的垂向运动更趋向于反相位运动。

3.   结　论

本文采用势−黏流双向耦合与重叠网格方法，对规则波作用下不同尺度双浮体的垂荡运动进行了研究。首先对固定双浮体进行模拟，并与试验数据进行对比，验证了求解器的准确性。随后，开展了不同尺度双浮体的垂荡数值模拟，设置较小浮体迎浪与较大浮体迎浪两种工况。研究得到以下结论：

1） 上游较大浮体对下游较小浮体存在显著的遮蔽效应。

2） 当入射波频率较大时，两浮体的垂荡运动倾向于反相位运动；当入射波频率较小时，倾向于同相位运动。