在陆地资源开发利用已趋于饱和的今天，海洋资源作为地球自然资源的主要组成部分得到了世界各国的青睐，海洋工程结构物的开发及应用的竞争日趋激烈。因海洋结构物及其组成部分之间多存在着窄缝，如：靠泊的船体与码头、并排的船与船、浮式生产储油卸油装置(floating production storage and offloading，FPSO)的月池及箱式超大型海洋结构物各模块间隙等。在入射波浪的持续砰击下，间隙内的波面抬高呈往复运动。在某些特定波频率下，其振幅要远大于入射波幅，对结构物的作用力也随之显著增大，且力的变化极为剧烈。这种现象被称为窄缝共振。

窄缝共振关系到海洋工程结构物的正常生产作业安全，对结构使用寿命也存在影响。研究窄缝共振问题，不仅要深入分析现象，更要明晰其产生的机理，对于多浮体海洋结构物及FPSO等系统的设计具有积极指导意义，同时也对规避海域作业时窄缝共振的不利影响具有重要的工程应用价值。

线性理论应用到求解两箱体间的窄缝共振问题^[1-8]，学者们基于线性理论进行了变参数研究和共振模态研究。基于半解析方法和势流理论的方法虽可高效地对双箱共振问题进行求解，但因忽略了流体的粘性作用，共振区间内预报双箱间窄缝处的波面抬高常常偏差较大，且多为固定双箱，未考虑双箱在入射波作用下的运动响应。

随着计算机技术的不断发展计算流体力学(computational fluid dynamics，CFD)方法更多的应用到窄缝共振问题的研究中^[9-13]。CFD方法考虑流体粘性，可更直观地展示流场细节，便于对精细流场进行分析。Saitoh等^[14]做了不同窄缝宽度与浮体吃水的固定双箱模型实验，对不同入射波浪作用下窄缝共振波高进行了研究，发现窄缝内流体的共振波数、共振波高受箱体吃水深度及箱体间距影响。Iwata等^[15]随后采用与Saitoh相似的实验方法，研究了3个毗邻排放箱体间隙内水体的共振现象。物理模型实验的方法更加直观，与数值模拟相结合的方法能更加深入地对窄缝共振现象进行分析。

由上述文献研究可见，对于双箱共振问题的研究大多基于双箱固定不动的假定^[6-15]，与多浮体海洋结构物实际运动响应的情况相差较大，仍需进一步研究。本文对自由垂荡双箱在不同频率入射波激励下的运动响应及双箱间窄缝处水体共振问题，采用实验和数值相结合的方法研究，深入分析了窄缝水体共振机理，给出了该类物理实验改进建议。自由垂荡双箱的窄缝共振问题新颖，涉及入射波砰击、箱体自由垂荡响应、窄缝处水体活塞运动之间的耦合作用，现象复杂且影响因素多，本文对于进一步研究多浮体结构物垂荡及窄缝水体共振问题具有积极作用。

1 实验模型及设置

本模型实验水槽长68 m，宽1.0 m，深1.0 m，水槽前后两端均设有消波装置。两箱体尺寸完全一致，如图 1所示，箱体宽度B=0.5 m, 双箱间窄缝宽度B_g=0.05 m, 水深h=0.5 m, 箱体吃水d=0.153 m。

实验所用浪高仪为TS-DWG型电容式数字信号浪高仪，绝对误差≤0.1 mm，浪高仪放置于双箱间窄缝的中心，如图 1所示。浪高仪的采样频率设为50 Hz，采用电磁式非接触六分量传感器LIBERTY 240\8系统来测定箱体位移。实验方箱由壁厚10 mm的PVC硬板制成，方箱沿水槽方向宽度为967 mm，箱体刚度及水密性满足实验要求。水槽内壁粘贴有机玻璃板，两箱体沿固定在玻璃板上的导轨做垂荡运动，该导轨由耐腐蚀、低噪音的SGR10内置双轴心滚轮直线结构组成如图 1所示。

实验中对箱体位移和窄缝处的波面抬高进行实时监测。图 2为实验中的经典工况kh=2.40(共振点)时窄缝处波面抬高和箱体A垂荡运动的时历曲线，k为波数、H₀为入射波高、h_g为窄缝处波面抬高、y_A为箱体A的垂向位移。波面抬高周期性较好，波面稳定，波峰处较为平缓。箱体A垂荡运动呈周期性变化，运动达到相对稳定状态。图 3为共振频率kh=2.40条件下，一个波周期内实验的照片，图中长实线为静水面位置，短虚线为窄缝内波面，窄缝内水体的运动幅值大、变化剧烈。

2 数值模型 2.1 有限体积法

STAR-CCM+使用有限体积法求解雷诺平均纳维斯托克斯(Reynolds-averaged Navier Stokes，RANS)方程，将流体域离散化为有限个控制体积，采用压力耦合方程组的半隐式方法(semi-implicit method for pressure-linked equations，SIMPLE)实现压力和速度之间的隐式求解。积分形式的控制方程为：

质量守恒方程：

$ \frac{\mathrm{d}}{\mathrm{d} t} \int\limits_{V} \rho \mathrm{d} V+\int\limits_{S} \rho\left(\boldsymbol{v}-\mathrm{v}_{\mathrm{b}}\right) \cdot \boldsymbol{n} \mathrm{d} S=0, $

(1)

动量守恒方程：

$ \begin{gathered} \frac{\mathrm{d}}{\mathrm{d} t} \int\limits_{V} \rho \boldsymbol{v} \mathrm{d} V+\int\limits_{S} \rho\left(\boldsymbol{v}-v_{b}\right) \cdot \boldsymbol{n} \mathrm{d} S= \\ \int\limits_{S}(\boldsymbol{T}-p \boldsymbol{I}) \cdot \boldsymbol{n} \mathrm{d} S+\int\limits_{S} \rho \boldsymbol{b} \mathrm{d} V \end{gathered} $

(2)

空间守恒方程：

$ \frac{\mathrm{d}}{\mathrm{d} t} \int\limits_{V} \mathrm{~d} V-\int\limits_{S} v_{b} \cdot \boldsymbol{n} \mathrm{~d} S=0 $

(3)

式中：v是流体速度矢量；v_b表示控制体积表面的速度；n是控制体积曲面的单位法向量；控制体积的面积为S；体积为V；流体密度为ρ；时间为t；压力为p；流体应力以张量T的形式表示；I是单位张量；物体力的矢量表示为b。

通过流体体积(volume of fluid，VOF)多相模型引入体积分数来计算出多相流中的自由液面位置；使用高分辨率交界面捕捉(high resolution interface capturing，HRIC)离散格式来精确捕捉自由液面的流体变化。

2.2 DFBI模型

动态流体相互作用(dynamic fluid body interaction，DFBI)模型可用来求解刚性体的流固动态耦合，使用6自由度(6-DOF)求解器计算结构物的6自由度运动，6-DOF求解器采用二阶精度的梯形方案，可以实现6个自由度同时求解。在总体坐标系下，物体的平动方程为：

$ m \frac{\mathrm{d} \boldsymbol{v}_{\mathrm{c}}}{\mathrm{d} t}=\boldsymbol{f} $

(4)

式中：m是物体质量；f是作用在物体上的合力；v_c是质心的速度。

在随体坐标系下物体的转动方程：

$ \boldsymbol{M} \frac{\mathrm{d} \boldsymbol{\omega}}{\mathrm{d} t}+\boldsymbol{\omega} \times \boldsymbol{M} \boldsymbol{\omega}=\boldsymbol{M}_{\boldsymbol{n}} $

(5)

式中：M是惯性矩张量；ω是刚体的角速度；M_n是作用在物体上的合力矩。

转动惯量的张量展开式为：

$ \boldsymbol{M}=\left[\begin{array}{ll} M_{x x} & M_{x y} \\ M_{x y} & M_{y y} \end{array}\right] $

(6)

2.3 计算域布置和边界条件

计算域的边界条件设置如图 4所示，箱体A和B尺寸与实验箱体完全一致，流域速度入口边界通过场函数设定自由液面位置和一阶波的速度，压力出口边界给出了自由面条件，出口处的压力被设定为波浪的静水压力。在计算域的入口和出口参照实验条件使用力消波(wave forcing)进行消波处理。

2.4 网格划分

网格划分采用了切割体、拉伸体、网格重叠等技术组合，保证网格质量。具体网格结构如图 5所示，箱体周围网格逐步加密(图 5左图)以精确捕捉作用在箱体上的波浪力和箱体运动响应，远离箱体区域使用渐疏网格，减弱远方边界处的波浪反射并减少网格数量(图 5右图)，同时自由液面附近网格适当加密捕捉波面分布。在箱体周围的网格为随体运动网格，该网格在背景网格上重叠平动，计算过程中各网格形状无变化，精确预报箱体的运动响应，计算效率较高。

3 结果与讨论 3.1 单箱垂荡运动收敛性验证

本文基于该数值水池对单个箱体在波浪中的运动响应进行数值模拟，数值模型的收敛性验证见文献[16]。本节进一步研究了无因次波数kb在0.4~1.2的典型工况(b=0.5B)。为对比分析，本文数值波浪水池和箱体参数与Rodríguez的实验^[17]保持一致，所得箱体的响应幅值算子(RAO)如图 6所示：kb=0.4, 0.5, 0.8, 1.0, 1.2工况与实验值对比吻合极好，相对误差仅在5%左右；kb=0.6, 0.7这2个工况的相对误差稍有增大；数值方法与实验方法对于共振频率的预报一致。综合考虑，可认为本文所用数值模型在求解箱体垂荡运动响应问题上有较好的可靠性与准确性。

3.2 双箱自由垂荡运动数值结果与实验结果对比

针对固定双箱问题，文献[15]中给出了计算共振频率的详细推导过程和计算公式：

$ \omega_{g}=\sqrt{\frac{g}{B B_{g} /(h-d)+d}} $

(7)

式中：ω_g为共振频率；g为重力加速度。由式(7)结合色散公式可计算得出，本文中固定双箱发生共振时的无因次波数kh为2.27，文献[14]中通过物理模型实验得到的共振波数为kh=2.154，二者之间相对误差为5.37%。

因2个箱体均考虑了垂向的自由度，双箱共振问题变得更加复杂；不仅存在窄缝处水体共振，还涉及窄缝处水体和箱体A、B垂荡响应之间相互耦合问题。本文采用实验方法进行研究；重复性良好，实验结果取平均值。实验中的入射波频率kh包括2.00、2.20、2.30、2.35、2.40、2.45、2.50、2.60和2.80，实验结果见图 7、8。数值模拟工况与实验保持一致，入射波高均为H₀=0.024 m，在k_h在2.0~2.8内均匀取点，并在共振点附近进行加密取点来精确寻找共振频率。图 7为双箱间窄缝处的波面抬高，将本文实验及数值结果与固定双箱实验^[14]相对比分析，发现：箱体可自由垂荡运动时，窄缝处波面抬高的共振点、共振区间和共振点处的幅值均较双箱固定条件下有较大的差异；固定双箱共振点幅值要远大于自由垂荡双箱，为实验结果的1.85倍，数值结果的2.17倍；自由垂荡双箱的实验和数值结果虽存在一定的差异，但2条曲线的变化趋势一致，两者对共振频率的预报基本一致。双箱间窄缝处的水体有2种共振形式：活塞式和晃荡式。双箱间窄缝宽度较大时共振为晃荡式，宽度较小时为活塞式^[18]；本文双箱的窄缝宽度较小，水体共振形式表现为活塞式共振。

自由垂荡双箱在共振时的窄缝波面抬高均小于固定双箱情况，因波面抬高是影响箱体所受波浪力的重要因素，因此，若在多模块海洋结构物的共振问题研究中，若不考虑结构物的垂荡响应，会在错误预报共振区间及共振频率的同时，高估双箱共振现象的危害。

图 8为箱体A和箱体B在波浪砰击下的垂荡运动的RAO，数值结果与实验结果二者变化规律相同，随kh增加，二者在数值上也趋于一致，尤其kh在2.5~2.8内吻合较好。表 1给出了在差异较大的区间内，不同kh条件下实验与数值结果的比值，易知该区间内二者比值十分接近，故而认为实验中存在一个固定比例的阻尼力或能量耗散，该现象将在下一小节中进行详细分析。

图 9为箱体A和箱体B分别在水平和垂直方向的无量纲波浪力，波浪力的无量纲化为：

$ \begin{gathered} F_{x}^{\mathrm{A}}=F_{\mathrm{A} x} / \rho { gheH }_{0}, F_{y}^{\mathrm{A}}=F_{\mathrm{A}y} / \rho { gheH }_{0}, \\ F_{x}^{\mathrm{B}}=F_{\mathrm{B} x} / \rho { gheH }_{0}, F_{x}^{\mathrm{B}}=F_{\mathrm{B} x} / \rho { gheH }_{0} \end{gathered} $

式中：ρ为流体密度；e为计算域厚度；H₀为入射波高；F_Ax、F_Ay、F_Bx、F_By分别为箱体A和箱体B所受水平波浪力和垂直波浪力。

由图 9可知：箱体A在垂向和水平向所受波浪力均分别大于箱体B；箱体B所受垂直波浪力随kh的增大而趋于0，即入射波长越小，箱体B在背浪侧所受的波浪砰击影响越小；当波长增大时，箱体A和箱体B在水平方向所受波浪力趋于一定值；双箱水平波浪力和垂直波浪力的变化趋势分别一致。

为对自由垂荡运动双箱的水动力特性进一步分析，图 10给出了双箱各参数无因次化后的时历曲线，其中y_A为箱体A垂向位移，y_B为为箱体B垂向位移。由图 10可知：箱体B所受水平波浪力与窄缝处波面抬高的时历曲线相位相同，窄缝处的波面抬高是影响箱体B水平波浪力的主要因素；箱体A所受水平波浪力则由于箱体自由垂荡运动与窄缝处的波面抬高的耦合影响，窄缝处的波面抬高不再是其主要影响因素，这点与固定双箱条件下有着明显的不同^[12]；箱体A和箱体B在垂直方向所受波浪力则主要受其本身的垂向位移影响，其时历曲线相位相反。

4 误差分析

由3.2节可知，本实验与数值模拟的结果之间存在一定误差，而数值模型已经过系统的有效性与可靠性验证^[16]，因此基于数值结果对模型实验进行误差分析。因实验中对每个工况都进行了多次实验，且重复性良好，因此排除偶然结果的影响。由表 1可知，在一定波数下，实验值与数值解之比较为稳定，而箱体A和箱体B垂荡RAO要小于数值结果，窄缝处的波面抬高的实验值要大于数值结果。综上可知：共振点处双箱运动幅值越大，窄缝处的波面抬高越低，符合能量守恒原理。

数值模拟与实验结果存在偏差有诸多因素，如：导轨存在对流场的影响、导轨安装(摩擦与偏向)的影响、箱体与水池边壁间缝隙的影响等。下面对几个关键因素进行分析与讨论。

4.1 导轨安装误差的影响

安装误差见示意图 11，是夸张化的导轨安装误差示意图。如图 11所示，安装误差可能使得箱体在导轨中运动时存在较小阻尼，而且导轨与滑轮本身又存在些许的摩擦。由于水槽宽度仅为1.0 m，人员的作业环境和作业条件较差，在导轨安装过程中容易产生偏差；在导轨偏向上的任何一丝偏差，对结果会产生较大的影响。

4.2 箱体-水槽的间隙影响

在安装箱体时，为保证水槽侧壁的安全，在箱体与水槽之间留有33 mm的缝隙，如图 12所示，图 12为实验装置的俯视图。在波浪的砰击下，会有一部分的波浪从该缝隙中穿过。为探究该缝隙的存在对于实验的具体影响，再次通过数值方法1∶1模拟了3个典型工况(kh=2.20、2.40和2.60)下，窄缝存在时的双箱算例，结果如图 13、14所示。由于实验中池壁与箱体的间距较小，此时相当于把水槽沿长度方向分成了3个分水槽，3个分水槽在结构物周围有一定的相互干涉作用，但总体来说：箱体-水槽的间隙对计算结果影响较小。

4.3 导轨存在对流场的影响

导轨会引起其周边流场的剧烈变化，从而影响双箱的运动响应。本节1 ∶1还原了有导轨时的三维双箱共振算例，仍然选取典型工况kh=2.20，2.40，2.60，其结果如图 15所示。由图 15可知，导轨的存在对双箱共振特性的影响十分明显，尤其是窄缝处的波面抬高和箱体B垂荡RAO，有导轨时窄缝处波面抬高的RAO与箱体B垂荡RAO均介于实验值与数值结果之间，可见导轨的存在对流场的干涉程度较大，对垂荡双箱的共振响应影响较大。

5 结论

1) 本文建立的二维粘流数值波浪水池可正确求解波浪中单箱和双箱的垂荡运动问题，且具有良好的稳定性与准确性。

2) 在研究双箱共振问题时，如果忽略两毗邻结构的垂荡影响，假设为固定不动条件，会在预报共振区间与共振频率时出现较大偏差。

3) 考虑箱体的垂荡运动后，各箱体在窄缝共振发生时的水动力会发生变化。共振点处，箱体A的水平波浪力受到窄缝处波面抬高和箱体垂荡运动的耦合影响；箱体B的水平波浪力受窄缝处的波面抬高影响；箱体A所受力的幅值都大于箱体B。

4) 本文所用实验方法能准确预报可垂荡双箱的窄缝共振频率；有利于自由垂荡双箱在不同入射波频率下窄缝处波面抬高、箱体垂荡位移等参数的变化趋势，对于研究箱间水体的窄缝共振问题具有积极意义。