海洋立管是连接海底和海洋平台的细长柔性结构。在来流作用下，立管两侧的涡脱落将导致立管在横流和顺流2个方向产生涡激振动，当涡脱落频率和立管振动频率锁定时，涡激振动会被放大，从而加剧立管的疲劳破坏。因此，有必要对涡激振动问题开展深入的研究。

回顾过去近半个世纪的时间里，对于涡激振动的研究已经取得了诸多成果，具体可参阅文献[1-4]。弹性支撑圆柱体是研究涡激振动的典型结构，许多学者采用该结构开展了圆柱体单自由度^[5-6]和双自由度^[7-8]的涡激振动实验研究。随着计算机和计算流体力学(CFD)方法的不断发展，越来越多的学者开始关注基于CFD的涡激振动研究，该方法直接求解N-S方程，可以有效捕捉流体力和涡脱落模态等^[9-12]。

在实际海洋环境中，波浪和顶端平台运动带动立管与流体之间产生的相对振荡流会导致涡激振动的产生。Sumer等^[13]、Kozakiewicz等^[14]和王俊高等^[15]采用实验的方法对振荡流中圆柱体的涡激振动特性开展研究，发现了间歇性振动和模态跃迁等现象，推动了振荡流中涡激振动数值研究的发展。邓迪等^[16]采用切片法对王俊高的实验模型进行了数值模拟，有效捕捉到了振荡流中涡激振动的现象。邓跃^[17]选取KC数为0.56和2.51的振荡流，研究了均匀流和振荡流叠加下圆柱体的涡激振动，分析了圆柱体的振动幅值和锁定区间。Zhao等^[18]研究了振荡流KC数为10和20时圆柱体的涡激振动，发现当最大折合速度小于8时，单一频率振动，而当最大折合速度大于8时，多频率振动，且振动周期会受到涡脱数量和方向的影响。整体来看，针对振荡流中圆柱体涡激振动的研究采用的KC数较小，且缺少涡脱落模态和流体力之间关系的分析。

本文基于Fluent软件建立了振荡流中圆柱体横流涡激振动的二维数值模型。该模型通过两端边界的简谐运动实现振荡流，并采用了RANS方程和SST k-ω湍流模型模拟流场。采用该模型对已公布的实验模型进行分析，验证了数值模型的有效性。基于该模型开展了振荡流KC数为25和50这2种情况下圆柱体的横流涡激振动模拟，重点研究了圆柱体的激发模态、流体力和涡脱落特征。

1 振荡流数值方法

尽管在湍流状态下，涡激振动的涡脱落具有一定的三维特征，但是已公布的研究表明二维数值模型可以较好的模拟流体力和涡脱落模态^[9]，因此本文采用二维模型进行模拟。图 1为二维圆柱体涡激振动的计算流域，圆柱体直径为D，流向长度为80D，横向宽度为25D，圆柱体位于计算流域的中心位置。通过左右两侧边界的来回运动产生振荡流，其速度表达式为：

$ U\left( t \right) = {U_m}{\rm{cos}}(2{\rm{ \mathsf{ π} }}{f_{os}}t) $

(1)

式中：U_m为边界往返运动的最大速度；f_os为往返的频率；t<为时间。振荡流的KC数定义为${\rm{KC}} = {U_{m}}/({f_{os}}D) $。

对于振荡流中流场，本文采用非定常不可压缩的RANS控制方程模拟。在考虑圆柱体运动带来的网格变化时，RANS方程的表达式为：

$ \frac{{\partial {U_i}}}{{\partial {x_i}}} = 0 $

(2)

$ \frac{{\partial {U_i}}}{{\partial t}} + {C_j}\frac{{\partial {U_i}}}{{\partial t}} = - {\rm{ }}\frac{1}{\rho }{\rm{ }}\frac{{\partial P}}{{\partial {x_i}}} + {\rm{ }}\frac{\partial }{{\partial {x_j}}}(2\nu {S_{ij}} - \overline {{u^\prime }_i{u^\prime }_j} ) $

(3)

式中：U_i为流体速度；ρ为流体密度；x_i为空间坐标；t为时间；P为压力；ν为运动粘度；$ {C_j} = ({U_j}-{\overline U _j})$ 为流体相对于网格的速度项；U_j为网格运动速度；S_ij为平均应变张量；$ - \overline {{u^\prime }_i{u^\prime }_j} $ 为雷诺应力项即湍流效应，分别表示为：

$ \begin{array}{*{20}{c}} {{S_{ij}} = {\rm{ }}\frac{1}{2}{\rm{ }}\left( {\frac{{\partial {U_i}}}{{\partial {x_j}}} + \frac{{{\rm{ }}\partial {U_j}}}{{\partial {x_i}}}} \right)}\\ { - \overline {{u^\prime }_i{u^\prime }_j} = 2{\nu _t}{S_{ij}} - {\rm{ }}\frac{2}{3}{\delta _{ij}}k} \end{array} $

(4)

式中：v_t为湍流粘性系数；k为湍流动能；δ_ij为克罗内克尔符号。由于雷诺应力项的存在，本文结合SST k-ω 湍流模型求解方程(3)，该湍流模型广泛应用于涡激振动的研究，且预报结果理想^[10]。

圆柱体两侧的涡脱落产生的升力可以激励圆柱体横流方向的振动。圆柱体振动的控制方程可表示为:

$ m\ddot y\left( t \right) + c\ddot y\left( t \right) + ky\left( t \right) = {F_y}\left( t \right) $

(5)

式中：y是圆柱体垂直于来流方向上的位移；m、c、k分别为单位长度圆柱体的质量、结构阻尼与弹簧刚度；F_y表示作用在圆柱体上的升力，通过圆柱体表面的压力积分得到。本文采用Newmark-β方法求解方程(5)。

圆柱体和左右边界的运动通过UDF传入到Fluent，两者的运动会引起网格的变化，本文采用扩散光顺的网格光顺方法，通过求解扩散方程得到网格节点的位移，从而实现网格运动。网格节点的计算位移公式为：

$ \nabla \cdot (\gamma \nabla {\mathit{\boldsymbol{U}}_j}) = 0 $

(6)

式中γ为扩散系数，本文取值1.5。

2 验证分析

本节采用振荡流中的涡激振动实验模型^[13]对数值模型进行验证，并与Zhao等^[18]的数值结果进行比较。实验模型的相关参数为：KC=10，D=0.02 m，f_N=0.77 HZ。计算流域内设置了3个速度监测点A⁺、B⁺、C⁺，来分析两端边界的简谐振动产生的流场，如图 1所示。图 2为部分区域计算网格，圆柱体表面有200个节点，边界层厚度为0.001D，保证无量纲壁面距离y⁺小于1.0，无量纲步长τ为0.002，这里定义τ=(tU_m/D)。

图 3(a)是U_r=10，U_m=0.154 m/s时，3个速度监测点沿流向的速度随时间变化曲线，U_r为最大折合速度，定义为U_r= U_mf_N/D。可以发现A⁺、B⁺、C⁺的速度变化相同且最大值为0.154 m/s，与边界U_m相同；另外，从图 3(b)可以看出，除了圆柱体的附近区域，流场的速度均匀分布。因此，边界往返运动可以有效地产生振荡流场。

图 4是本文预报的平均幅值比(A_y/D)与Zhao等^[18]模拟的数值结果以及Sumer^[13]的实验结果的曲线图。可以看出，在U_r小于5时，本文验证结果与Zhao等^[18]的结果相近，均高于实验结果；但在U_r大于5时，本文结果介于文献[18]和实验结果之间。与实验结果相比，数值预报产生的锁定区间偏大，这不排除实验误差的原因。因此，需要后续的实验研究和更有效的三维数值模型开展进一步验证。

3 算例分析与讨论

本节对KC=25和KC=50，最大折合速度U_r从2~15的情况进行了数值模拟，并对位移、激发频率、流体力和涡脱落特征进行了分析。圆柱体直径D为0.04 m，质量比为2.2，静水中的固有频率f_N为0.98 Hz，雷诺数Re由最大流速确定，取值范围为3 185~23 208。

3.1 位移和频率

图 5是KC=25、最大折合速度U_r为5、8和15时, 圆柱体在y方向的时间位移历程曲线。可以看出在U_r=5和8时，振动较为规律；而在U_r=15时，振幅曲线上部出现2次峰值，下部幅值相对较大，同样的现象在Zhao等^[18]模拟的KC=10的情况也有出现。

图 6是KC=50、最大折合速度U_r为5、8和15时，圆柱体在y方向的时间位移历程曲线。观察到U_r=5时曲线变化与KC=25、最大折合速度U_r为8的情况相近；而在U_r=15时，2次峰值不再只出现于正y方向，而是上下交替出现。综合KC=25和50的情况来看，在U_r较小时，圆柱体在振荡流中的涡激振动比较规律。

为了分析圆柱体的振动频率，对KC=25和KC=50的圆柱体位移进行快速傅里叶变换，得到幅值谱。由图 7可以发现当KC=25、U_r=5时，圆柱体的振动频率主要固定在1倍的固有频率，还包含0.6倍和1.4倍的固有频率，当U_r=15时，圆柱体的振动频率主要固定在1.1倍、2.3倍和3.5倍的固有频率；在图 8中可以看到，当KC=50、U_r=5时，圆柱体的振动频率主要固定在1.06倍固有频率，还包含0.87倍和1.26倍的固有频率，而在U_r=15时，圆柱体的振动频率的数目明显增多。因此，振荡流中圆柱体的涡激振动属于多模态振动。另外，比较图 7和8可以看出，KC数越大多模态特性越明显。

下面主要讨论随着最大折合速度的增加，平均振幅比A_y和各激发模态的振幅及对应频率比的变化情况。每一个U_r取3个主要激发模态，振幅定义为A_y1，A_y2和A_y3，对应的频率依次为f₁、f₂和 f₃。当KC=25时，结合图 9(a)中最大激发模态的振幅和频率比变化可以发现，锁定从U_r=4开始，在U_r =10时结束；当KC=50时，从图 10(a)中可以看到，锁定区间为 U_r 从4~12。可见，随着KC数的增大，锁定区间有变大趋势。从图 9(a)和10(a)的幅值比包络线可以看出，KC数对涡激振动中的最大幅值影响较小，均在0.8倍直径附近。从图 9(b)和10(b)中频率比的变化可以发现，当处于锁定区间时，幅值最大模态的频率(主导频率)锁定在固有频率附近；锁定结束后，当KC=25时，主导频率随着最大折合速度的增大而增大，但是当KC=50时，主导频率仍然在固有频率附近，其原因可能是大KC数对应大振荡周期，所以低流速持续的时间较长，容易产生低频涡脱落，从而产生低频流体力，如图 11(b)中虚线圈出部分，而离固有频率较近的低频流体力更容易激发圆柱的振动；另外，KC数较大时涡脱落产生的高频流体力的频率间隔要小于KC数较小的情况，如图 11所示，而相近的高频流体力会相互干扰，不利于锁定。

3.2 涡脱落和流体力

圆柱体的受力和振动可以通过观察圆柱体的涡脱落和流场变化来解释，由于云图较多，现结合圆柱体在y方向的受力和位移曲线选取几组典型涡量变化等值线图进行分析，图中X=x/D，Y=y/D，横向箭头表示流向，竖直箭头表示圆柱体运动方向。

图 12和图 13是KC=25、U_r=5时，不同时刻的涡量等值线图、圆柱体的流体力和位移曲线图。从图 12(a)和12(b)中可以发现，图中涡1和涡2是上一周期脱落的漩涡，在流向改变后，分别从圆柱体的上方和下方经过，此时对应图 13中在t=81.6 s和82.08 s左右出现流体力峰值的时候。从图 12(c)和12(d)中可以发现，流速经过峰值(或谷值)后，圆柱体的涡脱模式发生明显改变，前者为流速谷值之前，后者为谷值之后，此时对应图 13中A区域受力曲线出现小幅波动，相同的情况如图中的虚线箭头。

当KC=50，U_r=5时，图 14中出现了“2P”形式的涡脱模式，但这种涡脱模式并不是对称的，如图 14(b)和14(c)所示，这种改变引起了图 15中流体力曲线在B区域出现持续的小幅波动。变化较稳定的A区域，有较大的流体力，原因同图 13一样，也与脱落的涡在流向改变后经过圆柱体有关。另外还可以看出，图 15中位移曲线的峰值略迟于力曲线的峰值，但与流速的最大值差不多同时出现，在图 13中也出现同样的现象。

在U_r=15时，KC=25振荡流的流速及对应的圆柱体受力和位移见图 16。在图 17(b)中，振荡流方向开始变化，即将脱落的涡1沿着圆柱体上表面滑动到另一侧，导致图 16中受力曲线A区域的小幅波动。KC=50时的涡量云图及流体力特征与KC=25时相似，图 18 (a)的涡量云图与图 19中受力曲线A区域的小幅波动对应。由图 16和图 19中可以发现，流体力的最大值和流速的峰值基本同时出现，与上述U_r=5的结果不同。

在U_r=15时，尽管KC=25和50对应的涡脱落方式相似，但KC=50时涡数量偏多，这是因为大KC数振荡流产生涡脱落的流速时间跨度较大。在流向反转之后，这些脱落的涡反向运动经过圆柱体，影响圆柱体的受力和振动。这里主要讨论KC=50的情况，如图 18(b)，上半周期脱落的涡2主要从圆柱体上侧“相交”而过，并被圆柱体分为2个大小不等的涡，涡3主要从圆柱体下方“相切”而过，涡4从圆柱体上方“相切”而过。这些涡经过圆柱体的时刻对应流体力的幅值1、2和3，均高于后续出现的幅值，综合之前的情况说明流向反转后经过圆柱体的涡可以增加圆柱体的流体力。由于圆柱体的受力出现1、2和3及后续的高频幅值，导致涡激振动位移曲线出现B和C 2个区域的高频振动，振动幅度较小的原因是高频力的频率远离静水中的固有频率。

4 结论

1) 振荡流能激发圆柱体的多模态振动，且KC数越大多模态特性越明显；另外，处于锁定区间时，主导模态的幅值明显高于其他模态。

2) 随着KC数的增大，振荡流中的锁定区间有变宽趋势，一个振荡流周期内从圆柱体脱落的涡数也会增加，而KC数的增大似乎对涡激振动幅值比包络线的最大值影响较小。

3) 在流速方向发生改变时，圆柱体已脱落的涡会随着流体从圆柱体的上下两侧交替经过，增加圆柱体的流体力。

4) 在大KC数和大U_r 的情况下，最大流速附近出现多个周期的高频流体力，并引起圆柱体的小幅高频振动。