涡激振动是一种典型的流固耦合问题，它是由于钝体绕流产生的流体力导致的结构振动。当旋涡脱落频率接近结构的自振频率时会发生共振现象，结构振动幅度大幅增大，其长期作用会导致结构物的疲劳破坏，因此对于涡激振动及其抑制效应的研究是必不可少的。本文针对三维变截面圆柱的涡激振动进行模拟，发现通过改变结构形式，可以有效抑制变截面圆柱的涡激振动响应。

旋涡的脱落可以使圆柱产生横向及流向的振动。双自由度振动情况下，流向振动对圆柱横向振动产生的影响，是众多学者热切关注的问题。目前学者们已从物理模型和数值模拟两方面对涡激振动问题进行了深入研究。Williamson等^[1-5]对低质量阻尼比圆柱的涡激振动进行了实验研究，为以后的研究提供了重要指导。J.M. Dahl等^[6-7]通过物理模型实验对频率比为1 < f_x/f_y < 1.9的圆柱的涡激振动进行了实验研究，发现圆柱的振动轨迹形态与振动频率比有关。随后Sanchis^[8]在J.M. Dahl基础上对流横向频率比小于1的涡激振动进行了研究，发现频率比小于1时振动特性与f_x/f_y=1时基本相同。

计算机技术的发展为涡激振动问题的数值模拟提供了研究基础。ARD Silva等^[9]采用IBM法对二维单自由度圆柱涡激振动进行数值模拟，分析了涡脱形态及质量比对圆柱振幅的影响。Antoine等^[10]采用任意拉格朗日欧拉法，对Re=100情况下2维单自由度涡激振动进行数值模拟分析，发现旋涡脱落模态与圆柱振动频率有关。Li等^[11]采用有限元法及新显式积分法对Re=200情况下圆柱涡激振动情况进行分析，并对比单自由度及双自由度涡激振动特性，发现流向振动对横向振动有较大影响。Yan等^[12]对Re=150情况下双自由度单柱及串列双柱进行数值模拟分析，并考虑横向与流向不同自振频率比对圆柱涡激振动的影响。

本文将模型简化为质量-弹簧阻尼模型，采用四阶Runge-Kutta法求解运动方程，得到圆柱的涡激振动响应。在二维情况下通过改变结构的自振频率来改变约化速度，从而分析不同约化速度下双自由度圆柱涡激振动情况。在三维情况下模拟A.Sanchis^[13]的模型试验，通过对比分析进一步验证程序的准确性，并在此基础上进行变截面圆柱涡激振动模拟分析，考虑了柱体轴向截面的改变对柱体振动的影响。

1 控制方程与动网格技术 1.1 流场控制方程及参数

对于不可压缩粘性流体的控制方程为连续性方程及动量方程(N-S方程)，表示为

$ \nabla \cdot\mathit{\boldsymbol{u}} = 0 $

(1)

$ \frac{D}{{Dt}}\mathit{\boldsymbol{u}} = \frac{\partial }{{\partial t}}\mathit{\boldsymbol{u}} + {\rm{ }}(\mathit{\boldsymbol{u}}\cdot\nabla )\mathit{\boldsymbol{u}} = - \frac{1}{\rho }\nabla \mathit{\boldsymbol{P}} + \upsilon {\nabla ^2}\mathit{\boldsymbol{u}} $

(2)

式中：u、P分别代表速度和压力矢量，ρ、υ分别代表流体密度与运动粘度，对于水取ρ=1 000 kg/m³，υ=1.0×10^-6 m²/s，$ \nabla $为梯度算子。

本文采用SST k-ω湍流模型考虑湍流的影响，对于k、ω，其输运方程为

$ \frac{{\partial k}}{{\partial t}} + \frac{{\partial (k{\mathit{\boldsymbol{u}}_j})}}{{\partial {x_j}}} = {\rm{ }}{\mathit{\Gamma} _k}\frac{{{\partial ^2}k}}{{\partial {x_j}\partial {x_j}}} + {\rm{ }}{G_k} - {Y_k} + {S_k} $

(3)

$ \frac{{\partial \omega }}{{\partial t}} + \frac{{\partial (\omega {\mathit{\boldsymbol{u}}_j})}}{{\partial {\mathit{\boldsymbol{x}}_j}}} = {\rm{ }}{\mathit{\Gamma} _\omega }\frac{{{\partial ^2}\omega }}{{\partial {\mathit{\boldsymbol{x}}_j}\partial {\mathit{\boldsymbol{x}}_j}}} + {\rm{ }}{G_\omega } - {Y_\omega } + {D_\omega } + {S_\omega } $

(4)

式中：k、ω分别表示湍动能和比耗散率，Γ、G、Y、S分别表示各自的有效扩散率、生产项、耗散项和源项，D_ω表示交叉扩散项。

1.2 运动方程及参数

对于双自由度刚性圆柱，其运动方程可以表为

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

(5)

$ m\frac{{{\partial ^2}x}}{{\partial {t^2}}} + c\frac{{\partial x}}{{\partial t}} + kx = {F_D} $

(6)

式中：m、c和k分别代表结构的质量、阻尼和刚度，x、y分别为结构的流向位移、横向位移，F_D和F_L分别表示阻力和升力，通过联立式(1)、(3)、(4)沿浮体湿表面积分求出的流体压力。上述运动方程采用四阶Runge-Kutta法求解获得下一时刻物体运动的位置和速度。

1.3 动态层铺动网格技术

动态层铺模型针对单立柱涡激振动问题具有网格生成快且网格质量较好的优势。其原理是依据运动边界表面首层网格高度的变化，使网格进行分裂或者合并，从而实现网格运动。在FLUENT中，通常需要定义一个初始边界层网格高度h_i，分裂因子α及合并因子β。通过初始边界层网格高度h_i，分裂因子α及合并因子β计算得到分裂高度h_s及合并高度h_m。当边界层网格高度超过分裂高度h_s或小于h_m，边界层网格即进行分裂或者合并。

2 二维计算分析 2.1 计算模型及边界条件

将圆柱振动简化为质量-弹簧-阻尼系统，如图 1所示。网格划分采用结构化网格划分方式，在圆柱表面网格采用渐进加密方式划分。计算区域网格划分结果如图 2所示。初始状态下，一共存在63 240个网格。随着计算的进行，网格重新划分，网格数量会稍有增加。根据SST k-ω湍流模型，要求近壁面网格质量达到y+≈1。需要在粘性底层内布置足够的精细网格进行求解。本文对近壁面采用渐进加密方式划分。y+在0~1，符合网格划分要求。

边界条件设置：左端为流体入口，采用速度入口边界条件；右端为出口，采用完全发展出流条件；上下边界采用对称边界条件；圆柱体表面采用光滑无滑移壁面边界条件。

动网格采用动态层铺模型，为了减少网格变形引起的结构物表面流体计算的误差，圆柱体周围设置为随体运动区域，网格不发生变化。在10D×4D范围在流向方向与圆柱保持相同运动，保证圆柱及其随体运动区域相对于矩形区域10D×4D仅做单自由度方向运动。

流场计算采用FLUENT求解，压力与速度的耦合方式采用SIMPLEC算法，对流项采用二阶迎风格式。湍流模型选用SST k-ω湍流模型。将四阶Runge-Kutta方法代码写入用户自定义函数(UDF)求解运动方程。在求解过程的每一时间步内，首先求解流体控制方程，提取流体作用力；将流体力代入运动控制方程，采用四阶Runge-Kutta法进行求解，获得圆柱运动的速度和位移。

2.2 程序验证

为验证程序准确性，参考文献[12, 14]的两组数值模拟结果进行对比验证。选取圆柱的质量比m^*=2.0，阻尼比ζ=0，对雷诺数Re=150，约化速度在3 < U_r < 8的涡激振动进行数值模拟，振幅响应实验结果如图 3所示。

从图中可以看出本文数值模拟结果与文献[12, 14]结果吻合很好，在4 < U_r < 7圆柱达到共振状态，在U_r=4(U_r=U/(f_nwaterD), U为来流速度，f_nwater为圆柱在静水中的自振频率，D为圆柱直径)时得到圆柱最大无因次振幅A/D≈0.57。因此，本部分数值模拟可以验证程序的准确性。

2.3 计算工况

数值模拟可以通过改变圆柱的自振频率改变约化速度，即通过式U_r=U/(f_nwaterD)来改变约化速度U_r，实现在固定雷诺数下对圆柱涡激振动问题很好的研究。本文在二维模拟部分，针对固定雷诺数下，低质量比的双自由度圆柱涡激振动进行模拟分析。选取圆柱尺寸D为0.025 m，流速为0.045 m/s，即雷诺数为Re=1 125。结构的阻尼比选取为ζ=0.003，质量比为m^*=2。对约化速度2 < U_r < 9的涡激振动响应进行模拟分析。

2.4 振幅响应分析

对雷诺数Re=1 125情况下双自由度圆柱的涡激振动幅值进行数值模拟分析。图 4为圆柱涡激振动幅值的时间历程曲线。U_r=2时圆柱的横向振幅及流向振幅均为较小值，此时圆柱处于初始支阶段。在U_r=3和U_r=4时流向振幅和横向振幅均存在明显的周期性的跳跃，即发生了“拍”的现象，这是由于多个振动频率叠加的结果。并且此时圆柱的振幅峰值发生了较大改变。U_r=3时，横向振幅峰值最大值为0.31，存在着明显的“拍”现象。当U_r=4时，横向振幅峰值最大值增长到0.65，此时发生了三支响应中“初始支”到“上支”的跳跃。当U_r=5时，横向无因次振幅稳定在0.7左右，流向无因次振幅也达到了最大值0.1，此时“拍”现象也随之消失。当约化速度增长到U_r=9时，横向无因次振幅降低到0.42，也就是出现了从上支到下支的跳跃。

2.5 振动轨迹及涡脱形态

从圆柱的振动轨迹(如图 5)中可以看出，双自由度圆柱振动轨迹大致呈“8”形，这是因为流向控制频率大致为横向控制频率的2倍^[15]。当U_r=5时，圆柱振动轨迹呈现为“8”字形的下半部分，说明此时圆柱的横向振动频率与流向振动频率基本一致，即流向振动的控制频率与圆柱的自振频率接近，此时圆柱的流向振动幅度也达到了最大值0.1。当U_r=3和U_r=4时圆柱的振动轨迹呈现杂乱无章的状态，这是由于过渡区内多种振动频率叠加导致的结果，同时在振幅上表现出了“拍”的现象。当约化速度继续增大，圆柱的振动轨迹又表现为“8”字形。

图 6为不同约化速度下的涡脱模态图，从涡脱模态图中可以看出，U_r=2和U_r=3时涡脱模态表现为2S模态，圆柱两侧各有单个旋涡脱落，此时圆柱的振动幅度较小。当U_r=4时，涡脱模态呈现不规则状态，这与此时出现“拍”的现象以及圆柱不规则的振动轨迹相对应。当U_r=5和U_r=7时，旋涡脱落模态呈现2P模态，振动周期内圆柱两侧各有一对反向旋涡脱落，且在远场区域能量较小的旋涡逐渐耗散消失。当U_r继续增大时，旋涡涡脱模态又转变为2S模态，此时圆柱的振动幅值降低。

2.6 振动频率

本文通过改变圆柱的自振频率实现不同约化速度下的涡激振动模拟。根据雷诺数Re=1 125，可得对应的斯特罗哈数St≈0.21，由斯特罗哈数与式St=f_St·U/D可确定固定圆柱的涡脱频率为f_St=0.378。从图 7可以看出，随着约化速度的增大，涡脱频率与结构的自振频率逐渐接近。在5≤U_r≤9，圆柱的涡脱频率锁定到振动频率上，并接近自振频率，此时发生共振现象，圆柱的横向振幅明显增大。当U_r脱离锁定范围，圆柱的涡脱频率在斯特罗哈频率f_St左右波动，并且振动频率与涡脱频率保持一致，即此时圆柱的振动频率由斯特罗哈频率f_St控制。

3 三维计算分析 3.1 验证分析

考虑双自由度可以更真实地模拟模拟涡激振动，因此在三维情况下对双自由度刚性圆柱的涡激振动进实验。圆柱直径为D=0.08 m，长度为0.475 m，质量比m^*=1.04，阻尼比ζ=0.046，圆柱在静水中自振频率为f_nwater=0.42 Hz，附加质量系数Ca=1。模型实验中由于实验设备的限制，只进行了约化速度在4.8 < U_r < 6.9的涡激振动实验，但此振幅已足够接近三支响应振幅峰值^[13]。

数值模拟结果如图 8所示，在低约化速度下振幅模拟结果有所偏差，高约化速度情况下吻合较好。由于低约化速度情况下振幅较小，其模拟结果稍有偏差但在可接受范围内。

3.2 变截面圆柱振幅响应分析

在A.Sanchis^[13]实验基础上，进行变截面刚性圆柱的涡激振动数值模拟，分析变截面圆柱的涡激振动情况。变截面圆柱的几何模型及柱体表面网格划分如图 9所示, 圆柱的大直径尺寸为小直径尺寸2倍，圆柱的小直径及长度与实验参数一致，即小圆柱直径为D=0.08 m，大圆柱直径为D=0.16 m圆柱总长度为0.475 m。柱体表面网格划分依旧按照y+≈1进行网格划分。y+在0~3符合网格划分要求。变截面圆柱的质量按照与原模型相同质量比进行计算得到。变截面圆柱的阻尼比ζ及自振频率等按原参数选取。

图 10给出部分约化速度下变截面圆柱的横向振幅的时间历程曲线。约化速度为U_r=5时变截面圆柱的横向无因次振幅约为0.1。在U_r=6时圆柱的横向振幅达到最大值约为0.35，但此时变截面圆柱的横向振幅也小于传统圆柱的横向振幅。而后随着约化速度的增大，圆柱的无因次振幅逐渐较小，在U_r=9时，圆柱的振幅约为0.175，当U_r=12时，圆柱的振幅降低到0.125左右。

将变截面圆柱振幅与典型圆柱振幅进行对比，结果如图 11所示。变截面圆柱的振幅响应在U_r=6时达到最大值，而后变截面圆柱的振幅又逐渐降低。变截面圆柱的振幅响应远小于原始圆柱的振幅响应。由于本文选取与原始圆柱相同质量比，此类现象的产生应该与变截面处的突变有关。为讨论变截面圆柱的运动响应，还需从涡脱形态及振动频率等方面进行讨论分析。

3.3 变截面圆柱涡脱形态

图 12给出部分约化速度下圆柱的涡脱形态图。从图中可以看出，在以下约化速度中圆柱的涡脱形态均为“2S”形态，这与较低的振动幅度相对应。在变截面处，旋涡脱落形成的涡柱存在明显的被“割裂”的现象，这与圆柱的振动幅度降低及“2S”涡脱形态有密不可分的联系。此外这种现象的产生，应该与变截面圆柱的覆盖率及变截面部分圆柱的直径与原始圆柱的直径的比值有关。

3.4 变截面圆柱振动频率分析

对以上约化速度下变截面圆柱的振动响应进行频域分析，得到变截面圆柱的振动频率如图 13所示。U_r=5时，变截面圆柱的横向振动频率集中在0.406 9 Hz, 变截面圆柱在静水中的自振频率为0.42 Hz，此时圆柱的振动频率比较接近圆柱的自振频率，但此时圆柱的振动幅度较低，可能是由于在较低的约化速度下，涡柱被割裂对圆柱的振动响应影响较大。在U_r=6和U_r=9时，变截面圆柱的振动频率均锁定在0.447 6 Hz，但由于U_r =9时频谱峰值较低，其振幅响应较低。随着约化速度的继续增大，变截面圆柱的振动频率继续增大，圆柱脱离锁定区域，U_r=12时圆柱的振动频率为0.488 3 Hz，其振幅降低到0.13左右。

4 结论

1) 在Re=1 125情况下，二维圆柱涡激振动在“初始支”到“上支”的过渡阶段，出现了“拍”的现象，并且伴随着振动轨迹不规则和涡脱形态混乱的现象。

2) 当约化速度U_r=5时，流向振动频率与横向振动频率相一致，圆柱的流向振幅达到最大值，并且此时圆柱的振动轨迹表现为“8”字形的下半部分。

3) 在三维情况下，对变截面刚性圆柱涡激振动进行数值模拟，并从振幅、涡脱形态及振动频率进行分析。结果显示，变截面圆柱的振动幅度较原始圆柱的振动幅度有所降低。变截面圆柱的突变处会对涡柱产生“割裂”的影响，这与变截面圆柱的振动响应降低有不可分割的联系。