0 引　言

深海平台是世界上普遍使用的深海油气田开采装备，这种海洋平其整体尺寸大，服役海域深。将海洋平台置于一定来流中，会在平台两侧交替的产生脱离结构物表面的漩涡，在周围流场压力作用下，平台容易发生涡激运动(VIM)现象^[1]。1988年，Williamson^[2]对振荡柱体涡旋泄放结构进行了较为完整全面的实验研究，为之后的涡激振动和涡激运动的尾涡泄放形式研究提供了依据和参考。白治宁等^[3]将CNOOC设计的新型半潜式平台作为研究对象，采用模型试验方法验证动网格技术模拟平台涡激运动的可靠性，得到了较为理想的结果。Rijken^[4]在现场和模型试验中观察到半潜式平台涡激运动(VIM)现象，使用CFD分析获得不同位置数据，显示大部分激励发生的位置。谷家扬等^[5]采用雷诺平均法对纳维叶－斯托克斯方程进行求解，将湍流模型和CFD动网格技术结合起来，研究分析了FDPSO平台的涡激运动流场特性。

本文利用Fluent仿真软件中的重叠网格技术，结合涡激运动特征方程求解子程序模拟来流速度在0.02～0.20 m/s（每隔0.02取一个值）共计10种工况下Spar平台的涡激运动情况，捕捉到了涡激运动情况下的流场特征，研究不同约化速度下，影响平台的关键参数和涡激运动过程中发生的锁定现象。

1 理论基础及计算模型 1.1 涡激运动特征方程

实验研究发现，VIM的运动方向主要为流向和横向，流向表现为纵荡运动，而横向为横荡运动。将平台涡激运动响应简化为弹簧-阻尼-质量系统，如图1所示。

纵荡和横荡方向上的运动方程可表示为：

$ m\frac{{{{\rm{d}}^2}x}}{{{\rm{d}}{t^2}}} + {C_x}\frac{{{\rm{d}}x}}{{{\rm{d}}t}} + + {K_x}x = {F_d}(t) \text{，}$

(1)

$ {\text{m}}\frac{{{{\text{d}}^2}y}}{{{\text{d}}{t^2}}} + {C_y}\frac{{{\text{d}}y}}{{{\text{d}}t}} + + {K_y}y = {F_l}(t) \text{。}$

(2)

式中： $ t $ 为时间，s； $ m $ 为柱体质量，kg； $ x，y $ 为柱体在纵荡和横荡方向上的瞬时位移，m； $ {F_l}(t) $ 和 $ {F_d}(t) $ 为结构物受到的随时间变化的升力和阻力，N； $ {K_x} $ 和 $ {K_y} $ 是系泊系统刚度； ${C_x} = 4 {\text{π}} {\xi _x}{f_{nx}}$ ， ${C_{\text{y}}} = 4{\text{π}} {\xi _y}{f_{ny}}$ ，分别为纵荡和横荡方向运动的阻尼系数。其中 $ {f_{nx}} $ 和 $ {f_{ny}} $ 为平台纵荡和横荡下的固有周期，s； $ {\zeta _x} $ 和 $ {\zeta _y} $ 分别为平台在纵荡和横荡方向的阻尼比。

DNV-RF-F103 规定，由于顺流向振动幅值的最大值只有横向振动幅值的20%，所以可以忽略平台纵荡运动的影响，仅考虑平台的横荡运动^[6]。利用四阶Runge-Kutta方法对式(2)进行求解，求出每一个时间内流场力和平台的运动位移，经化简后求解，可得横荡涡激运动微分方程的表达式为：

$ \left\{ \begin{gathered} {y_{n + 1}} = {y_n} + hy_n' + \frac{h}{6}({m_1} + {m_2} + {m_3})\text{，} \\ y_{n + 1}' = y_n' + \frac{1}{6}({m_1} + 2{m_2} + 2{m_3} + {m_4})\text{，} \\ \end{gathered} \right. $

(3)

$ \left\{ \begin{gathered} {m_1} = {h_g}({t_n},{y_n},y_n')\text{，} \\ {m_2} = hg\left({t_n} + \frac{h}{2},{y_n} + \frac{1}{2}hy_n',y_n' + \frac{{{m_1}}}{2}\right) \text{，} \\ {m_3} = hg\left({t_n} + \frac{h}{2},{y_n} + \frac{1}{2}hy_n' + \frac{1}{4}h{m_1},y_n' + \frac{{{m_2}}}{2}\right)\text{，} \\ {m_4} = hg\left({t_n} + \frac{h}{2},{y_n} + hy_n' + \frac{1}{2}h{m_2},y_n' + \frac{{{m_3}}}{2}\right)\text{，} \\ n = 0,1,2, \cdots \text{。} \\ \end{gathered} \right. $

(4)

式中：h为时间步长；m₁，m₂……为函数 ${y^{{'}}}(t)$ 在 $ {t_n} $ 点处的近似斜率； $ {y_n} $ 和 ${\text{y}}_n'$ 分别表示第n次迭代中的瞬时位移和速度。根据式(3)编写平台涡激运动预报子程序，仿真软件在求解流体与平台相互作用的流固耦合过程中，先利用自带的求解器计算平台的流体力，并将计算好的升阻力系数作为已知量调入编写好的涡激运动程序中计算平台的运动响应位移幅值。

1.2 计算模型

Spar平台的主体部分（硬舱）为大直径的圆柱体，垂立在水中，在海流作用下，柱体两侧会发生旋涡脱落，引起涡激运动。Koopman^[7]和Grifin^[8]研究了柱体结构物在振荡流中的响应特征，结果发现当圆柱体的振动频率接近于其自然频率时，圆柱体的涡激运动会抑制尾部涡旋脱落的三维效应，涡旋发放还是二维的，所以采用二维模型足以得出准确的结果。首先，在平台原有尺寸的基础上按缩尺比1:100建立Spar平台及矩形流场区域的二维有限元缩尺模型，平台设计参数参见表1。其中横荡运动固有周期通过静水衰减试验得出，模型质量比m^*≈1。

图2为整个计算域模型简图及网格划分模型，将整个流场区域设置为15 D×40 D的矩形，边界条件为：左边为速度入口，右边界为压力出口，上下为Symmetry边界条件，平台简化为D=0.368 m的圆柱模型，采用无滑移壁面边界条件^[9]。

表2给出了10种计算工况，来流速度范围为0.02~0.2 m/s，经过表中的公式计算出对应的雷诺数Re和约化速度 $ {U_r} $ 的范围分别为7360~73600和1.11~11.1。后续结果只对 $ {U_r} $ 分别为1.11，3.34，7.8，10.2时的平台升阻力变化情况进行具体分析。

1.3 平台升力F_L和阻力F_D

在涡旋的交替泄放的作用下，结构物两侧产生了脉动升力 $ {F_L} $ 和阻力 $ {F_D} $ 。脉动升力的方向与来流速度方向垂直，脉动阻力的方向与来流速度方向平行，二者是平台产生涡激运动的主要激励力。柱体上端所受的压力要小于下端，垂直于来流速度方向上的升力由此产生，由此涡旋上下交替泄放，脉动升力表现为交变形式， $ {F_L} $ 是Spar平台涡激横荡运动的主要激励力，脉动升力原理如图3所示。

在对平台涡激运动的激励力进行研究时，一般用无量纲量来表示，其中升力 $ {F_L} $ 无量纲化之后为升力系数 $ {C_L} $ ，相应地， $ {F_D} $ 经过无量纲化之后为阻力系数 $ {C_D} $ 。涡旋泄放的过程中，升阻力具有脉动特性，考虑将升阻力系数分成时均和脉动2种情况，时均下升力系数 $ {C_L} $ 一般为0，用脉动升力的均方根来计算得出脉动下的升力系数，时均阻力系数，脉动阻力系数^[10]如表3所示。

1.4 网格无关性验证

取来流速度0.1 m/s，Re=36800情况，采用SIMPLIC求解器对网格数量分别为Mesh1，Mesh2，Mesh3，Mesh4和Mesh5等5种情况下平台的绕流状态进行有限元仿真，并将计算结果依据施特鲁哈尔数St相互比较，结果如表4所示。可知，当网格数增加到Mesh3(95360)时，其升阻力及斯特朗哈尔数几乎不再变化，为节省计算资源，选取Mesh3网格数作为本次数值模拟的网格数。

2 计算结果分析 2.1 Spar平台尾流场特征

图4为约化速度U_r为1.1～11.1时，平台发生涡激运动时某一时刻的尾流场。从图4(a)中可以看出，涡旋泄放的模式为2S模式，与静止绕流状态下的形式相差不大，涡旋都呈上下相互交替泄放的形态；图4(b)中平台后方的涡泄形式与图4(a)中的形式发生了本质的变化，涡旋出现在平台中心线上下对称位置，形态由图4(a)中的“单排”转变为到“双排”，这主要是由于平台横荡运动幅值增加所导致的。该涡旋泄放的模式为2 P模式，随着约化速度的逐渐增大，尾流场后方出现了一些小型的碎涡，且碎涡的出现位置不固定，说明了平台涡激运动具有随机性；在平台后方的湍流流动得到充分发展后，湍流强度慢慢增加，尾流区域的涡旋形状由开始的规则形状被破坏，变得小且不规则。

2.2 平台静止情况下 $ {C_D} $ 、 $ {C_L} $ 变化情况

图5给出了平台在不同约化速度下升阻力系数时历变化曲线和升力系数经过傅里叶变换之后得到的FFT图谱。从升阻力系数时历曲线图可以看出，升阻力系数均呈现出周期性变化的趋势， $ {C_D} $ 周期大约为脉动 $ {C_L} $ 周期的一半，来流速度的增加使得旋涡脱落频率增加，升阻力系数的变化周期也相应缩短； $ {C_L} $ 均值趋近于0，这是因为流场流过平台产生的涡旋具有对称性且上下涡旋脱落频率大致相同所导致的结果。FFT分析图谱中给出了脉动升力的频率分布图， $ {C_L} $ 频率为单一的窄带状，在平台静止绕流中，涡旋泄放的频率与脉动升力的频率一致，且随着约化速度的增加，平台最大升力时对应的频率也在逐渐增大。

图6(a)～图6(c)分别给出了不同来流速度下最大升力系数幅值和阻力系数均值曲线以及涡旋泄放频率图，可以看出本文所采用的计算模型得出来的结果与文献所计算的差别不大。在约化速度 $ {U_C} $ 为0.02~0.2 m/s范围内，随着 $ {U_C} $ 的增加，升力系数幅值 $ {C_L} $ 和阻力系数 $ {C_D} $ 均值相应的减小，不过在 $ {U_C} $ =0.1 m/s时， $ {C_L} $ 和 $ {C_D} $ 减小的速度下降，不像之前速度变化快，而涡旋泄放频率 $ {f_{st}} $ 与 $ {U_C} $ 成正比，呈现上升的趋势。为了更加直观展示出 $ {C_L} $ 与 $ {C_D} $ 的关系，图6(d)给出了来流速度在0.1m/s下的 $ {C_L} $ 与 $ {C_D} $ 的轨迹线，可以看其轨迹线呈歪“8”字形。这也从流体激励的方面解释了涡激运动中横荡与纵荡的轨迹线呈“8”字形的根本原因。

2.3 $ {C_D} $ 和 $ {C_L} $ 及A变化情况

图7(a)～图7(d)为不同约化速度下，平台发生涡激运动时的升阻力系数随时间变化的曲线和无量纲振幅A及其相对应的FFT分析（最大无量纲振幅A的值取涡激运动达到稳定之后的值）。在约化速度为1.11～10.2范围内，不管是哪一个约化速度下，无量纲振幅比A的值都跟升力系数的幅值紧密相关。在达到一定时间之后，升阻力系数都会呈现出周期性的变化规律，升力系数变化的幅度较大。随着约化速度的增大，涡激运动下的升阻力系数变化规律与静止绕流状态下的变化规律一致，都是出现减小的趋势， $ {U_r} $ 从3.34开始， $ {C_L} $ 幅值相较于静止绕流状态出现大幅度减小的现象，且无量纲振幅比A的值时而增大时而减小。这表明，平台由静止绕流状态进入涡激运动状态，平台受到涡激运动的影响，横向升力幅值减少，平台涡激运动开始进入锁定状态。当 $ {U_r} $ >3.34时，频率图谱中出现了2个能量比较集中的频率带：一个为平台的振荡频率，另一个为平台的自然频率；当 $ {U_r} $ =7.8时，A频率的2个极值最明显，能量最为集中，平台振动幅度也相对较大；当 $ {U_r} $ >7.8时，无量纲振幅比A急剧减小，说明平台脱离锁定区。

2.4 频率锁定现象

在涡泄频率接近平台固有频率时，平台涡激运动的频率比( $ {f_{vim}}/{f_n} $ )被锁定在约等于1的位置，在一定的来流速度范围内，涡泄频率不再随着 $ {U_r} $ 的增加而增加，这种现象被称为频率锁定。在均匀流场中的柱体发生涡激振动或者涡激运动的状态可分为锁定状态、过渡状态以及非锁定状态3种。其中，过渡状态又可根据与锁定状态的关系分为上阶段和下阶段，从非锁定状态进入到锁定状态这一阶段称为下过渡阶段，经历过锁定状态到离开锁定状态这一范围称为上过渡阶段^[11]。图8和图9给出了利用重叠网格技术计算出来的不同约化速度下的频率比和无量纲振幅比A的最大值。通过与文献[12]模型试验的结果比较，大部分数据点的误差都在10%之内，间接证明了利用Overset Mesh技术对涡激运动中流场信息进行更新的可行性和准确性。图8中， $ {f_{st}} $ 为平台静止状态下的涡泄频率， $ {f_{vim}} $ 为平台在涡激运动下的涡泄频率， $ {f_n} $ 为平台的自然频率。

从图8可以看出，平台静止绕流状态下的频率比与约化速度之间呈现出线性递增关系，并且涡泄频率 $ {f_{st}} $ 随着约化速度的增加而增大。结合图8和图9分析得出：当平台发生涡激运动且约化速度1.11< $ {U_r} $ <3.34时，下过渡阶段平台的运动幅值非常小，大约为0.05D；当 $ {U_r} $ 达到3.34之后，涡泄频率 $ {f_{vim}} $ 与平台自然频率 $ {f_n} $ 相接近，频率比大约为1，锁定现象开始发生，无量纲振幅比A大幅度增加；当 $ {U_r} $ =6.68时，无量纲振幅比A达到了最大值，为0.7D左右；当 $ {U_r} $ =3.34～7.8时，锁定区域内的涡泄频率 $ {f_{vim}} $ 不再随着约化速度 $ {U_r} $ 的增加而增加，而是被锁定在平台的自然频率上；当7.8< $ {U_r} $ <11.1时，上过渡阶段的涡泄频率开始高于平台的自然频率，平台涡激运动脱离锁定区域，无量纲振幅比A也急剧减小。所以在进行平台结构设计时，要充分了解平台所处的海洋环境，避免或者减小平台在所处环境海流作用下的频率锁定概率，减小其振动幅度，保证平台安全性。

3 结　语

本文利用四阶Runge-Kutta方程求解Spar平台涡激运动中横荡运动微分方程，根据表达式编写涡激运动子程序，结合重叠网格技术对涡激运动流场进行更新。通过数值模拟方法捕捉到平台在来流速度 $ {U_r} $ =1.11～11.1范围内的流场特征，得到以下结论：

1) 通过与文献[12]模型试验的结果相比较，采用Overset重叠网格技术计算得到的无量纲振幅A跟频率比的误差都控制在10%之内，捕捉到的涡旋泄放图像较为清晰。表明在模拟平台涡激运动时，重叠网格技术完全可以用于模型的网格划分，并且能够克服动网格重构问题；

2) 平台涡激运动历经一定时间之后，阻力系数 $ {C_D} $ 和升力系数 $ {C_L} $ 呈现出周期性变化的规律，并趋近于稳定。随着约化速度的增加，无量纲振幅比A的值时而增大时而减小。当 $ {U_r} $ 大于3.34时，频率图谱中出现了2个能量比较集中的频率带；当 $ {U_r} $ =7.8时，A的频率2个极值最明显，能量最为集中。在平台设计时，要考虑外部环境可能引起的涡激运动，通过改变平台的结构参数减少涡激运动的影响；

3) 涡旋在Spar平台的硬舱两侧呈现上下相互交替泄放的形态，且关于平台中心线对称。随着约化速度的增加，平台后方的旋涡泄放由2 S模式逐渐演变为2 P模式，形态由“单排”转变为“双排”，尾流场后方出现了一些小型的碎涡，碎涡的出现位置不固定，说明了平台的涡激运动具有一定的随机性。