亚临界雷诺数下单圆柱涡激振动研究

周力，邱中秋，袁亚帅，宗智

大连理工大学 船舶工程学院，辽宁 大连 116024

0 引 言

海洋立管结构大多以柔性管形式存在，在洋流的作用下自由振动。然而由于海底地形复杂多变，铺设在海底的某段立管两端容易受到约束而形成一定的跨度，在洋流作用下振动并产生破坏。研究中通常把这段立管作为弹性安装的刚性管考虑，利用自激振动或强迫振动的方法对其相关水动力特性进行探讨[1-2]。

Williamson 等通过强迫振动实验研究对振动立管尾涡脱落的结构模式进行了分类，并依据涡脱模式划分了不同的区域[3]。Peppa 等研究分析了低雷诺数下来流和立管结构之间的能量传递关系，发现在较低无量纲振幅比时尾涡脱落模式为“2S”（S 表示单个旋涡脱落），较高无量纲振幅比时尾涡脱落模式逐渐变得复杂[4]。Meneghini 等将立管的振动频率和幅值作为变量，对Re＝200 时的横向振动立管尾涡脱落频率与“锁定”现象进行了研究，确定了该雷诺数下发生“锁定”时的最大振幅值及振动频率的变化范围[5]。Morse 等[6-7]通过圆柱的强迫振动和自激振动组合实验发现了新的尾涡脱落模式，位于“2S”和“2P”模式的过渡区（P 表示一对旋转方向相反的旋涡脱落），并将其定义为“2P0”模式。

王凯鹏等[8]基于紧致插值法对均匀流(Re=200)中圆柱横向强迫振动问题进行了二维数值模拟研究。朱永健等[9]通过对定常流中横向振动圆柱的升力突变现象研究，发现随着强迫振动振幅的增加，圆柱脉动升力系数出现了突变衰减，圆柱静止情况下的涡脱频率与受迫振动频率控制的波动出现了相位逆变，使得圆柱与流体的能量传递出现了逆向改变。

从海洋工程应用角度来看，弹性支撑立管强迫振动研究的主要目的是揭示自激振动更深层次的规律，为海洋立管涡激振动预测模型提供可靠的实验数据[10]。事实上，大多数涡激振动（vortexinduced vibration，VIV）预测程序是基于强迫振动实验数据形成数据库后建立相应的预报模型，以满足工程应用需要，如VIVANA 和SHEAR7[11-12]。

自激振动能更直观地观察圆柱受到流体力作用发生VIV 时的现象，但强迫振动能更深层次地揭示流体力与振动相互作用下的水动力特性[13]。强迫振动和自激振动的流场变化及振动响应存在很大的差异，目前鲜有文献对2 类实验的内部联系进行深入研究。基于此，本文将结合强迫振动数值计算与自激振动实验2 种方法，将强迫振动的动力响应、涡脱模式的转变与该圆柱自激振动出现“锁定”时的最大振幅联系起来，期望为深入讨论二者之间的联系提供参考。

1 数值模型

1.1 控制方程

以二维不可压缩黏性流体模拟流场运动，其控制方程Navier-Stokes 方程组在坐标系下表示为：

本文采用的求解器为 STAR-CCM+ 软件中的Realizablek-ε 湍流模型。标准k-ε 湍流模型的主要优点是计算快速、稳定，计算结果合理，适用于高雷诺数的流动，但不建议用于分离流动。Realizablek-ε 湍流模型是标准k-ε 湍流模型的一个改进模型，可以较好地模拟带有强压力梯度的边界层流动和分离流问题，由Shih 等[14]提出，其流湍动能及其耗散率输运方程为：

1.2 计算域和边界条件

STAR-CCM+ 软件可以有效地实现圆柱在均匀来流下发生横向强迫振动的模拟。首先对圆柱作减运算，对包含减运算在内的重叠网格施加不同的场函数，可实现模拟圆柱在不同频率和振幅下的运动。流经圆柱的均匀流速度为U0。

计算域如图1 所示，来流速度方向与x轴平行，圆柱横向振动方向与y轴平行，圆柱直径D=0.05 m，圆柱中心上下两侧20D处为对称边界；左端距离圆柱中心20D处为速度进口，右端距离圆柱中心50D处为压力出口；圆柱表面为无滑移壁面条件。

图1 圆柱受迫振动计算域示意图Fig. 1 Schematics of the forced-oscillation computational domain

圆柱垂直于来流方向沿y轴做横向强迫振动，瞬时位移：

式中：A为圆柱横向振动振幅；fe为圆柱强迫振动的激振频率。

圆柱强迫振动过程中沿y轴方向的瞬时速度：

圆柱开始振动时的瞬时速度vy(0)=2πfeA。

2 计算可行性验证

2.1 固定圆柱验证计算

网格划分及其数量会直接影响到计算结果的准确性。为此，采用多边形和棱柱层的网格划分方法，对圆柱横向振动的延展区域和尾部区域进行了局部加密，在圆柱周围区域设置包裹圆柱运动的运动重叠网格，重叠区域大小为2D×5D。网格划分如图2 所示，计算模型是具有较高精度的结构化网格。采用壁面函数法进行圆柱体近壁面处理，圆柱体表面第1 层网格满足y+≈1 条件。由Re=20 000，计算得到来流速度U0=0.356 m/s。流体的密度取为997 kg/m3，其动力黏度取为0.895×10−3kg/(s·m)，则可计算得到紧邻壁面第1 层网格的厚度为∆s=0.04 mm。设置棱柱层数n=10，增长率q=1.3，由δ=(∆s(1−qn))/(1−q)，则可计算得到边界层总厚度为1.67 mm。

图2 固定网格划分示意图Fig. 2 Schematics of fixed mesh division

对网格、时间步长进行了无关性验证，在保证计算精度的同时，应节约计算成本，最终选取网格数为39 万，时间步长为0.01 s。为了验证计算模型的可靠性，对雷诺数Re=20 000 均匀流中固定圆柱和振动圆柱的绕流问题进行了数值计算，如图3 所示，通过观察升力系数的能量谱密度曲线，确定圆柱涡街脱落的特征频率为1.46 Hz。

图3 Re=20 000 固定圆柱升力系数能量谱密度Fig. 3 Energy spectral density of lift coefficient of fixed cylinderwith Re=20 000

根据公式：

其中，fov为圆柱尾涡脱落频率，可得到该固定圆柱在Re=20 000 时的斯特劳哈尔数St=0.204 8，与文献[15]中St=0.2 基本保持一致。

2.2 圆柱强迫振动计算

在相同来流条件下对圆柱进行强迫振动数值计算，将振幅比A/D=0.5，强迫振动频率比范围0.5≤fe/fs≤1.5（fs为斯特劳哈尔频率）时的数值计算结果与喻晨欣等[16]的数值结果和Sarpkaya[17]的实验结果进行对比验证。从图4 可以看出，本文计算结果与文献[16-17]的最大升力系数幅值变化趋势保持一致，随着频率比的增大，均呈现不断增大的趋势，并在fe/fs=1 附近出现了跳跃[18-19]。由于本文与参考工况对应的雷诺数存在差异，对应的频率比也有所不同，导致结果存在偏差，与实际情况相符。以上算例验证了本文模型可用于求解亚临界雷诺数下圆柱的强迫振动问题。

图4 A/D=0.5 时最大升力系数Clmax 随fe/fs 变化曲线Fig. 4 Clmax variation versus fe/fs with A/D=0.5

3 计算结果分析

3.1 强迫振动数值结果分析

本文采用外加激励的方法，使圆柱在均匀来流下做类似自激振动的运动。通过控制振幅A和激振频率fe，可以实现使圆柱位移时域曲线呈简谐运动的规律，无量纲振幅比A/D和激振频率比fe/fn如表1 所示（fn为结构固有频率，考虑附加水质量）。涡激振动达到“锁定”时对应的振幅存在阈值[3]，最小振幅比取A/D=0.3。通过自由衰减振动实验测得fn=0.714 2 Hz。

表1 无量纲频率比及振幅比范围Table 1 Non-dimensional frequency ratio and amplitude ratio

考虑到强迫振动过程中圆柱受到流体力和激振力的双重作用，假定在数值计算中得到的升力为Fm，完全由流体力提供的升力为Fl，强迫振动对圆柱升力的影响用惯性力Fi表示，三者关系可用式（9）来表示：

当流体力Fl与惯性力Fi方向相同时，Fm>Fl；当Fl与Fi方向相反时，Fm

图5 为频率比fe/fn=0.83 时不同振幅比下，单周期内圆柱振动速度vy曲线及相应升力系数Cl曲线。从图中可以看到，振动速度曲线始终遵循简谐运动规律，但由于受到强迫振动带来的惯性力影响，单周期内升力系数曲线不再规则变化。当振幅比A/D足够大时，可观察到升力系数曲线单周期内包含多个峰值。

图5 单周期内振动速度及升力系数变化曲线Fig. 5 Oscillating speed and lift coefficient varies within a single period

同时可以发现，随着振幅比A/D的增大，在圆柱上升阶段（0～0.5T），最大振动速度（0.25T）对应的升力系数值会由正值转变为负值；而在圆柱下降阶段（0.5T～T），最大振动速度（0.75T）对应的升力系数值由负值转变为正值。强迫振动研究中通常取圆柱固定（A/D=0）时的升力系数为正值，所以在接近A/D=0 一端，升力系数应取正值。在本文圆柱强迫振动数值计算中，将在上升阶段圆柱振动速度达到最大时（0.25T）的位置选为升力系数监测位置。

图6 为VIV 预测程序SHEAR7 中的保守模型，各无量纲频率比下的Cl-A/D曲线通常由强迫振动实验所得数据点（0，Cl0），（（A/D）*，Clmax）、（（A/D）max，0）通过二次函数拟合得到，其中（A/D）*和（A/D）max分别表示升力系数为零和最大值时对应的无量纲振幅比[19]。该模型呈现升力系数Cl随着振幅比A/D的增大先上升后下降、由正升力系数逐渐转变为负升力系数的趋势。

图6 升力系数随无量纲振幅比变化曲线Fig. 6 Variation of lift coefficient with dimensionless A/D

流体力具有增大圆柱振幅趋势的作用，而振动则会限制圆柱的振幅继续增大[20]。在两者的相互作用下，当升力系数为正时，说明振动对流体力的限制能力有限，此时若继续增大来流速度，立管振幅仍有增大趋势，因此涡激振动未达到“锁定”状态时圆柱振幅会随着来流速度的增大而增大。随着振幅比的增大，振动对流体力的限制能力也进一步增强，直至平衡位置处的升力系数为零，二者达到平衡，振幅达到最大且不再有增大的趋势。值得一提的是，在涡激振动中，振动的限制作用不会超过流体力的激励作用，即升力系数不可能取负值，所以达到“锁定”状态时，在一定流速范围内，即使增加流速，圆柱的振幅仍能保持最大值不变。

圆柱在外加激励作用下发生强迫振动，通过改变振幅比A/D，观察监测位置处的升力系数，可拟合出各激振频率比fe/fn下的升力系数随振幅比A/D的变化曲线。对各工况下数值计算所得结果进行三阶多项式拟合，并控制拟合条件Cl0=0.168（Re=20 000 时固定圆柱的升力系数），可得拟合曲线。在VIV 预测程序SHEAR7 中输入相关参数条件，同样可得预报的拟合曲线。将同一频率比下的2 种拟合曲线进行对比分析，结果如图7所示。

由图7 可见，各激振频率下数值计算结果拟合曲线与SHEAR7 拟合曲线的变化趋势高度吻合，均呈现出升力系数随着振幅比A/D增大先上升后下降的变化规律。传统强迫振动计算中升力系数的取值方法往往未呈现该规律，说明选取上升阶段最大振动速度时的升力系数，对于建立VIV 预测模型是适合的。

为了进一步了解数值计算结果拟合曲线与涡激振动振幅响应之间的联系，从频率比fe/fn=1 附近的零升力系数点着手分析。如图8 所示，可观察到圆柱平衡位置处的零升力系数点均保持在振幅比A/D=0.8 附近。通过类比SHEAR7 预测模型进行分析，可知在亚临界雷诺数下该圆柱涡激振动出现“锁定”状态时最大振幅应保持在0.8D左右。由此可对该圆柱在自激振动频率比fov/fn≈1时的最大响应振幅做出预测，其中fov为自激振动时圆柱的尾涡脱落频率。

图8 各工况下Cl 随振幅比A/D 变化曲线Fig. 8 Variation of lift coefficient with A/D under all conditions

同时由图8 可见，在同一振幅比下，随着频率比的增大，升力系数逐渐增大。当频率比fe/fn=1.42 时，相对于其他组，升力系数大幅升高，发生了突增现象[12]，相关机理有待进一步探讨。

3.2 自激振动实验对比分析

强迫振动数值计算结果表明，监测位置处升力系数为零时对应的振幅比均保持在A/D=0.8 附近。在自激振动工况下，预测弹性安装时该圆柱在亚临界雷诺数下发生涡激振动的最大振幅为0.8D左右，且此时处于“锁定”状态（自激振动频率比fov/fn≈1）。为了验证预测结果，设计并进行了圆柱自激振动实验。

图9 为自激振动实验装置图，振动圆柱嵌合在振动框架中，并弹性安装于固定装置。弹簧具有一定的的弹簧系数，振动框架可带动圆柱沿两侧的滑轨做垂向往复运动。为了尽量减小三维效应，固定圆柱两端的U 型架端口尽可能采用厚度较薄、质量较轻的端板。实验所采用圆柱的直径与数值计算中保持一致，弹簧系数K=230.48 N/m，均匀来流速度的取值范围为0.13～0.34 m/s，对应的约化速度Ur取值范围为3.7～9.4 m/s。

图9 圆柱涡激振动实验装置Fig. 9 Experimental device setup of cylindrical vortex-induced vibration

通过改变来流速度，监测得到圆柱振动幅度及相应的自激振动频率比fov/fn，可观察到振动圆柱发生“锁定”时的最大振幅，以及对应的频率比范围。图10 为自激振动实验中振幅比A/D和频率比fov/fn随约化速度Ur的变化情况。由图可见，在fov/fn接近1 时，振幅比大幅升高；且在一定的fov/fn范围内，振幅比A/D不发生大的改变，说明发生了“锁定”现象。同时可观察到“锁定”范围对应的约化速度使圆柱的最大振幅达到0.8D左右。

图10 圆柱振幅比及自激振动频率比随约化速度变化图Fig. 10 Changes in the amplitude ratio and self-induced vibration frequency ratio of the cylinder at different reduced speeds

3.3 涡脱模式转变分析

Williamson 等[3]对振动圆柱尾流中漩涡脱落模式进行了实验研究，并按尾涡脱落时的结构形式划分了不同的区域，发现“锁定”区域附近的主要涡脱落模式为2S，2P 和P+S。随着无量纲振幅比的增加，在涡脱模式转变临界线附近，通常会从2S 模式转变为P+S 或2P 模式。实验研究表明，在低雷诺数下（Re<300），P+S 模式会取代2P模式，以至于整个区域只出现P+S 模式，而未出现2P 模式。而在较大雷诺数时，P+S 和2P 模式之间的边界与临界线基本一致。可见，随着雷诺数的改变，同一区域内涡脱模式发生转变后不确定呈现2P 模式或P+S 模式，但转变临界线基本保持不变。

在同一强迫振动频率比下，选取振幅比范围为A/D=0.6～1.0 内同一时刻的涡量图。如图11 所示，各频率比工况中的振幅比由A/D=0.6 至A/D=1.0 依次增大。可以发现，在振幅比小于0.8 时，如A/D=0.7 时红色线圈所标记涡对，圆柱尾涡脱落模式为P+S 模式；当振幅比大于0.8 时，如A/D=0.9 时红色线圈所标记涡对，圆柱尾涡脱落模式为2P 模式。说明在圆柱强迫振动中A/D=0.8 附近是尾涡脱落模式发生转变的临界线，监测位置零升力系数点与尾涡脱落模式转变同时发生。

图11 各频率比下涡量随振幅比变化图Fig. 11 Variation of vorticity with amplitude ratio at various frequency ratios

4 结 论

本文研究了圆柱强迫振动时升力系数响应值的选取位置，并进行相关机理分析，对圆柱发生涡激振动“锁定”时的最大振幅进行了预报和对比验证，对平衡位置升力系数为零时的涡脱模式转变进行了探讨。通过分析，得到如下结论：

1） 亚临界雷诺数下圆柱强迫振动时的升力系数响应与自激振动相比存在较大差异，当强迫振动中圆柱速度达到最大时，监测到的升力仅由流体力提供，将圆柱在上升阶段的振动速度达到最大时（0.25T）对应的升力系数作为计算结果是合理的。

2） 数值模拟计算结果拟合曲线与VIV 预测程序SHEAR7 的相关拟合曲线吻合较好，且变化规律保持高度一致，升力系数均随着振幅比的增大，呈现出先增大后减小的趋势，提供了一种以强迫振动数值计算结果建立Cl-A/D模型用于涡激振动响应振幅预测的方法。

3） 数值计算结果表明，圆柱在强迫振动过程中平衡位置处的零升力系数点均保持在振幅比A/D=0.8 附近，采用自激振动实验方法对预测结果进行了对比验证。实验结果表明，当自激振动频率比fov/fn接近于1 时，振幅比大幅升高，且在一定的自激振动频率比fov/fn范围内，振幅比A/D没有较大改变，说明发生了“锁定”现象。同时确定了圆柱在“锁定”时对应的最大振幅为0.8D左右，证明了该预报方法具有可行性。

4） 对平衡位置处的零升力系数与相应的圆柱尾涡脱落模式进行了探讨。计算结果表明：在振幅比小于0.8 时，圆柱尾涡脱落模式为P+S 模式；当振幅比大于0.8 时，圆柱尾涡脱落模式为2P 模式。强迫振动圆柱Cl-A/D拟合曲线零升力系数对应的振幅比与圆柱在涡激振动中“锁定”时的最大响应振幅比基本一致，且圆柱尾涡脱落模式在此振幅比下发生了转变。

未来可通过对不同雷诺数工况下的圆柱振动进行计算，对圆柱在均匀来流下发生“锁定”的频率比区间进行预报。