平动与转动受迫谐振圆柱的水动力特性分析

刘名名， 唐国强， 吕 林， 滕 斌

(大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024)

平动与转动受迫谐振圆柱的水动力特性分析

刘名名， 唐国强， 吕 林， 滕 斌

(大连理工大学 海岸和近海工程国家重点实验室，辽宁 大连 116024)

通过求解二维不可压缩黏性流体的Navier-Stokes方程，对低雷诺数下圆柱同时具有周期性旋转和横流向振动的流固耦合问题开展了数值分析研究。重点考察了不同强迫振动频率(f*([0.5, 2.0])以及转动与平动之间的相位差(φ=±180°，±120°，±60°及0°)下，圆柱结构的受力特性、锁定模式及尾迹涡流场演化特征。数值分析结果表明，在周期性转动和平动的联合激励下，流动过程除存在常规的基本锁定(Primary Lock-in)现象外，还存在非线性准周期锁定(Quasi-periodic Lock-in)现象；旋转运动与横流振动之间的相位差φ对锁定形式及锁定的频率区间有重要影响，对称的正负相位差可导致高度一致或截然不同的锁定形式；两自由度谐振圆柱的水动力系数在锁定条件下会发生明显跳跃，总体上随强迫振动频率的增大而增大；对于特殊的准周期锁定情况，尾迹区内存在涡旋结构的多级次非线性分叉行为，尾涡演化过程遵循准周期特性。

圆柱；旋转振动；横流向振动；准周期锁定；尾涡模式

当流体经过非流线型结构时，通常会在其尾流区内产生交替的涡旋脱落，导致作用在结构上的流体作用力表现出周期性特征，这是诱发结构发生涡激振动以及驰振等非线性动力响应(甚至是疲劳破坏)的主要原因之一。圆柱结构在水利工程、海洋工程、化学工业以及航空航天等领域都有十分广泛的应用。理解和认识圆柱结构与流体的相互作用，特别是多自由度振动圆柱结构的受力以及相应的流动特性，无论是对解决工程实际问题(包括流动控制)，还是认识自然规律，都具有重要的科学研究价值。

对结构施加一定运动规律的强迫振动是研究流固耦合问题的重要方法，这样允许针对特定的运动形式，揭示结构的受力、流动特性以及流固耦合作用的机理。在以往的研究中，根据圆柱结构受迫振动方式的不同，可划分为平动(包括横流向及顺流向振动)和转动(包括单向和双向)两种基本形式。对于只具有横流向单自由度受迫简谐振动的问题，Williamson等[1]通过物理模型试验，总结了不同的尾涡脱落模式(如2S、2P及P+S等)及其与圆柱振动幅值和振动频率之间的依赖关系。在后来的研究工作中，很多学者关注了锁定现象的发生，以及锁定区间与圆柱受力及尾涡脱落模式之间的关系[2-6]。横流向受迫振动的研究成果目前已经成为细长柔性结构涡激振动经验预报模型的基础[7-8]。对于圆柱发生顺流向强迫振动的问题，人们同样也针对不同振幅和频率下的结构受力、锁定特性以及尾迹模式等开展了大量研究，并取得了丰硕的研究成果[9-12]，这对于认识水流和波浪联合作用下的结构受力和稳定性等都具有直接的借鉴意义。相比较而言，针对圆柱的受迫旋转运动问题，相关的研究工作要少一些。通过强迫圆柱的旋转运动，可改变边界层内的流动分离特性，从而达到流动控制的目的。Du等[13]发现通过旋转的方法可以抑制圆柱涡激振动。Bai等[14-15]对旋转圆柱的流固耦合问题开展了数值模拟研究，揭示了圆柱旋转速率对圆柱受力及流场特性的影响作用。Lu等[16]则提出了通过圆柱的主动性双向旋转来减小圆柱受力的流动控制方法。此外，学者们针对均匀流中圆柱发生受迫双向旋转的问题，也给出了圆柱受力和锁定区间等对旋转幅值和受迫振动频率的依赖关系[17-18]。

上述的这些研究工作主要是针对具有单自由度的圆柱受迫振动问题开展的，而对于圆柱在水流中同时具有多个自由度强迫运动的研究工作则相对较少。Blackburn等[19]在横流向受迫振动圆柱的基础上引入圆柱旋转运动，开展了流动控制研究。Nazarinia等通过模型试验，研究了单向均匀流中，圆柱同时具有横流向振动与转动时的受力特性，并观察到了非线性准周期锁定的现象。但是，由于受到试验条件和观测方法的限制，Nazarinia等[20]没有详细给出发生准周期锁定时的流场演化特征，同时也没有进一步研究圆柱的振动频率对其受力和锁定特性的影响作用。因此，针对同时具有平动和转动的多自由度圆柱受迫振动问题的研究还有待进一步加强，这一物理过程的背后具有丰富的非线性流固耦合现象，理解和认识其中的物理机理，也是今后开展更为深入的流动控制研究的基础。本文将通过数值模拟，对同时具有横流向振动和旋转振动圆柱的水动力特性开展分析，重点考察圆柱受迫振动频率以及横流振动与旋转振动之间的相位差对圆柱结构受力、锁定特性以及涡流场演化特征的影响作用。

1 数值模型

1.1 流体运动控制方程

描述均匀不可压缩黏性牛顿流体运动规律的基本控制方程为连续性方程和Navier-Stokes方程。在任意拉格朗日-欧拉(ALE)观点下，可写成如下的无量纲形式

(1)

(2)

式中：xi表示笛卡尔坐标(对于二维问题，规定x1=x，x2=y)；ui表示流体速度在xi方向的分量；p为压力；t表示时间；Re=UD/υ为雷诺数(υ为流体运动学黏性系数，U为均匀来流速度，D为圆柱直径)；cj表示ALE观点下的j方向的网格运动速度分量。

1.2 圆柱强迫运动方程及网格更新方法

设圆柱同时发生周期性的横流向振动及绕轴心的旋转运动，其无量纲的平动位移y(t)及角位移θ(t)分别为

(3)

(4)

式中：φ2和φ1分别为旋转运动与横流向振动的初始相位，φ=φ2-φ1表示旋转运动与横流向振动之间的相位差；f*=fe/fs为无量纲的强迫运动频率，其中fe表示圆柱的振动频率，fs为相同雷诺数下固定圆柱的涡脱落频率；A*=A0/D(A0为圆柱平动位移幅值)和θ*=θ0U/2πfe(θ0为圆柱角位移幅值)分别表示无量纲的平动幅值与转动幅值。

流体与结构之间的相互作用通过动网格方法进行模拟，相应的网格坐标更新，通过求解如下的Laplace方程来获得

(5)

式中：Δ=∂2/∂x2+∂2/∂y2为Laplace算子；q1和q2分别表示网格发生运动后对应的x和y坐标。边界条件为已知的圆柱表面计算节点坐标。

1.3 数值求解方法

采用三步有限元方法对流动动量方程进行数值离散，具体形式如下[21]

(6)

(7)

(8)

由于网格更新是在一个完整时步完成后进行的，因此，式(6)～(8)中的网格运动速度cj均采用了n时刻的运动速度。数值计算结果表明，上述的线性化处理方法能够在保证数值精度的前提下，获得较高的计算效率。此外，由式(8)可知，要获得n+1时刻的速度场uin+1，首先需要求得n+1时刻的压力pn+1。通过对式(8)两端同时取散度，并应用n+1时刻的连续性方程∂uin+1/∂xi= 0，可得到如下的压力泊松方程：

(9)

本文所采用的三步有限元方法是一种高阶迎风Taylor-Galerkin格式，并且流速和压力变量通过投影方法进行解耦，因此在利用标准有限元方法进行空间离散时，速度和压力可采用同阶插值函数。同时，三步有限元方法满足统一的CFL(Courant-Friedrichs-Lewy)稳定性条件，本文采用如下的动态时间步长：

(10)

式中：Se为网格面积；ue为网格中心点速度；Min(*)表示对所有计算网格取最小值；β=0.2为安全系数。对圆柱受到的压力和黏性剪切应力沿着圆柱表面进行积分，将合力分别投影到x,y方向上并通过除以0.5ρDU2对其进行无因次化处理，可以得到单位长度的圆柱受到的升力系数和拖曳力系数分别为[22-23]

(11)

(12)

式中：p∞为无穷远处的参考压力。

2 计算模型及验证

2.1 计算域及定解条件

图1给出了圆柱发生转动和横流向两自由度受迫振动时的示意图及边界条件。坐标原点与圆柱圆心重合，整个流体计算域的范围为90D×145D，圆柱距离来流端60D，距离出口端85D，距离上下边界分别为45D。

数值模型的边界条件如下：在入口处指定无量纲速度u=1，垂向速度v=0；侧壁采用对称边界条件∂u/∂y=0，v=0，∂p/∂y=0；出口处指定自由出流边界条件∂u/∂x=0，∂v/∂x=0及相对压力p=0。圆柱表面施加不可滑移边界条件：u=-rωsinθ，v=dy/dt+rωcosθ及∂p/∂n=0，其中，ω=dθ/dt为旋转角速度。在初始时刻，流场中的速度及压力分布均设为零。

图1 计算域及边界条件

2.2 网格收敛性验证

表1 网格收敛性验证结果

2.3 数值模型验证

本文首先针对圆柱仅发生横流向受迫振动的问题开展了数值模拟，并与文献[6]的数值结果进行对比验证。计算中，Re=185，A*=0.2，f*=0.8、0.9、1.0、1.1、1.12和1.2。图2给出了本文计算所得流体力系数与文献[6]结果的对比情况。从图2可以看出，本文的数值结果与已发表论文结果吻合良好，说明本文所建立的数值分析模型能够正确模拟圆柱横向受迫振动流固耦合问题。

为进一步说明本文数值模型的可靠性，对Re=100，f*=2.423，θ*=1.114条件下，受迫旋转圆柱在均匀流场中的绕流问题开展了数值模拟。图3给出了本文计算得到的升力系数时变曲线与Choi等[18]数值结果的对比。从图3可以看出，二者具有良好的一致性，进一步说明了本文模型的正确性。

图2 流体力系数随圆柱横流向受迫振动频率的变化

图3 受迫旋转圆柱的升力系数时间历程对比

3 数值计算结果与分析

在以下的计算中，将圆柱的旋转振动幅值和横流向平移振动幅值分别取为θ*=1.0和A*=0.5(此时圆柱表面的最大切向速度与圆柱的最大平动速度具有相同的数值)，重点考察强迫运动频率f*(fe/fs)以及旋转振动与横流向振动之间的相位差φ的影响作用。计算中，雷诺数Re=100保持不变，相应雷诺数下固定圆柱涡脱落频率为fs=0.165 1。

3.1 圆柱受力特性及锁定分析

图4以相位差φ=-60°，强迫振动频率f*=0.55、0.80、1.20及2.00为例，给出了升力系数和圆柱横流向振动位移的时间历程线，对应的升力系数频谱分析结果，如图5所示。对比图4结果可以看出，随着圆柱振动频率f*的增大，升力系数的幅值显著增加，说明圆柱振动频率是影响圆柱受力的一个重要因素。从图4(a)中可以看出，当圆柱的振动频率f*=0.55时，升力系数的时间历程线不规则。相应的频谱分析结果图5(a)表明，此时的升力频率成份同时包含了圆柱的振动频率fe及固定圆柱的涡脱落频率fs以及由二者之间相互作用所产生的新频率成份(fs-fe)/2及(3.0fe+fs)/2等。当f*=0.80时，从图4(b)中可见，此时圆柱所受升力变化基本满足简谐振动的形式。图5(b)的升力系数频谱分析表明，相应的涡脱落频率(也即升力主控频率)与圆柱的振动频率相等，这意味着圆柱后方的涡旋脱落主要由圆柱的振动所控制，此时尾流场发生锁定现象。此外，虽然图5(b)的频谱中包含了3.0fe频率成份，但其对升力贡献很小，因此，此时的升力振荡主要表现为周期性的简谐振动的形式。对于f*=1.20的情况，图4(c)的结果表明，此时的升力系数具有明显的多频振动特性，并且每5个圆柱横流向振动周期对应1个完整的升力系数振荡周期，即圆柱后方发生一次完整的尾涡脱落周期为圆柱振动周期的5倍。在本文中，将尾涡脱落周期与圆柱振动周期相等的工况定义为基本锁定(Primary Lock-in)，这与传统的圆柱涡激振动锁定的定义是一致的[24]；而将完整的尾涡脱落周期为N(N>1)倍圆柱振动周期的情况定义为准周期锁定(Quasi-periodic Lock-in)。对应于前述的图4(c)的情况，准周期锁定发生在5T的周期上。与图4(c)的情况相类似，当f*=2.0时，尾流场也存在准周期锁定现象。如图4(d)所示，此时每3个横向(旋转)振动周期，升力系数出现周期性重复特征，对应的准周期锁定周期为3T。图5(c)和(d)分别给出了与图4(c)和(d)相对应的升力系数的频谱分析结果。从图4中可以看出，对于这两个发生准周期锁定的情况，相应的涡脱落频率分别为fe/5和fe/3，即fe/N(N=3, 5，…)，而不是圆柱的振动频率fe，但升力振荡的主控频率均对应圆柱的振荡频率。同时，在升力变化以及涡旋脱落的过程中，对应于N倍准周期锁定，还存在(2M-1)fe/N的次谐波频率成份(M=1, 2, …)。

根据非线性振动力学的基本理论，仅依靠频率关系来判断是否发生锁定是不全面的。为进一步证实上述的基本锁定和准周期锁定现象，图6给出了与图4相对应的升力系数与圆柱横流向振动位移的相图。从图6(a)可以看出，由于圆柱升力系数振动的不规则性(未锁定，参见图4(a))，升力系数和圆柱振动位移的相图呈现出杂乱无章的特征，构成典型的非周期性混沌运动状态。图6(b)给出了发生基本锁定时的相图，从图中可以看出，由于升力系数的振动周期与圆柱振动周期一致，相图为严格封闭的单环结构，这充分证实了基本锁定现象的发生。图5(c)和(d)分别给出了在fe/5和fe/3频率下，发生准周期锁定时的相图。从图6中可以看出，此时的相图也表现为封闭的圆环。但是，由于升力振荡周期分别为圆柱振动周期的5倍和3倍，因此对应的相图中也同时包含了5个和3个封闭环。

(a) φ=-60°, f*=0.55

(b) φ=-60°, f*=0.80

(c) φ=-60°, f*=1.20

(d) φ=-60°, f*=2.00

(a)

(b)

(c)

(d)

(a) φ=-60°, f*=0.55

(b) φ=-60°, f*=0.80

(c) φ=-60°, f*=1.20

(d) φ=-60°, f*=2.00

Fig.6 Phase diagram of lift coefficient and transverse displacement

另外，对于图6(b)所示的在f*=0.80情况下的基本锁定情况，由于相图的主轴左向倾斜，圆柱横流向振动位移和升力系数之间的相位差Ψ介于π/2和π之间，这表明流体力起到阻尼的作用。对于一般的涡激振动问题而言，当Ψ∈(π/2, π)时，流体作用力将起到抑制圆柱涡激振动响应的作用。但是，随着圆柱振动频率的进一步增大，如图6(c)所对应的f*=1.20及图6(d)所对应的f*=2.00的情况，圆柱横流向振动与升力系数之间的相位差Ψ介于0到π/2之间(相图主轴右向倾斜)，此时流体向圆柱输入能量，流体力起到激振力的作用[3]。应该注意的是，对于上述的f*=1.20和2.00的情况，二者均对应准周期锁定。

本文通过综合考虑升力系数的频谱特性以及升力振荡与圆柱振动之间的相图来判断基本锁定和准周期锁定的发生与否。表2全面给出了在本文所研究的31个受迫振动频率(f*=0.5～2.0，频率间隔0.05)以及7个圆柱转动与横向振动之间相位差(φ=±180°, ±120°, ±60°和0°)下，共计217组数值计算所得到的锁定判别结果。此外，为了便于对比分析，表中也一并列出了圆柱在相同雷诺数(Re=100)、相同振动幅值(A*=0.5)以及相同振动频率下，仅发生横流向受迫振动的锁定区间判别结果。应该说明的是，对于圆柱仅发生横流向振动的情况下，只存在基本锁定的现象，相应的锁定频率区间范围为[0.7, 1.2]，这与以往的研究成果是一致的[8]。

表2 不同相位差φ及振动频率f*下的锁定区间

Tab.2 Lock-in regime at various phase differencesφand forced oscillatory frequencyf*

f∗θ∗=0φ(θ∗=1)180°120°60°0°-60°-120°-180°0．50×××××13T7T×0．55××××4T×××0．60×××■××■×0．65×■×■■■■■0．70■■■■■■■■0．75■■■■■■■■0．80■■■■■■■■0．85■■■■■×■■0．90■■■■■×■■0．95■■■■■×■■1．00■■■■■××■1．05■■■■■×11T■1．10■×■■■×××1．15■6T■■■××6T1．20■×■■■5T11T×1．25××■■■5T2T×1．30××■■■×2T×1．35××■■■×2T×1．40×5T■■■2T×5T1．45××■■■2T11T×1．50××■■■19T××1．55××■■■×7T×1．60××■■■×××1．65××■■■×××1．70×2T■■■7T3T2T1．75××■■■×3T×1．80×17T■■■×3T17T1．85××■■■×3T×1．90××■■■×3T×1．95××■■■3T3T×2．00××■■■3T××注：×，非锁定；■，基本锁定(尾涡脱落周期与圆柱的振动周期T一致)；NT，准周期锁定(完整的尾涡脱落周期为圆柱振动周期T的N倍)

从表2的结果中可以看出，当引入圆柱的旋转运动后，不但会诱发准周期锁定现象，同时，圆柱转动与横流向振动之间的相位差φ对锁定形式以及锁定的频率区间都有十分重要的影响作用。表2的分析结果表明：① 当0°≤φ<180°时，锁定现象主要以基本锁定形式为主；当-180°<φ<0°时，发生准周期锁定的概率显著增加；② 当φ=±180°时，无论是基本锁定，还是准周期锁定，正负两个相位差下所得到的数值计算结果保持严格一致；③ 当φ=±60°和±120°时，正负两种相位差下的锁定形式和锁定区间存在较大的差别。

对于φ=60°和120°的情况(正的相位差)，基本锁定区间从低频到高频连续分布，未出现准周期锁定现象，其中φ=60°的情况，取得最大的基本锁定区间f*∈[0.6, 2.0]；当φ=-60°和-120°时(负的相位差)，基本锁定区间主要集中在低频区，而准周期锁定主要发生在高频区域。特别的，对于φ=-60°和-120°这两种情况，即便圆柱的振动频率与相同雷诺数下固定圆柱的涡脱落频率相等，即在f*=fe/fs=1.0的条件下，也未发生锁定。当φ=-60°时，基本锁定区间具有最小的频率范围，f*∈[0.65, 0.8]。对于同时具有旋转与平移耦合振动的圆柱结构而言，其锁定的发生条件以及锁定形式都与这两种运动之间的相位差密切相关，特别是在对应的正负相位差下所表现出的显著性差异，可能是非线性迟滞(hysteresis)效应的体现，其中的物理机理还有待于进一步的深入研究。

(a) 升力系数均方根

(b) 平均拖曳力系数

对于相位差为φ=-120°和-60°的特殊情况，由于它们在低频区内都存在锁定区间的明显变化，因此在相应的频率范围内，也都出现了升力幅值的跳跃现象。可见，锁定形式的转换是导致升力发生变化的原因之一。在各相位差下，平均拖曳力系数随振动频率的变化情况相对比较复杂。从图7(b)中可以看出，当φ=120°时，平均拖曳力系数在各频率下几乎都是最小的，且在整个频率范围内变化不大。而与之对应的φ=-120°，其平均拖曳力则一直维持在较高的水平。从图7(b)中可见，相对于仅有横流向振动的情形，当φ=±180°，-120°，-60°及0°时，无论是锁定区间还是非锁定区间，引入圆柱旋转运动后，平均拖曳力系数均会有所增大。参照表2的结果可以发现，图7(b)中有关平均拖曳力系数发生跳跃变化的频率几乎都对应着锁定模式的改变。

3.2 基本锁定与准周期锁定下的涡脱落模式

图8以φ=-60°，f*=0.8的工况为例(对应图4(b))，给出了在基本锁定情况下(即涡脱落频率与圆柱的振动频率相等)，一个完整涡脱落周期内涡量场的变化情况。图中，白色表示正涡量，黑色表示负涡量；N0时刻对应圆柱的横流向振动位移正处于由负转正的平衡位置，各图之间的时间间隔为1/4横向振动周期。

从图8中可见，在N0时刻，圆柱背流侧的下方和上方分别形成正涡A和负涡B，其中正涡A的尺度和强度明显大于负涡B。由于圆柱旋转自由度的存在，在圆柱的右下表面可见明显的剪切层。随着时间的推移，正涡A在N1时刻发生脱落，并在N2时刻形成新的贴体正涡C。而负涡B在圆柱运动与来流的共同作用下，发生显著扭曲，同时在圆柱右上表面形成明显的负涡量剪切层。N2时刻的涡量场与N0时刻构成反对称结构。在N3时刻，负涡B发生脱落，并形成了新的贴体负涡D。同样，N3时刻的涡量场与N1时刻也呈反对称结构。之后，经过T/4时间后，在N4时刻，涡流场又重新恢复到初始时刻N0的状态，完成一个完整的尾涡脱落过程。从以上的分析中可以看出，对于本例的基本锁定工况，在一个圆柱谐振周期内，会释放方向相反的两个涡(即正涡A和负涡B)，这是典型的2S尾涡脱落模式。值得注意的是，由于本文的圆柱运动还同时涉及到绕轴心的转动，转动运动分量的出现会使得圆柱表面出现额外的剪切作用，从而形成附加的剪切层，并影响分离点的位置。在图8中，由旋转自由度引起的附加剪切层可明显见于N2和N4时刻，其剪切作用的强弱与相位差有直接关系，并可影响到圆柱的受力以及锁定等宏观特性，这与简单的单自由度横流向受迫振动情形是不同的[24]。

图8 φ=-60°，f*=0.8时一个涡脱落周期的涡量图

图9进一步以φ=-60°和f*=1.2条件下的准周期锁定情况为例(涡脱落周期为圆柱振动周期的5倍，对应图4(c))，给出了涡流场的演化特性，其中各图之间的时间间隔为1/2圆柱的振动周期。从图9中可以看出，在初始时刻N0，圆柱的背流侧存在一个由负涡A和正涡B构成的涡对。经过半个圆柱振荡周期后，在N1时刻，圆柱右上方的负涡A分裂为一个贴体负涡A1和一个独立负涡A2。由于A2是单独脱落的负涡，在此将其定义为S-形式。此后的N2和N3阶段，圆柱尾流区内未见明显的涡旋脱落，伴随着A2不断向下游输运，A1和B涡都得到了进一步的发展。在N4时刻，即圆柱完成2个完整的振荡周期时，上述的贴体负涡A1分裂成一个新的贴体涡A11和一个独立的分离涡A12，其中后者仍以先前的S-单涡形式发生脱落。与此同时，正涡B也一次性分裂为B1、B2和B3，其中B1以贴体涡的形式存在，而B2和B3以成对的正涡形式同时脱落，本文将其定义为P+形式。应该指出的是，在N4时刻同步出现的S和P涡脱落模式是以往在仅有横向振动的情况下所未见的，这体现了旋转自由度的重要作用。在N5时刻，贴体涡B1和A11都得到了进一步的发展，同时B2与B3在对流作用下，他们之间的分离也越发明显。至N6时刻，贴体涡B1也发生分裂，形成B11和B12，其中B12是以S+的形式脱落，而B11仍附着在圆柱表面。N7和N8是尾流场的进一步发展阶段，贴体涡B11和A11得到了充足的拉伸，并具备了进一步分裂的条件。在N9时刻，再次出现了S和P的同步分裂模式。B11分裂为贴体涡B111和单体脱落涡B112，后者构成S+的脱落形式。而A11则一次性分裂为贴体涡A111和一个负的脱落涡对A112与A113，即P+的形式。此后，在经过5个完整圆柱振荡周期后，流场于N10时刻重新恢复到初始的N0时刻状态，此时的贴体涡旋B111和A111分别对应N0时刻的B涡和A涡。

图9 φ=-60°，f*=1.2时完整涡脱落周期的涡量图

纵观上述整个涡旋脱落过程，可以发现，N0与N5、N1与N6、N2与N7、N3与N8以及N4与N9，互为反对称涡量场(涡旋结构)。同时，整个尾迹区的演化过程包括了涡旋的多次二级分裂和三级分裂，这些分裂在本质上是对应着非线性问题中的分叉行为。而分叉的发生正是流动不稳定性的体现，也是层流向湍流转捩的重要形式之一[25]。

针对上述的5倍振动周期的准周期锁定情况，相应的尾迹区涡旋分裂过程总结于图10中。对于初始时刻的正涡A，先后在N1和N4时刻发生两次二级分裂，之后于N9时刻进一步发生三级分裂。对于正涡B，则是首先在N4时刻发生一次三级分裂，之后于N6和N9时刻分别发生两次二级分裂。这样，在一个尾迹演化周期内，由初始的一对涡旋，共分裂生成14个次级涡旋。而就涡旋的脱落过程而言，历经了{S-}，{S-，P+}，{S+}，{S+，P-}，共形成2正2负4个单涡以及1正1负2个涡对的脱落。本文在此仅给出了以5T为代表的准周期锁定情况下的尾流场结构演化特征，对于其它具有更低频率的准周期锁定工况，其尾涡模式也更为复杂。但就总体特征而言，都具有复杂的多级次涡旋分叉过程。

4 结 论

本文通过高阶Taylor-Galerkin迎风有限元方法求解二维不可压缩Navier-Stokes方程，并结合任意拉格朗日-欧拉动网格方法，针对圆柱在均匀流中，同时发生旋转运动及横流平动的两自由度谐振问题开展了数值模拟研究，重点讨论了圆柱谐振频率f*以及转动与平动之间的相位差φ对圆柱受力、锁定特性以及尾涡脱落模式的影响作用。本文研究工作的主要结论如下：

(1) 对于基本锁定工况，尾涡脱落频率与圆柱的振动频率一致；而对于准周期锁定的情况，完整的涡脱落过程所对应的频率为fe/N(N=3, 5，…)。此外，随着流场的演化，升力系数中还存在频率为(2M-1)fe/N的次谐波成份(M=1, 2, …)。

(2) 圆柱旋转运动与横流向振动之间的相位差φ对锁定区间有重要影响作用。当0°≤φ<180°时，锁定主要以基本锁定为主，而-180°<φ<0°时，基本锁定所对应的频率范围比较狭窄，在高频区广泛分布准周期锁定现象；当φ=±180°时，无论是基本锁定还是准周期锁定，两个相位差下的数值计算结果完全一致；而当φ=±120°及±60°时，正负相位差条件下的基本锁定区间及准周期锁定的分布均存在较大的差别，体现出流体与结构相互作用的非线性特性。

(3) 圆柱旋转运动与横流向振动之间的相位差φ对圆柱的受力也存在较大的影响作用。总体上，在各相位差条件下，升力系数的均方根均随着圆柱振动频率的增大而增大；相比于圆柱只发生横向强迫振动的情况，当引入圆柱旋转振动后，锁定区间内的升力系数均方根有不同程度的增大；平均拖曳力系数的特性随圆柱振动频率的变化则相对复杂。当φ=120°时，平均拖曳力系数最小，并且随圆柱振动频率的增大几乎不发生变化；相比于圆柱只发生横向强迫振动的情况，除φ=120°及60°，其他相位差条件下所对应的平均拖曳力系数均增大。

(4) 对于基本锁定情况，其尾迹区内为典型的2S涡脱落模式；本文以φ=-60°，f*=1.2为例，给出了准周期锁定下的涡脱落模式。数值计算结果表明，在一个完整的尾涡脱落周期内，会存在涡旋结构的多次二级分裂及三级分裂，这些分裂的过程在本质上对应着非线性问题中的Hopf分叉行为，同时这也是流动不稳定性的具体表现。

[1] WILLIAMSON C H K, ROSHKO A. Vortex formation in the wake of an oscillating cylinder[J]. Journal of Fluid and Structures, 1988, 2(4): 355-381.

[2] LU X Y, DALTON C. Calculation of the timing of vortex formation from an oscillating cylinder[J]. Journal of Fluid and Structures, 1996, 10(5): 527-541.

[3] BLACKBURN H M, HENDERSON R D. A study of two-dimensional flow past an oscillating cylinder[J]. Journal of Fluid Mechanics, 1999, 385: 255-286.

[4] ANAGNOSTOPOULOS P. Numerical study of the flow past a cylinder excited transversely to the incident stream. Part 1: Lock-in zone, hydrodynamic forces and wake geometry[J]. Journal of Fluid and Structures, 2000, 14(6): 819-851.

[5] CARBERRY J, SHERIDAN J, ROCKWELL D. Forces and wake modes of an oscillating cylinder[J]. Journal of Fluid and Structures, 2001, 15: 523-532.

[6] GUILMINEAU E, QUEUTEY P. A numerical simulation of vortex shedding from an oscillating circular cylinder[J]. Journal of Fluid and Structures, 2002, 16(6): 773-794.

[7] SARPKAYA T. A critical review of the intrinsic nature of vortex-induced vibration[J]. Journal of Fluid and Structures, 2004, 19(4): 389-447.

[8] LEONTINI J S, STEWART B E, THOMPSON M C, et al. Wake state and energy transitions of an oscillating cylinder at low Reynolds number[J]. Physics of Fluids, 2006, 18: 067101.

[9] ONGOREN A, ROCKWELL D. Flow structure from an oscillating cylinder. Part 2. Mode competition in the near wake[J]. Journal of Fluid Mechanics, 1988, 191: 225-245.

[10] CETINER O, ROCKWELL D. Streamwise oscillations of a cylinder in a steady current. Part 1. Locked-on states of vortex formation and loading[J]. Journal of Fluid Mechanics, 2001, 427: 1-28.

[11] XU S J, ZHOU Y, WANG M H. A symmetric binary-vortex street behind a longitudinally oscillating cylinder[J]. Journal of Fluid Mechanics, 2006, 556: 27-43.

[12] LEONTINI J S, JACONO D L, THOMPSON M C. Wake state and frequency selection of a streamwise oscillating cylinder[J]. Journal of Fluid Mechanics, 2013, 730: 162-192.

[13] DU L., SUN X. Suppression of vortex-induced vibration using the rotary oscillation of a cylinder[J]. Physics of Fluids, 2015, 27: 023603.

[14] BAI W. Numerical simulation of flow past a rotating and rotary oscillating circular cylinder on unstructured meshes[J]. Coupled System Mechanics, 2013(2): 191-241.

[15] MITTAL S, KUMAR B. Flow past a rotating cylinder[J]. Journal of Fluid Mechanics, 2003, 476: 303-334.

[16] LU L, QIN J M, TENG B, et al. Numerical investigations of life suppression by feedback rotary oscillation of circular cylinder at low Reynolds number[J]. Physics of Fluids, 2001, 23: 033601.

[17] BAEK S J, SUNG H J. Numerical simulation of the flow behind a rotary oscillating circular cylinder[J]. Physics of Fluids, 1998, 10(4): 869-876.

[18] CHOI S, CHOI H, KANG S. Characteristics of flow over a rotationally oscillating cylinder at low Reynolds number[J]. Physics of Fluids, 2002, 14(8): 2767-2777.

[19] BLACKBRUN H M, ELSTON J R, SHERIDAN J. Bluff-body propulsion produced by combined rotary and translation oscillation[J]. Physics of Fluids, 1999, 11(1): 4-6.

[20] NAZARINIA M, JACONO D L, THOMPSON M C, et al. Flow behind a cylinder forced by a combination of oscillatory translational and rotational motions[J]. Physics of Fluids, 2009, 21: 051701.

[21] JIANG C B, KAWAHARA M. A three-step finite element method for unsteady incompressible flows[J]. Computational Mechanics, 1993, 11: 355-370.

[22] 郭晓玲, 唐国强, 刘名名, 等. 低雷诺数下串联双圆柱涡激振动机理的数值研究[J]. 振动与冲击, 2014, 33(4): 60-69.

GUO Xiaoling, TANG Guoqiang, LIU Mingming, et al. Numerical investigation on vortex induced vibration of twin tandem circular cylinders under low Reynolds number[J]. Journal of Vibration and Shock, 2014, 33(4): 60-69.

[23] 刘为民, 谷家扬, 陶延武, 等. 低质量比圆形四立柱涡激运动特性研究[J]. 振动与冲击, 2015, 34(19): 175-180.

LIU Weimin, GU Jiayang, TAO Yanwu, et al. Vortex induced motion characteristics of four circular columns with a low mass ratio[J]. Journal of Vibration and Shock, 2015, 34 (19): 175-180.

[24] WILLIAMSON C H K, GOVARDHAN R. Vortex-induced vibration[J]. Annual Review of Fluid Mechanics, 2004, 36: 413-455.

[25] 朗道 Л Д, 栗弗席兹 E M. 流体动力学[M]. 5版. 李植 译, 北京: 高等教育出版社, 2013.

Hydrodynamic characteristics of laminar flow over a circular cylinder having forced rotational and transverse harmonic oscillations

LIU Mingming, TANG Guoqiang, LÜ Lin, TENG Bin

(State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 1160241, China)

Based on the two-dimensional finite element solution to incompressible viscous Navier-Stokes equations using the Arbitrary Lagrangian-Eulerian (ALE) method, the hydrodynamic characteristics of laminar flow over a circular cylinder having forced rotational and transverse harmonic oscillations were examined. The effects of oscillatory frequency (f*[0.5, 2.0] with an interval of 0.05) and phase difference between the rotation and translation motions (φ=±180°, ±120°, ±60° and 0°) on fluid force, lock-in mode and wake mode were analyzed. The numerical results showed that two distinct lock-in modes are identified under the present two forced oscillation excitations, they are the ordinary primary lock-in and the nonlinear quasi-periodic lock-in; the phase difference between the rotation and translation motions has an important influence on the types of lock-in and the corresponding lock-in frequency intervals, symmetric positive and negative phase differences can cause highly consistent or completely different lock-ins; fluid forces exerted on the circular cylinder generally increase with increase in oscillating frequency, the obvious hydrodynamic coefficient’s jumps are observed under the condition of lock-in; for the specific quasi-periodic lock-in cases, there are multi-stage nonlinear bifurcations of vortex structure in the field of wake flow, the wake evolution process follows quasi-periodic features.

circular cylinder; rotational oscillation; transverse oscillation; quasi-period lock-in; wake mode

国家重点基础研究发展计划(973计划)(2014CB046803);国家自然科学基金(51409035；51279029)

2016-01-12 修改稿收到日期：2016-04-12

刘名名 男，博士生，1986年9月生

吕林 男，副研究员，1976年生

P751

A

10.13465/j.cnki.jvs.2017.11.006