低雷诺数下直圆柱和波型圆柱受迫振动的数值研究

平 焕， 张 凯,2， 周 岱,3,4， 包 艳， 朱宏博， 韩兆龙

(1.上海交通大学 船舶海洋与建筑工程学院，上海 200240； 2.横滨国立大学 都市创新学院，横滨 2408501；3.高新船舶与深海开发装备协同创新中心，上海 200240； 4.上海交通大学 海洋工程国家重点实验室，上海 200240)

自从Bishop等[1]开创性地开展直圆柱的受迫振动实验以来，均匀来流下直圆柱的受迫振动问题成了涡激振动领域一个聚焦的热点[2-4]。相对于自激振动，受迫振动能更多地反映出尾涡与柱体运动之间的相互作用以及能量转移的关系，因此具有重要的研究意义。

在受迫振动的研究中，发现了许多重要的现象，例如：Ongoren等[5]通过实验研究发现，当直圆柱的振动频率(fex)与自然泻涡频率(fs)之比fex/fs从0.85增至1.17时，初始泻涡形成的相位会发生近180°的转变，即在圆柱一侧形成的涡在圆柱运动到了同一侧的极限位置处脱落的现象转变为在圆柱运动到了另一侧极限位置处才脱落。Gu等[6]研究了直圆柱在Re=185和Re=5 000下，振幅为0.2D(D为直圆柱直径)的受迫振动，也发现了类似的现象；受迫振动下，直圆柱尾流区并不是单一的泻涡模式，Williamson等[7]通过大量研究统计，观察到了2S，2P和P+S三种主要的泻涡模式，并在A*-fex/fs平面上进行了划分，同时给出了临界曲线。

另外一个备受关注的现象就是“锁定(Lock-in)”，即当运动圆柱的振动频率接近自然泻涡频率时，圆柱后侧的泻涡频率(fcl)不再由斯特劳哈尔数St确定，而是锁定到振动频率上。但对于泻涡频率锁定到振动频率的具体判定，不同的学者有不同的定义，例如：Karniadakis等[8]研究了Re=100下直圆柱的受迫振动。他们将尾流区各个监测点上的流向速度能量谱的主频等于振动频率的情况视之为锁定；Meneghini等则将升力系数能量谱的主频等于振动频率定义为锁定；Kumar等[9]提出锁定必须同时满足：①升力系数能量谱的主频等于fex；②主频以外的其他频率等于fex的整数倍这两个条件；另外，Koopmann[10]发现振动频率在一定范围时，直圆柱后侧由倾斜泻涡变为平行排列泻涡时，就说明发生了锁定。不同的定义对数据结果的呈现产生了一定的差异，特别是对锁定区间的划分。

近些年研究发现，波型圆柱(直圆柱横截面直径沿轴向按正弦曲线形式变化得到的一种圆柱体)具有较好的减阻和抑制升力的效果，展现出了潜在的工程应用价值，故开始引起学者们的关注。Ahmed等[11]通过实验研究了不同波长波型圆柱的表面压力分布，发现最大截面处的压力系数要大于最小截面处。邹琳等[12]通过数值模拟研究了Re=3 000下不同表面波动幅值的波型圆柱的物理特性，指出较小的表面波幅并不能达到减阻的目的。但随波幅的增加，平均阻力系数开始减小，最大减阻量可达16%。Lam等[13]则系统地研究了Re=100时不同波长和波幅的静止波型圆柱绕流，发现在最优波长和波幅下卡门涡街能够完全得到抑制，同时定义了A，B，C三种泻涡模式。

综合以上国内外研究现状可以看出，目前受迫振动的研究对象都是直圆柱，几乎没有关于波型圆柱受迫振动的研究，同时对两者在受迫振动下物理特性比较的文献也较稀少。因此，本文在低雷诺数下，对均匀来流中作横向正弦受迫振动的直圆柱和波型圆柱进行了数值模拟。以Re=150为例，对比研究直圆柱和波型圆柱在不同振幅和振动频率下的阻、升力和尾涡的变化特征以及锁定区间的确定。

1 数值模型

1.1 几何模型

图1给出了波型圆柱的模型示意图，其直径沿轴向Z的变化由下式给出

Dz=Dm+2acos(2πz/λ)

(1)

Dm=(Dmax+Dmin)/2

(2)

式中：Dz为波型圆柱沿轴向的局部直径，Dm为平均直径，Dmax，Dmin分为沿轴向最大和最小截面所在处的直径，a为波型圆柱表面正弦曲线的幅值，λ为波长。由于a=0时，Dz=Dm=D，波型圆柱即成为直圆柱，故这部分只以波型圆柱为例作为说明。本文采用λ/Dm=2，a/Dm=0.18的波型圆柱研究其在受迫振动下的特性。

图1 波型圆柱的模型示意图

1.2 流动控制方程

三维不可压缩均匀黏性流体运动的基本控制方程为连续性方程和Navier-Stokes方程，在直角坐标系下，分别可以表示为如下的无量纲形式

(3)

(4)

式中：下标i为坐标分量(三维情况下，i=1，2，3为x，y和z三个坐标分量)，ui为i方向的流体运动速度分量，p为压力，Re=U∞Dm/υ为雷诺数，U∞为均匀来流速度，υ为流体的运动学黏性系数。对于圆柱运动引起的网格变形问题则通过求解Laplace方程来解决

Δ·(γ▽dm)=0

(5)

式中:dm为网格位移，γ为到圆柱表面的逆距离控制的扩散系数。本文利用CFD软件OpenFOAM内置求解器pimpleDyMFoam求解上述方程。具体地说，采用有限体积法求解N-S方程：对时间项采用隐式积分方法；对流项采用二阶迎风离散格式；对控制方程中的速度和压力耦合则是采用PIMPLE算法来处理。PIMPLE算法是SIMPLE和PISO算法的结合体，有关该算法的具体介绍，请参阅相关文献[14-15]。

1.3 计算区域及边界条件

1.3.1 计算区域

采用圆柱体型计算区域，坐标原点位于圆柱底面的中心，计算区域半径为30Dm，高度为2Dm。Lam等[16]的研究表明，无论流动是处于层流还是湍流，选取一个波长的计算区域高度可以得到准确的模拟结果。另外，本文对最大振动频率下的工况(A*,f*)=(0.25，0.3)进行了验证，发现采用两个波长的计算区域高度得到的结果与采用一个波长时一致。基于上述理由，本文选取2Dm(一个波长)的计算区域高度能够保证获得准确的计算结果。沿x-y平面，网格呈非均匀分布，靠近圆柱表面网格密集，远离圆柱，网格逐渐变稀疏。而沿圆柱的轴向Z，网格均匀分布。波型圆柱周围的网格分布，如图2所示。

图2 波型圆柱周围的网格分布

1.3.2 边界条件

1.4 无量纲物理量的定义

本文中用到的一些重要的无量纲物理量定义，如表1所示。

2 网格及算法可靠性验证

对直圆柱分别采用较为稀疏的网格模型mesh I、中等疏密的网格模型mesh II和较为稠密的网格模型mesh III进行网格可靠性验证，相应的网格参数，见表2。

表3给出了在Re=100下静止直圆柱绕流计算得到的升力系数均方根值CL(rms)、平均阻力系数CD(mean)、斯特劳哈尔数St以及与其他已有参考文献结果的比较。可以发现，使用网格模型mesh II和III得到的计算结果与文献结果较为接近，说明这两种网格模型均满足计算精度的要求。考虑到计算效率，对于直圆柱，选用网格模型mesh II作为后续计算的网格模型。

表1 无量纲物理量的定义

表2 直圆柱网格模型参数

表3 直圆柱网格模型计算得到的结果和已有文献结果的比较

为进一步验证算法的正确性，本文对直圆柱作横向受迫振动的情况进行了验证。研究了Re=185，A*=0.2下6种不同振动频率的工况。图3为本文数值结果与Guilmineau等[20]结果的对比，其中CD(rms)为阻力系数均方根值。可以看出，两者结果相吻合且具有一致的变化趋势，说明本文采用的算法适用于低雷诺数下的受迫振动。

另外，我们对波型圆柱计算网格模型的可靠性进行了验证。同样采用三种不同疏密程度的网格模型mesh IV、mesh V、mesh VI，相应的网格参数见表4。

图4则给出了波型圆柱在(A*,f*)=(0.25, 0.15)时，通过三种网格计算得到的阻、升力系数的时历曲线。可以看到，三条时历曲线基本一致。计算表明，由mesh IV获得的CL(rms)、CD(rms)和CD(mean)相对于mesh VI，其误差分别为4.9%、1.9%和0.6%；而由mesh V获得的结果相对于mesh VI，误差则分别为1.9%、1.8%和0.2%，说明计算结果在三种网格上趋于网格收敛。考虑到计算效率，对于波型圆柱，选用网格模型mesh V作为后续计算的网格模型。

图3 受迫振动下直圆柱计算得到的结果和已有文献结果的比较

Fig.3 Comparison with published results for uniform circular cylinder undergoing forced oscillation

表4 波型圆柱网格模型参数

图4 不同波型圆柱网格模型计算得到的阻、升力时历曲线

Fig.4 Time histories of drag and lift coefficients with three different grid resolutions for wavy cylinder

3 计算结果及分析

3.1 静止波型圆柱绕流的数值模拟

图5为Re=150时，模拟静止波型圆柱绕流得到的阻、升力系数时历曲线。从图5可知，阻力系数无波动，升力系数为零，说明涡脱落现象被完全抑制，卡门涡街型泻涡频率消失。图6则为计算稳定后某一时刻的三维涡结构图。从图6可知，两侧对称分布的自由剪切层稳定地朝下游延伸发展，没有卷起形成涡对。故本文采用的波型圆柱在静止状态下能完全抑制泻涡的发生。

图5 静止波型圆柱绕流的阻、升力时历曲线

图6 静止波型圆柱绕流的三维涡结构图(ωz)

Fig.6 Instantaneous vortical structure for flow past a stationary wavy cylinder

3.2 受迫振动下直圆柱和波型圆柱绕流的数值模拟

3.2.1 阻、升力系数

图7给出了不同振幅下，直圆柱和波型圆柱的CL(rms)，CD(rms)以及CD(mean)随振动频率的变化情况。对于直圆柱，研究了区间为 [0.05,0.5]的10种不同振幅；对于波型圆柱，研究了A*=0.05, 0.25, 0.5三种有代表性的振幅，分别对应着小振幅、中振幅和大振幅的情况。

从图7(a)可知，在f*<0.18的低频范围内，随着振动频率的增加，直圆柱的升力系数均方根值曲线先下降后上升，但不同振幅的曲线开始下降的点所对应的振动频率不同，其随振幅的增加而减小。通过对比图8可知，该振动频率即为该振幅下系统开始进入锁定状态的频率，即升力系数均方根值的下降是由系统开始进入锁定状态引起的。在f*>0.18的高频范围内，升力系数均方根值对振幅的依赖性较大，振幅越大引起的升力系数幅值也越大。另外，值得注意的是，在振动频率区间[0.18,0.2]内，各个振幅的曲线有一个突变，通过对比图8可知，这是由于系统脱离锁定状态引起的。对于波型圆柱，从图7(d)可知，A*=0.5与A*=0.25时的曲线走势相似。初始阶段由于振动频率小，振动十分缓慢，接近于静止状态，而静止状态下波型圆柱的泻涡是被完全抑制的，故升力系数均方根值不变且较小，升力是由圆柱的受迫振动引起。接着，由于泻涡的发生，升力系数均方根值出现了小的陡增。随着振动频率的进一步增加，系统开始进入锁定状态(对比图9)，导致曲线有一个明显的下降。最后，伴随着系统脱离锁定状态，升力系数均方根值出现突增。A*=0.05时的情况也类似，只是曲线对应阶段的变化幅度比前两者小。但有所不同的是，在f*=0.25, 0.26和0.3处波型圆柱后侧的泻涡被完全抑制，导致其升力值明显小于相邻值。对比CL(rms)纵坐标值可以发现，波型圆柱较直圆柱在低频段体现出了抑制升力的效果。

振幅和振动频率对阻力系数均方根值的影响示于图7(b)和图7(e)。对于直圆柱，总体上，振幅越大，平均阻力系数也越大。每条曲线经历了两次下降过程。第一次下降代表系统进入锁定状态；第二次下降代表系统脱离锁定状态；从图7(e)可知，波型圆柱的情况与直圆柱类似。A*=0.5与A*=0.25时，曲线的两次下降过程同样分别对应着系统进入和脱离锁定状态。A*=0.05时的曲线在高频段的表现则略有差异：随着振动频率的增大，曲线并没有呈现上升趋势，反而由于在f*=0.25，0.26和0.3处由于泻涡被抑制，导致出现下凹。对比CD(rms)纵坐标值可以发现，波型圆柱在高频段体现出了减阻的效果。

图7(c)和图7(f)则展示了平均阻力系数与振幅和振动频率之间的关系。图7(c)曲线的总体走势和图7(b)大致相似，但初始段不同振幅对平均阻力系数几乎没有影响；对于波型圆柱，通过对比图7(f)和图7(e)可知：平均阻力系数曲线少了第一次下降过程。

(a) 直圆柱CL(rms)

(b) 直圆柱CD(rms)

(c) 直圆柱CD(mean)

(d) 波型圆柱CL(rms)

(e) 波型圆柱CD(rms)

(f) 波型圆柱CD(mean)

3.2.2 锁定现象

采用Kumar等对锁定的定义，即：①升力系数能量谱的主频等于fex；②主频以外的其他频率等于fex的整数倍。不同时满足上述两个条件的情况则视为非锁定。

图8给出了直圆柱的锁定区间和A*-f*平面上4个特定工况的升力系数时历曲线及对应的能量谱。分界线下方区域代表非锁定，上方区域代表锁定。对于一个固定振动频率，当振幅超过某个临界值时，系统发生锁定，临界振幅随振动频率远离直圆柱的自然泻涡频率而增加，这与Koopmann、Anagnostopoulos[21]的研究结果相似。

在f*<0.18的区间，分界线变化较平缓。图8(a)给出了(A*,f*)=(0.05, 0.14)的升力系数时历曲线及其对应的能量谱。从图8(a)可知，升力并不是按规律的正弦形式波动，而是存在一个“差拍”现象。对应的能量谱显示，主频并不等于外加振动频率，而是出现在自然泻涡频率附近，故处于非锁定状态。当振幅增加到0.25时，此时主频已经固定到了外加振动频率处，且其他频率等于外加振动频率的整数倍，这从图8(b)可知，系统由之前的非锁定状态转变到了锁定状态，这时升力系数时历曲线也呈现出了规律的正弦波动形式。

在f*>0.18的区间，分界线较陡。图8(c)和图8(d)对应着(A*,f*)=(0.05, 0.19)和(A*,f*)=(0.25, 0.19)的工况，分别处于非锁定和锁定状态。值得注意的是，由于(A*,f*)=(0.05, 0.19)的工况靠近过渡带，故能量谱出现了两个幅值相当接近的频率(见图8(c)局部放大图)，但由于主频与外加振动频率并不重合，故仍处于非锁定状态。

(a) (A*, f*)=(0.05，0.14)

(b) (A*, f*)=(0.25，0.14)

(c) (A*, f*)=(0.05，0.19)

(d) (A*, f*)=(0.25，0.19)

图9给出了波型圆柱在振幅为0.05、0.25和0.5时的锁定区间以及A*-f*平面上4个特定工况的升力系数时历曲线和对应的能量谱，同时对比了相同条件下直圆柱的结果。从图9可知，在图示频率区间内, 大部分波型圆柱的锁定区间被包络在直圆柱的锁定区间内。波型圆柱左侧分界线的走势和直圆柱相似，然而右侧却不同，在A*=0.25处有一个明显的转折。

(a) (A*, f*)=(0.05，0.15)

(b) (A*, f*)=(0.25，0.15)

(c) (A*, f*)=(0.05，0.18)

(d) (A*, f*)=(0.25，0.18)

图9(a)、图9(b)、图9(c)、图9(d)分别对应(A*,f*)=(0.05,0.15)，(0.25,0.15)，(0.05,0.18)和(0.25,0.18)四个特定工况。升力系数能量谱显示，图9(b)、图9(d)处于锁定状态，图9(a)、图9(c)处于非锁定状态。与直圆柱的情况相同，在锁定状态下，波型圆柱升力系数时历曲线为规律的正弦波动形式，而非锁定状态下则捕捉到了明显的“差拍”现象。除了(A*,f*)=(0.25,0.15)的情况，波型圆柱较直圆柱起到了良好的抑制升力的效果，在(A*,f*)=(0.05,0.15)，(0.05,0.18)和(0.25,0.18)时升力系数均方根值分别降低了43.75%、49.82%和32.20%。

图10、图11分别为A*=0.25时，直圆柱和波型圆柱在各种振动频率下的升力谱。其中纵坐标mag代表对升力系数的时程信号做快速傅里叶变换(FFT)后，获得的以某种频率成分振动的幅度特性。

图10 A*=0.25时直圆柱在各种振动频率下的升力谱

图11 A*=0.25时波型圆柱在各种振动频率下的升力谱

对于直圆柱，大部分范围内，升力谱由外加振动频率和卡门涡街型泻涡频率共同控制。但在振动频率区间[0.14,0.19]内，代表卡门涡街型泻涡频率的竖线中断消失，由振动频率唯一控制，即发生锁定；对于波型圆柱，同样可以清楚地发现卡门涡街型泻涡频率的存在，这说明在静止波型圆柱绕流中消失的卡门涡街型泻涡频率在受迫振动中又得以显现。但不同于直圆柱，其发生了两处中断，分别在[0.14,0.18]和[0.05,0.08]区间上。前者代表着锁定，而后者则是由于不泻涡导致的。

3.2.3 尾涡特性

重点观察了锁定区间内圆柱尾流区的泻涡模式及尾涡形态。通过统计发现，直圆柱后侧的泻涡模式由振动频率控制，振幅则几乎没有影响。振动频率较小时(低频段)，泻涡呈现2S模式；振动频率较大时(高频段)，则呈现C(2S)模式，这种泻涡模式在Singh等[22]的研究中也有发现。以振幅A*=0.25为例，图12给出了直圆柱处于平衡位置并朝y轴正方向运动时刻尾流区的泻涡形式。

f*=0.14下，直圆柱尾流区呈现出典型的2S模式，即直圆柱在运动到最高点和最低点时各泻放一个旋转方向相反的涡，如图12(a)所示。

当振动频率增加到f*=0.19时,如图12(b)所示，在近尾流区，同f*=0.14一样，其呈现出典型的2S模式；但在离开直圆柱一定距离后，旋转方向相同的涡发生合并，呈带状向下流延伸，即C(2S)模式。同时尾迹区的侧向宽度随离直圆柱距离的增加而增加。

(a) 2S泻涡模式

(b) C(2S)泻涡模式

Fig.12 Instantaneous vorticity contours showing two kinds of modes for uniform circular cylinder in lock-in

对于所有处于锁定状态的波型圆柱，其后侧的尾涡呈现出了大致相同的特征：总体看来，和直圆柱泻涡形式类似，但在轴向出现了轻微的扭曲，开始显现出三维效应。最大和最小截面的泻涡形态并没有明显的不同，均呈现出2S模式。以(A*,f*)=(0.25,0.18)为例，图13分别给出了三维尾涡结构图以及最大和最小截面处的泻涡模式。这种尾涡呈现出来的特征和Lam等[13]提出的I类波径比下的泻涡模式A十分相似。

4 结 论

本文利用CFD软件OpenFOAM对比研究了Re=150下直圆柱和波型圆柱的横向受迫振动，特别是对波型圆柱受迫振动的研究具有创新意义。通过改变振幅和振动频率，观察阻、升力的变化特征，确定锁定区间的范围以及锁定状态下的尾涡特性，得到的主要结论如下：

(1) 总体上，不同振幅下的波型圆柱和直圆柱的升力系数均方根值、阻力系数均方根值以及平均阻力系数随振动频率变化的规律大致相同。波型圆柱较直圆柱体现出了一定的减阻和抑制升力的效果。减阻主要体现在高频段，抑制升力则主要体现在低频段。

(2)大部分波型圆柱的锁定区间被包络在直圆柱的锁定区间内。在锁定状态下，两者的升力系数时历曲线均呈现出规律的正弦波动形式，而非锁定状态下，均捕捉到了明显的“差拍”现象。

(a) 三维涡结构图(ωz)

(b) 最大截面处的泻涡模式

(c) 最小截面处的泻涡模式

(3)在静止波型圆柱绕流中消失的卡门涡街型泻涡频率又显现在受迫振动中。

(4)锁定区间内，在直圆柱后侧观察到了两种泻涡模式。泻涡模式由振动频率控制，与振幅几乎无关。振动频率较小时，呈现2S模式；振动频率较大时，呈现C(2S)模式。而对于波型圆柱，在锁定区间内只观察到一种泻涡模式，该模式和Lam等提出的I类波径比下的泻涡模式A十分相似。

(5)本文的工作为今后进一步开展均匀来流中波型圆柱涡激振动的研究作了铺垫与参考。