柱形结构物的涡激振动(或涡激振动)是一种在工程领域经常发生的现象，如海底管线、海洋立管、桥梁、高耸建筑等。因为出现这种现象会减少结构物的寿命，容易造成疲劳破坏，所以涡激振动问题已引起了工程界和学术界的广泛关注。

由于流动问题的复杂性，关于涡激振动的研究主要集中于实验和数值模拟，而其中又以对单柱体的研究为主。Feng^[1]的研究最早发现了柱体振幅随折合速度变化的两种模态，Williamson等^[2-4]则在水槽中对小质量比圆柱进行了一系列实验。Anagnostopoulos等^[5]对层流状态下的涡激振动进行了研究，雷诺数设置为90≤Re≤150，对数值模拟的结果具有非常大的参考意义。近期，随着计算机技术的发展，大量数值模拟方法被开发出来, 从不同的角度，研究了阻尼比、质量比、雷诺数、折合速度、计算域大小等参数对计算结果的影响^[6-8]。

多柱体涡激振动中，由于多个柱体同时与流体产生相互作用，其振动和流场状态要比单柱涡激振动复杂得多。Carmo等^[9]对串列双圆柱(上游圆柱固定)的涡激振动响应进行数值模拟发现串列情况下的锁定区间比单柱更宽。Prasanth等^[10]对串列和错列双圆柱进行了研究发现错列排布情况下流向振动更为剧烈而横向振动结果与串列类似。Borazjani等^[11]在串列双圆柱的涡激振动中发现“间隙流”现象，即上游圆柱脱落的涡从两柱间隙通过。陈文曲等^[12-14]发现双自由度下圆柱的振幅和出现振幅峰值的频率都比单自由度时要大。实验研究方面，Sanaati等^[15]针对串列双圆柱涡激振动的实验表明，与单圆柱相比，上游圆柱在小间距比下不会出现“上端分支”，而在大间距比下，上下游圆柱的振动响应均接近于单圆柱的情况。Huera-Huarte等^[16]则对并列双圆柱进行了研究，并针对一柱固定一柱自由振动和两柱均自由振动的情况分别进行了实验。

本文采用了自主研制的CIP-ZJU模型对雷诺数Re=150时串列双圆柱的涡激振动问题展开了数值模拟，该模型已被应用于固定柱体绕流^[17-18]和柱体强迫振动问题^[19]的领域，应用于涡激振动问题尚数首次。文中取两种情况进行研究，即上游圆柱固定和上游圆柱自由振动，讨论两者之间的区别。考虑到间距比对串列双柱的影响很大，本文的间距比范围取为2~8。

1 数值计算模型介绍

涡激振动问题是一种流固耦合问题，需考虑流场、结构物和耦合三个过程。流场模型考虑二维粘性不可压缩流体，控制方程为质量守恒方程和N-S方程，其张量形式如下

$ \frac{{\partial {\mathit{\boldsymbol{u}}_i}}}{{\partial {\mathit{\boldsymbol{x}}_i}}} = 0 $

(1)

$ \frac{{\partial {\mathit{\boldsymbol{u}}_i}}}{{\partial t}} + {u_j}\frac{{\partial {\mathit{\boldsymbol{u}}_i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} =-\frac{1}{\rho }\frac{{\partial p}}{{\partial {\mathit{\boldsymbol{x}}_i}}} + \frac{1}{\rho }\frac{{\partial {\mathit{\boldsymbol{S}}_{ij}}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} + {\mathit{\boldsymbol{f}}_i} $

(2)

式中：t、u_j、p，r分别表示时间、流速、密度和动水压力，f_i是力源项，S_ij是粘性项，如下:

$ {\mathit{\boldsymbol{S}}_{ij}} = \mu \left( {\frac{{\partial {\mathit{\boldsymbol{u}}_i}}}{{\partial {\mathit{\boldsymbol{x}}_j}}} + \frac{{\partial {\mathit{\boldsymbol{u}}_j}}}{{\partial {\mathit{\boldsymbol{x}}_i}}}} \right) $

(3)

式中μ表示粘性。

柱体结构可简化为一个质量-弹簧-阻尼系统，仅考虑垂直于流向方向的运动，其控制方程为

$ m = \frac{{{\partial ^2}y}}{{\partial {t^2}}} + c\frac{{\partial y}}{{\partial t}} + ky = {F_y} $

(4)

式中：m、c和k分别表示柱体的质量、阻尼参数和刚度，F_y表示垂直于流向方向柱体的合力。

模型在直角坐标系下建立，采用多相流理论将计算域分为固相、液相两个部分。定义体积函数ϕ_m来区分固相和液相(m=1，代表液相；m=2，代表固相)并捕捉固-液界面。ϕ_m须满足：

$ \frac{{\partial {\phi _{12}}}}{{\partial t}} + {\boldsymbol{u}_i}\frac{{\partial {\phi _{12}}}}{{\partial {\boldsymbol{x}_i}}} = 0 $

(5)

式中：ϕ₁₂=ϕ₁+ϕ₂，并在一个网格内满足ϕ₁+ϕ₂=1。网格内的流体特性可用下式来表示：

$ \lambda = \sum\limits_{m = 1}^2 {{\phi _m}{\lambda _m}} $

(6)

式中λ表示密度ρ或粘滞系数μ。

本文所采用的CIP-ZJU模型，核心部分在于对流体的控制方程式(1)和(2)进行离散，其主要过程分为三步：1)采用CIP方法^[20-21]求解对流项；2)采用中心差分法求解除压力项外的非对流项；3)采用SOR方法进行压力速度耦合。固-液边界采用浸入边界方法^[22]处理。

计算区域及网格划分如图 1所示。圆柱直径D=1，两圆柱间距L分别取2D、4D、6D、8D，入口边界和上下边界距离上游圆柱圆心均为10D，出口边界距离上游圆柱圆心35.5D。数值模拟的参数为：雷诺数R_e=150，折合速度U_r=U/f_nD=3~12 (U=1.0为均匀来流流速，f_n为圆柱自振频率)，圆柱的质量比m^*=4m/ρ_fπ D²=2.0，阻尼比ζ^*=0。为了更好的捕捉柱体的运动和其周围的复杂流动现象，对柱体近壁面进行网格加密，最小网格为Δx=Δy=D/40。计算流场的边界条件和初始条件设置如下：初始速度场设定为从左向右并保持均匀来流速度，u=U，v=0 (u、v分别为x、y方向的速度)；初始压力场设定为零；出口边界条件为开边界；柱体表面采用无滑移边界；两侧壁面采用自由滑移边界。

2 柱体绕流和单柱体涡激振动验证

为验证本文中数值模型的正确性，首先对单圆柱绕流和涡激振动的工况分别进行计算。

2.1 网格收敛性验证

为验证网格收敛性，对单圆柱绕流进行计算并比较阻力系数平均值(C_d·mean)和斯特劳哈尔数(St)，为便于与已有文献结果对比，选取雷诺数Re=100。选用三套网格(mesh₁, mesh₂, mesh₃)，对应的最小网格分别为：Δx=Δy=D/50, D/40和D/20。模拟结果如表 1所示，mesh₁和mesh₂与文献[7, 23-24]结果对比吻合良好，综合考虑精度及计算效率，本文采用mesh₂对接下来的工况进行计算。

2.2 单圆柱涡激振动

下面对单圆柱的涡激振动进行模拟。在本节中取质量比m^* =10.0，雷诺数R_e=100，并将振幅比(A^*=A/D，A为柱体振幅)随折合速度变化情况与文献[7, 26]的结果进行比较，如图 2所示。当U_r =5时，振幅比A^*=0.565，达到最大值，此时圆柱已进入频率锁定区域；而随着折合速度增大，振幅比逐渐减小，然而仍保持较大的振幅；当U_r =8.5时，振幅比开始变小接近于0，此时圆柱脱离频率锁定。对比可看出，本文模拟结果与文献[7, 26]结果吻合较好。

3 结果与讨论 3.1 上游圆柱固定

上游圆柱固定情况下，下游圆柱的振幅比如图 3所示。当间距比较小(L/D=2、4)，且折合速度U_r≤4时，下游圆柱振幅非常小，此时圆柱还没有进入共振区域；而L/D较大时(L/D=6、8)，即使U_r≤4，下游圆柱振幅也比小间距比情况下大很多。当U_r>4时，各间距比下振幅变化趋势相同，均急剧增大，并在6≤U_r≤7时达到最大，最大值为0.93，出现在L/D=4，U_r=6，远大于单柱的最大振幅；随着折合速度的增加，下游圆柱的振幅呈减小趋势，但与单柱相比，即使当U_r≥9，下游圆柱依然保持较大的振幅，这是由于受到上游圆柱尾涡影响的结果。

图 4所示为两圆柱阻力系数平均值随折合速度的变化曲线。对于上游圆柱，L/D和U_r均较小时，阻力系数平均值(C_d·mean)也较小，这是由于在这种情况下，下游圆柱振幅较小，剪切层从上游圆柱分离后附着在下游圆柱上，并不产生涡脱。而当L/D和U_r较大时，C_d·mean基本不随折合速度变化，且接近单圆柱绕流的结果(C_d·mean=1.32)。对于下游圆柱，C_d·mean最大值出现在U_r=5~6，基本和振幅最大值出现在同一位置；和单柱涡激振动相比，串列双柱情况下下游圆柱C_d·mean普遍偏小。

图 5给出了个间距比下下游圆柱频率比随折合速度的变化关系，可以看出，在上游圆柱尾流的影响下，下游圆柱几乎没有出现“锁定区间”(频率比并不锁定在1附近)。在小间距比下，当U_r≤6时，两种间距比下圆柱振动频率比变化趋势类似，且L/D=2时频率比稍大；当U_r＞6时，随着折合速度增大，圆柱振动频率比呈直线增加趋势，L/D=4时斜率为0.158，L/D=2则斜率为0.107。大间距比时，两种间距比下频率比变化基本相同，可见当间距比大于6时，间距比变化对频率比影响不大。以U_r=7为界，频率比有两种不同的增长幅度，U_r≤7时曲线斜率较大，约为0.174；U_r＞7时斜率约为0.043。值得一提的是，由于上游尾涡的影响，除L/D=4外，各个间距比下，斯特劳哈尔频率和圆柱固有频率在振动中似乎都不起主导作用。

双柱情况下，流场的涡量分布较单柱涡激振动更加复杂。图 6给出了不同间距比和折合速度下的6个典型的涡量图。图 6(a)所示为L/D=2，U_r =4时的涡量图，可见在间距比较小时下游圆柱和上游圆柱形成一个整体，上游圆柱并不产生涡脱，下游圆柱的位移非常小，其后方出现一条很长的涡带，然后形成涡脱。类似的情况在L/D=4时也会出现，见图 6(c)。间距比小而折合速度很大时，如图 6(b)所示，两柱不再形成一个整体，而是分别出现涡脱，上游圆柱产生的涡脱落后从两柱间隙通过，加剧了下游圆柱的振动，并在下游圆柱后方出现一列稳定的卡门涡街。间距比较大且折合速度很大时，下游圆柱在上游圆柱尾涡的作用下剧烈振动，且后方的涡呈现不规则的形态，如图 6(d)所示。图 6(e)、(f)所示为间距比较大而折合速度较小的情况，此时虽然下游圆柱振幅很大，但其后方形成的涡却比较规则，呈现两条涡街，间距比达到8时，涡的融合现象比较明显。

3.2 上游圆柱自由振动

上游圆柱自由振动情况下，上、下游两圆柱振幅比随折合速度的变化情况如图 7所示。上游圆柱在L/D=2时的振幅要大于单柱的情况，其振幅在U_r=7处达到最大值后，随着折合速度的增大而减小，但相对单柱仍然较大；而随着间距比的增大，上游圆柱的振动响应则非常接近于单柱，这是由于下游圆柱影响减弱造成的。下游圆柱在各个间距比下，U_r＜6时的振幅都较小，而当U_r≥6后振幅变得非常大，最大值为1.082，出现在L/D=2，U_r=7，其共振区域依然没有出现上限，在U_r=12时，振幅仍然很大。在U_r=4, 5时，上游圆柱振幅较大而下游圆柱振幅较小，这是由于此时上游圆柱产生两列涡，下游圆柱的振动被这两列涡抑制所造成的，相应的流场情况见图 10(d)、(e)，这与文献[12]的结果一致。值得注意的是，在L/D=4，U_r=4时后柱的振幅也很大，上游圆柱的尾涡虽然在一定程度上抑制了下游圆柱的振动，但在下游圆柱影响下，上游圆柱振幅减小，其形成的两列尾涡间距减小，从而削弱了对下游圆柱的抑制作用。

图 8给出了两圆柱阻力系数平均值的变化情况。从图 8(a)可以看出，随着折合速度的变化，C_d·mean最大值出现在振幅最大处。当U_r＜6时，各个间距比下C_d·mean均比单柱情况下要小，L/D>2时C_d·mean的变化趋近于单柱，而L/D=2时C_d·mean则比单柱要略大。对于下游圆柱，如图 8(b)所示，其变化趋势和上游圆柱固定时下游圆柱的情况类似。然而在L/D=6、8时，在小折合速度下，C_d·mean也很小，甚至在U_r=4时接近于0，这是因为这种情况下下游圆柱在上游圆柱脱落的两列涡之中运动，流速很小，所以阻力也很小。

图 9所示为两柱的频率比。可以看出，两柱均自由振动时，其频率比变化基本一致，且除L/D=2外，其他间距比下频率比接近于单柱情况。值得注意的是，在L/D=6、8，U_r=6时，频率比有一个较小的突变，这是由于此时下游圆柱的振动模式有微小的变化，同样的现象文献[10, 26]也有记载，文献[26]中称之为“soft lock-in”。

图 10所示为上游圆柱自由振动情况下流场的涡量图。在L/D=2，U_r=3时，如图 10(a)所示，上下游圆柱振幅都很小，上游圆柱不产生涡脱，这与图 6(a)展示的情况类似。当折合速度增大，U_r=9时，如图 10(b)、10(c)所示，上下游的振幅都很大，上游圆柱产生涡脱，涡从上游圆柱脱落后从两柱间隙穿过，形成间隙流，这也与上游圆柱固定的情况类似，但在上游圆柱振动情况下，两柱相对距离增大，涡更容易从间隙通过，使同样的间距比和折合速度下下游圆柱振动更加剧烈。图 10(d)、10(e)所示分别为U_r=4时L/D=6和8的情况，此时上游圆柱振幅较大，在上游圆柱后方形成两列涡并伴随着涡的融合现象，下游圆柱受这两列涡影响，上下两侧压力差较小，其振幅也较小。这种情况在上游圆柱固定的情况下则不会发生。图 10(f)所示为L/D=8，U_r=9的情况，此时间距比和折合速度都较大，上游圆柱受下游圆柱影响小，振动响应接近于单柱，振幅较小，下游圆柱的振动则类似于上游圆柱固定时的情况。

4 结论

1) 双圆柱涡激振动的阻力系数普遍比单圆柱涡激振动时要小，但在L/D=2，U_r＞6时阻力则要大于单柱的情况。圆柱振幅达到最大值所对应的折合速度处，阻力系数平均值也达到最大值。

2) 上游圆柱固定时，下游圆柱的振动频率几乎不由斯特劳哈尔频率控制，即使不在共振区间内。且此时的共振区间很不明显。上游圆柱自由振动时，除L/D=2外，上下游圆柱的振动频率基本相等，且非常趋近于单柱的情况。

3) 上游圆柱自由振动情况下，在4≤U_r≤5时，上游圆柱后方产生两列涡，使下游圆柱的运动范围限制在两列涡之间，对其振动产生了抑制作用，但在某些特定条件下(如L/D=4，U_r=4)这种抑制作用会减弱。