考虑旋转自由度的涡激振动数值模拟研究

康 庄， 倪问池

(哈尔滨工程大学 船舶工程学院， 哈尔滨 150001)

涡激振动现象在海洋工程中广泛存在，对涡激振动研究已成为研究的热点之一。在海洋工程中，柱形结构被广泛采用，且质量比都较低[1]。因此，对涡激振动的研究大多集中在低质量比圆柱结构上。

Khalak等[2]对弹性固定的圆柱结构进行了模型实验研究，发现当圆柱质量阻尼比较高时，涡激振动响应存在两个分支，分别对应于初始分支与下端分支，而当质量阻尼比较低时，涡激振动响应存在三个分支，分别对应于初始分支、上端分支和下端分支。Jauvtis等[3]还对m*=2.6的低质量阻尼比圆柱进行了双自由度涡激振动实验研究，发现在加速工况中，圆柱横向最大振幅可达1.5D左右，并且存在回滞现象。此外，还将圆柱涡激振动的尾涡形式总结为2S、2P、2T和2C模式。Khalak等[4]通过对双自由度圆柱的实验研究指出，当质量比较大时，圆柱共振频率将约等于固有频率。而当质量比较小时，圆柱共振频率将会略大于固有频率。进一步，Govardhan等[5]发现，当质量比小于0.54时，圆柱涡激振动响应的下端分支将会缺失。

根据现有的研究成果，当前涡激振动的研究对象主要为单自由度或双自由度弹性固定的圆柱，对于增加了旋转自由度的涡激振动，主要的研究形式为人为对结构物施加固定的旋转角速度，分析流场的特性。

Naik等[6]对雷诺数等于100，不同旋转角速度下的椭圆结构的流场特性进行了数值模拟研究，发现了旋转运动对流场特征产生了显著影响。在旋转效应下，尾涡模式发现了变化，观察到了“盘漩涡”以及多频现象。Lam[7]对旋转圆柱尾流场进行了实验研究，并发现当圆柱无量纲旋转角速度小于1.9时，圆柱的涡泄模式与固定圆柱类似，但随着旋转速度的增加，圆柱的涡泄模式会发生显著变化，并且偏向一边。

已有文献表明，当前，在对旋转结构尾流场的研究中，已经发现了一些新颖的现象。然而这些研究中，结构都是在人为作用下发生旋转，对于结构在涡致力作用下发生的旋转，已有的研究并不多。

Wilde[8]对自由站立式立管进行了涡激运动模型试验，发现在一定的来流速度下，浮筒会产生艏摇运动。李鹏[9]在对浮力筒的六自由度涡激运动实验中也发现了艏摇现象。

涡激振动的艏摇现象已经在浮力筒的实验与工程实际中被发现。这表明尾涡的周期性脱落，还会产生周期性的力矩，对于浮力筒这类系泊结构，将会造成艏摇现象，对于弹性或刚性固定的立管，则会造成周期性的扭矩，这将大大降低结构的疲劳寿命，甚至有可能使结构发生扭转变形。然而，当前对于该现象缺乏详细的分析。本文基于本文基于OpenFoam软件，对考虑旋转自由度下的低质量比圆柱涡激振动响应特性进行数值模拟研究。首先，参照Jauvtis等的实验，对双自由度低质量比圆柱进行数值模拟，然后在相同的边界条件下，对圆柱增加旋转自由度，进行数值模拟，并对两种自由度下涡激振动响应进行对比分析。

1 数值方法

1.1 湍流模型

本文数值模拟中，流体采用的是雷诺平均控制方程，其表达式为

(1)

(2)

式中：u,p为时均速度与压力；μ为动力黏度；ρ为流体密度。

为了对上述控制方程进行封闭，结合了k-ω与k-ε湍流模型优势的SST模型[10]被应用于本次数值模拟中。SST模型具有较高计算精度以及稳定性，在工程中被广泛应用。其无量纲的形式如下

涡黏性

(3)

湍动能

(4)

比耗散率

(5)

第二混合函数

(6)

湍动能k的生成项

(7)

第一混合函数

(8)

其中，

(9)

F1在近壁面处取1，在远壁面处取值为0，实现了在近壁面处采用k-ω模型，在远壁面处使用k-ε模型。相应地，输运方程中的系数可表示为

φ=φ1F1+φ2(1-F1)

(10)

SST模型中，各系数的取值如表1所示。

表1 SST湍流模型系数

固体可简化为质量-弹簧-阻尼系统，其运动的控制方程可表示为

Mx″(t)+Kxx(t)+Cxx′(t)=Fx(t)

My″(t)+Kyy(t)+Cyy′(t)=Fy(t)

式中：M为圆柱质量；K和C分别为弹簧刚度和阻尼系数；x(t)、y(t)与θ(t)分别为顺流向位移、横流向位移以及转角位移；F为流体力；Iθ为圆柱对瞬心的转动惯量；r为位置矢量。

1.2 计算域及网格

本次数值模中，圆柱参数参照Jauvtis等的实验设定参数，具体参数为：直径D=0.038 1 m，质量比m*=2.6,质量阻尼比m*ξ=0.013，静水中振动固有频率fn=0.4。流场的网格划分如图1所示。考虑到圆柱振动区域以及尾涡的影响范围，总体计算域范围取为：-12D

时间步长根据库朗数(Cn)确定[12]。库朗数的表达形式为

(11)

式中：u为流速；Δt为时间步长；Δx为网格尺寸。在不同工况中，时间步长根据流速作相应调整，使得计算过程中库朗数最大值为0.2左右。

图1 计算域整体网格

本文采用的网格以及数值模拟方法与文献[13]一致。文献[13]的计算结果表明，该网格及数值模拟方法对该工况下双自由度涡激振动的计算具有较高的精度。考虑到增加旋转自由度后，流场特征并未发生显著变化，因此，该网格对于增加旋转自由度后的涡激振动数值模拟仍然具有较高的正确性。

1.3 旋转自由度的设定

本文数值模拟中，对涡激振动艏摇现象的考察，通过对相应的OpenFoam限制文件中，增加圆柱的旋转自由度实现。运用后处理软件，对圆柱涡激振动稳定阶段数值模拟结果进行分析，可提取旋转角速度时历数据，圆柱的艏摇角度可通过对角速度进行积分求得。此外，圆柱的旋转阻尼设为0。

2 计算结果及分析

2.1 计算工况

根据Khalak等的实验结果，涡激振动响应在流场加速、匀速以及减速工况下会呈现出不同的响应。本文仅分析流场加速的工况，即流速从0匀速增加到目标值，然后保持匀速。

运用上述网格，分别计算稳定阶段流场约化速度范围从0～14，双自由度，以及增加旋转自由度的涡激振动响应，并进行对比。

2.2 振幅响应对比

两种自由度下，涡激振动随约化速度的横向及顺流向最大振幅响应如图2和图3所示。此外，Jauvtis等的实验结果也绘于图中作为参照。

图2 不同约化速度下横向最大振幅响应

图3 不同约化速度下顺流向最大振幅响应

如图所示，横流向以及顺流向的数值模拟的结果均成功捕捉到三个响应分支，分别为初始分支，上端分支与下端分支。上端分支与下端分支的跳跃点在约化速度Ur=7附近。数值模拟的结果与实验结果基本吻合，除了数值模拟响应的上端分支与下端分支的跳跃点略微向左偏移。在林琳等[14-16]的文章中也出现了类似的现象，这也许是实际实验过程中三维效应或者其他干扰因素造成的。

此外，还可以发现，无论是横流向还是顺流向，增加了旋转自由度的涡激振动振幅响应与双自由度涡激振动振幅响应具有相同的趋势，且振动幅值较为接近，但还是可以观察到，增加了旋转自由度后，涡激振动最大横向振幅响应幅值小于双自由度涡激振动幅值。根据机械能守恒原理，增加了旋转自由度后，有一部分振动的机械能转换为旋转机械能，使得最大振幅有所减小。

为了更直观地体现两种自由度下，涡激振动振幅响应差异，将两种自由度工况的振幅差异绘于图4中。其中，横坐标代表稳定阶段的约化速度，纵坐标代表最大振幅差异，用百分比表示(增加旋转自由度后最大振幅减小的百分比)。考虑到顺流向的振幅较小，工程中主要关心的是横流向的振动，并且横流向与纵流向的振幅响应变化趋势基本一致，因此本文仅对横向振动进行分析。

图4 增加旋转自由度后振幅减小比例

从图4中可以看出，两种自由度工况的振幅差异随着流速的变化，也呈现出三段分支，且三段分支对应的流速区间与振幅响应的分支区间基本一致。在此，为了便于对照，将振幅差异图中的分支也称为“初始分支”、“上端分支”与“下端分支”。在上端分支中，振幅的减小比例约为5%，在上端分支中，振幅的减小比例约为2%，而在下端分支中，减小比例达到了8%左右。

根据机械能守恒原理，在下端分支中，转换为旋转机械能所占的比例最大，而在上端分支中，比例最小。这显然与各分支中的流场特性密切相关。为了分析其机理，将各分支中的尾涡模式也绘于图4中。

在初始分支中，尾涡的模式为“2S模式”，即在每个周期内，两个独立的，旋转方向相反的尾涡均匀形成并从物体后方交替脱落。此时，圆柱上下侧交替产生均匀的压力差，在压力差的作用下，圆柱产生位移，并出现艏摇现象。在上端分支中，尾涡模式为“2T模式[17]”，此时，物体后方每个周期发放出两个涡对，每个涡对包含三个涡，并且三个涡的旋转方向不完全相同。在“2T模式”中，上下侧存在涡同时脱落的情况，此时上下侧压力差存在部分抵消的情况，转化为旋转机械能的比例也有所降低。在下端分支中，尾涡模式为“2P模式”，在每个周期内，物体后方发放出两个涡对，并且涡对中的两个漩涡旋转方向相反。该尾涡模式的上下侧压力变化情况与“2S模式”类似，交替产生均匀的压力差，但由于每次发放出两个涡对，相对压力差较大。此时，转化为旋转机械能的比例也有所增加。由此可以断定，艏摇运动是上下侧压力差造成的。

2.3 频率响应对比

两种自由度下，涡激振动频率响应随约化速度的变化趋势如图5所示。

图5 不同约化速度下横向振动频率响应

在圆柱涡激振动频率响应中，也可以观察到明显的三段分支，分别对应于初始分支、上端分支以及下端分支。在初始分支中，振动频率随着速度的增大而增大，并且接近fst，但比fst略小。在上端分支中，fy/fn约为0.9，即横向振动频率锁定于系统固有频率。而在下端分支中，仍然可以观察到锁定现象，但此时fy的频率不再等于固有频率，fy/fn的值稳定在1.25左右。数值模拟结果与实验结果基本吻合。

此外，还可以发现，两种自由度工况下，圆柱涡激振动频率响应结果完全重合。因此可以得出，增加旋转自由度不会改变圆柱振动频率。

2.4 艏摇响应分析

下面单独考察增加了旋转自由度的圆柱涡激振动数值模拟结果，提取圆柱的最大旋转角度，分析其艏摇响应特性。各约化速度下，圆柱艏摇幅度如图6所示。其中，横坐标表示流场约化速度，纵坐标表示艏摇幅值，在图中用角度表示。

从图6中可以看出，随着流场速度的增加，艏摇幅值经历了“增加→减小→增加→稳定→减小”的过程。当Ur=2～3时，艏摇幅值随速度增加而增大，最大可达3.2度左右。当Ur=3～4时，艏摇幅值迅速减小，并达到最小值2度左右。随后，艏摇幅度随约化速度增大而增大，当Ur=6时，艏摇幅值达到最大值，5.2度左右，随后在Ur=6～10内，艏摇幅值基本保持恒定。当约化速度大于10后，艏摇幅值开始减小。

此外，如图4所示，当约化速度等于10左右时，振幅衰减量远大于约化速度为6时的值，即总机械能中，转化为旋转机械能的比例远比约化速度为6时大。但Ur=10时的艏摇幅值仍然小于Ur=6时的幅值。总体上看，艏摇幅值的变化趋势与图2所示的振幅随约化速度的变化趋势基本一致，仍然呈现出三个分支，但上端分支与下端分支之间的跳跃没有振幅响应这么明显。由此可以得出，涡激振动艏摇运动幅值大小主要取决于振动的剧烈程度，即系统总机械能的大小，但也会受到尾涡形态等流场特性的影响。

图6 不同约化速度下艏摇幅值响应

艏摇运动的频率响应如图7所示。其中，fyaw/fn表示艏摇频率与固有频率的比值，fyaw/fy表示艏摇频率与横向振动频率的比值。特别需要说明的是，当艏摇运动出现多个主频率时，取最大的频率分量。

图7 不同约化速度下艏摇频率响应

从图7中可以看出，艏摇运动的频率响应也存在明显的三个分支。从fyaw/fy的响应曲线可以看出，艏摇频率与圆柱涡激振动横向振动频率基本相等，除了在每段分支之间的过渡阶段存在略微偏差。这也许是由于在过渡阶段中，频率较为杂乱，频率分析较为困难，存在一定的误差。因此，从频率分析中可以得出，圆柱涡激振动引起的艏摇频率与横向振动频率一致，可以推断出，两者具有相同的激励源，即艏摇运动也是由于漩涡交替脱落引起的压力差造成的。

3 结 论

本文基于OpenFoam软件，采用数值模拟方法对对考虑旋转自由度的低质量比圆柱涡激振动响应特性进行了数值模拟研究，可以得出以下结论：

(1) 当具有旋转自由度时，圆柱在涡激振动过程中，会出现艏摇现象。

(2) 圆柱产生艏摇运动的原因是部分振动机械能转化为旋转机械能，因此，艏摇运动的产生会使圆柱涡激振动的幅值略微降低，且能量的转化比例与尾涡特征密切相关。

(3) 艏摇运动的出现，不会改变圆柱涡激振动的振动频率。

(4) 艏摇运动幅值的变化趋势与横向振幅的变化趋势类似，且大小主要取决于振动的剧烈程度，即系统总机械能的大小，但也会受到尾涡形态等流场特性的影响。

(5) 圆柱艏摇频率与横向振动频率一致，可以推断出，两者具有相同的激励源，即艏摇运动也是由于漩涡交替脱落引起的压力差造成的。

参 考 文 献

[1] KANG Zhuang, NI Wenchi, SUN Liping. An experimental investigation of two-degrees-of-freedom VIV trajectories of a cylinder at different scales and natural frequency ratios[J]. Ocean Engineering,2016,126:187-202.

[2] KHALAK A, WILLIAMSON C H K. Dynamics of a hydro elastic cylinder with very low mass and damping[J]. Journal of Fluids and Structures, 1996, 10(5): 455-472.

[3] JAUVTIS N, WILLIAMSON C H K. The effect of two degrees of freedom on vortex-induced vibration at low mass and damping[J]. Journal of Fluid Mechanics, 2004, 509(509):23-62.

[4] KHALAK A, WILLIAMSON C H K. Motions, forces and mode transitions in vortex-vibrations at lowmass-damping[J]. Journal of Fluids and Structures, 1999, 13: 813-851.

[5] GOVARDHAN R N, WILLIAMSON C H K. Vortex-induced vibrations of a sphere[J]. Journal of Fluid Mechanics, 2005, 531(531):11-47.

[6] NAIK S N, VENGADESAN S, PRAKASH K A. Numerical study of fluid flow past a rotating elliptic cylinder[J]. Journal of Fluids & Structures, 2017, 68:15-31.

[7] LAM K M. Vortex shedding flow behind a slowly rotating circular cylinder[J]. Journal of Fluids & Structures, 2009, 25(2):245-262.

[8] WILDE J D. Model tests on the vortex induced motions of the air can of a free standing riser system in current[C]∥ Proceedings of the Deep Offshore Technology Conference. Stavanger：DOTC， 2007.

[9] 李鹏. 浮力筒自由涡激运动试验研究[D]. 哈尔滨:哈尔滨工程大学, 2015.

[10] MENTER F R. Two-equation eddy-viscosity turbulence models for engineering applications[J]. Aiaa Journal, 2012, 32(8):1598-1605.

[11] CLOETE S, JOHANSEN S T, AMINI S. Grid independence behaviour of fluidized bed reactor simulations using the two fluid model: detailed parametric study[J]. Powder Technology, 2015, 289:65-70.

[12] MORAGUES GINARD M, BERNARDINO G, VZQUEZ M, et al. Fourier stability analysis and local courant number of the preconditioned variational multiscale stabilization (P-VMS) for euler compressible flow[J]. Computer Methods in Applied Mechanics & Engineering, 2015, 301:28-51.

[13] KANG Z, NI W, SUN L. A numerical investigation on capturing the maximum transverse amplitude in vortex induced vibration for low mass ratio[J]. Marine Structures, 2017, 52:94-107.

[14] 林琳, 王言英. 不同湍流模型下圆柱涡激振动的计算比较[J]. 船舶力学, 2013,17(10):1115-1125.

LIN Lin, WANG Yanying. Comparison between different turbulence models on vortex induced vibration of circular cylinder[J]. Journal of Ship Mechanics, 2013, 17(10):1115-1125.

[15] 谷家扬, 杨建民, 肖龙飞. 两种典型立柱截面涡激运动的分析研究[J]. 船舶力学, 2014，18(10):1184-1194.

GU Jiayang, YANG Jianmin, XIAO Longfei. Study on vortex induced motion of two typical different cross-section columns[J]. Journal of Ship Mechanics, 2014,18(10):1184-1194.

[16] PAN Z Y, CUI W C, MIAO Q M. Numerical simulation of vortex-induced vibration of a circular cylinder at low mass-damping using RANS code[J]. Journal of Fluids & Structures, 2007, 23(1):23-37.

[17] WILLIAMSON C. A high-amplitude 2T mode of vortex-induced vibration for a light body in XY motion[J]. European Journal of Mechanics-B/Fluids, 2004, 23(1):107-114.