考虑迟滞效应的圆柱体涡激振动分析研究

康 庄， 张立健， 倪问池

(哈尔滨工程大学 深海工程技术中心， 哈尔滨 150001)

在深海工程领域，圆柱体结构例如海洋立管是工程上应用最多最广泛的一种结构物，当水流绕过圆柱体时，会发生结构物后漩涡脱落的现象，交替的漩涡会激发与来流方向垂直的流体作用力，致使结构物产生振动，当泄涡频率与固有频率接近时，将会发生“共振”现象，同时圆柱振动幅值将增大[1]。圆柱体涡激振动(Vortex-Induced Vibration, VIV)一直都是研究所关注的热点问题。国内外众多学者进行了弹性支撑圆柱涡激振动研究。黄维平等[2-4]基于流固耦合的涡激升力模型和Morison方程，提出了一个非线性的圆柱体涡激振动时域分析模型，并且采用ANSYS-CFX 软件对细长柔性圆柱体双自由度涡激振动进行了深入研究。Jauvtis等[5]通过实验发现，当质量比为7时，单自由度与双自由度涡激振动试验结果类似，其横向最大振幅为1.05D附近，当质量比为2.6时，双自由度横向振幅出现“超上端分支”，其幅值达到1.5D，弹性支撑双自由度涡激振动出现如此大横向振幅，此现象是以前从未发现的。谷家扬等[6]希望采用数值模拟方法研究“超上端分支”现象，然而即使将无因次阻尼比ξ设为0，模拟出的最大横向振幅也不超过1D。部分学者注意到, 试验结果或者计算结果跟初始条件有很大的关系，Guilmineau等[7]参照Jauvtis等的单自由度低质量比试验，计算了2D横向弹性支撑圆柱(m*=2.4,m*ξ=0.013)的涡激振动,数值模拟采用了三种初始条件:①速度保持均匀不变;②速度由最高逐渐减小;③速度由最低逐渐增加。计算结果表明：圆柱横向振幅与初始条件有关:仅在①和②条件下预测到了“下支”响应;在条件③时,观察到了“上支”。试验采用了匀加速和匀减速的形式，从而产生了“迟滞现象”，而数值模拟一般采用固定速度形式，迟滞现象对涡激振动响应有显著影响。

CFD(Computational Fluid Dynamics)模拟采用固定速度初始条件时，会忽略尾流模态的“惯性”对流场影响[8]。最近，Kim等[9-12]采用数值模拟方法对圆柱涡激振动开展研究，更深入阐述涡激振动的机理。本文在前人的基础上利用开源软件OpenFOAM采用匀加速，匀减速和匀速三种初始条件对质量比为2.6双自由圆柱涡激振动模拟。研究考虑迟滞效应与不考虑迟滞时，圆柱振动幅值，锁定区间，升力系数，阻力系数以及尾涡的变化。同时针对迟滞机理进行讨论，并研究了加速度的大小对迟滞现象的影响。

1 计算模型与模拟工况

1.1 计算模型建立

本文基于OpenFoam嵌入自编程序进行圆柱涡激振动数值模拟，计算域描述如图1所示。

图1 计算域几何示意图Fig.1 Digram of computational domain

计算区域为左边及上下边距圆心6D，右边距离圆心20D的矩形流域，采用URANS二维模型，在z向设置为“empty”边界条件。网格划分必须确保“低雷诺数”处内的网格质量，“低雷诺数”并不是指整体雷诺数较低，而是指黏性子层处的雷诺数较低，为了描述黏性子层内的流动引入无量纲数y+，y+与壁面剪切应力有关。y+的计算公式由式(1)得出

(1)

式中:τw为壁面剪切力。

网格划分必须确保第一层网格厚度能使得y+<1。为了保证圆柱面处网格的厚度同时保证整体网格不过多导致计算时间冗长，因此在网格划分时，需要将圆柱周围局部加密，将计算域划分为9个区域，圆柱周围4D的范围内采用O型网格。远场区域采用结构网格，其尺寸逐步递增，如图2 所示。

采用SST(Shear Stress Transport)k-ω湍流模型和非稳态一阶隐式进行求解，动量方程的压力速度耦合采用PIMPLE算法，动量、湍流动能、耗散率项均采用二阶迎风格式。

(a) (b)图2 网格计算区域划分示意图Fig.2 Grid computing region division

1.2 时间步长控制以及网格无关性

时间步长影响着收敛速度和稳定性。合理的调整时间尺度，能够在获得理想收敛速度和稳定性的同时缩短计算时间。库朗数(Courant Number)指的是时间步长和物理空间尺度的对应关系，对于湍流而言有着重要影响，库朗数由式(2)表示

Cn=u·Δt/Δx

(2)

式中:Cn为库朗数;u为来流速度;Δt为时间步长;Δx为网格尺寸。为满足Courant-Friedrichs-Lewy条件,必须要求库朗数足够小。使用数值模拟计算时候，设定库朗数最大值小于0.50，时间步长取为Δt=0.000 1。

通过Re=104的圆柱绕流为例验证网格无关性，采用第一层网格厚度为0.02 mm，近壁面第一层网格的长厚比设置为8∶1。根据网格加密与放疏原则，设置了不同的五种网格方案，不同网格方案的详细参数如表1所示。

表1 不同网格设计方案

将阻力与斯特劳哈尔数的结果随网格总体数量的关系绘制在图3和图4中，由图3和图4可知，当网格数达到以及超过51 502时，计算的阻力系数与斯特劳哈尔数均达到了比较好的收敛。而网格数低于51 502时，阻力与斯特劳哈尔数偏差较大，由图3可以知：网格3已经能够满足计算精度的需求。

图3 阻力均值随网格数量变化趋势Fig.3 The mean drag coefficient with different grid

图4 斯特劳哈尔数随网格数量变化趋势Fig.4 The mean drag coefficient with different grid

1.3 模拟工况

选取Javutis等(2004)所做的经典的m*=2.6的双自由度实验作为参考，圆柱直径为0.038 1 m，质量阻力比m*ξ=0.13，固有频率为0.4 Hz，流向与横向频率比fnx/fny=1,约化速度范围为2～14，约化速度间隔为0.2。

本文对速度的控制是在OpenFOAM的速度自定义函数中实现的，在数值模拟过程中，三种约化速度的实现方式如下定义：

(1)匀加速。约化速度从2逐步增加，加速度为每无量纲单元时间(Ut/D)约化速度增加0.025，增加至目标速度后保持该约化速速继续模拟30个周期，同时上一个折合流速模拟结束时刻的物理场是下一个折合流速模拟的初始条件。

(2)匀减速。约化速度从14逐步减少，加速度为每无量纲单元时间(Ut/D)约化速度减小0.025，减小至目标速度后保持该约化速速继续模拟30个周期，同时上一个折合流速模拟结束时刻的物理场是下一个折合流速模拟的初始条件。

(3)匀速。在数值模拟开始时，约化速度直接给定为目标约化速度，随后保持该约化速度继续模拟。

本文对速度的控制是在OpenFOAM的速度自定义函数中实现的，三种初始条件下圆柱涡激振动各工况参数如表2所示。

表2 三种初始条件下数值模拟工况

2 计算结果与分析

2.1 涡激振动幅值响应

图3给出了不同初始条件下双自由度圆柱涡激振动横向振幅随约化速度变化曲线，方块代表Jauvtis等的试验结果，圆点代表采用匀加速方法，正三角代表采用匀减速方法，倒三角代表采用匀速方法。由图5可知，三种方式均捕获到初始分支和下端分支，采用匀加速方式时，在Ur(约化速度)=6.8时，Ay/D达到最大值1.4D。在Jauvtis等的m*=2.6试验中，Jauvtis等发现在约化速度为8.3，横向达到最大振幅比1.5D，他定义该分支为“超长端分支”，因此可知匀加速数值模拟可以捕获到“超上端分支”现象，但是其值较为提前并且结果略小。当Ur>6.8时，振幅下降至0.75D步入下端分支。采用匀减速方式时，捕获到上端分支，上端分支与下端分支交界点为Ur=6时，其最大横向振幅为0.96D，当约化速度大于6时，振幅比Ay/D陡然下降到0.7D附近，进入“下端分支”。采用固定速度方式时，没有捕获到上端分支，初始分支与下端分支交界点为Ur=5，其最大振幅0.85D。在下端分支中，匀加速的振幅略大于试验而匀减速和匀速的振幅略小于实验值。采用匀加速与匀减速初始条件可以数值模拟出迟滞现象，如图5中箭头所示。然而值得注意的是，Jauvits等(2004)试验的迟滞区域为7.8～8.6，而数值模拟得到的迟滞区域为6～6.8。数值模拟的迟滞区域相对于试验提前了2左右，但是迟滞区域的范围大致都为0.8个约化速度。

图5 不同初始条件下圆柱横向最大振幅Fig.5 Transverse max amplitude in different initial conditions

图6给出了不同初始条件下圆柱涡激振动流向振幅随约化速度变化曲线。由图6可知，采用匀加速方式时，在Ur=6.8时，Ax/D达到最大值0.29D。Jauvtis等的试验中在约化速度为7.3时达到中最大流向振幅比为0.32D，数值模拟结果较为提前并且结果偏小。当Ur>6.8时，振幅下降至0.07D步入下端分支。采用匀减速方式时，在Ur=6时，Ax/D达到最大值0.11D，数值模拟得到的迟滞区域为6～6.8。

对比图5与图6可知，当Ur≤4时，涡激振动处于初始分支中，圆柱体两侧漩涡对称地脱落，产生的升力互相抵消一部分，此时横流向运动杂乱且运动幅值较小，而顺流向运动规则且幅值较大，轨迹呈扁平状。在初始分支内，三种速度方式的流向振幅均保持较高的振动幅值。当Ur>4时，涡激振动进入锁定区域，圆柱形圆柱体两侧漩涡开始交替脱落，此时，横向作用力使横向运动幅值变大，同时顺流向运动幅值也增大。进入锁定区后，采用匀加速与匀减速的流向振幅随Ur的增大而逐步增大至极值。采用匀速初始条件时，由于流向振幅进入锁定区后没有捕获到上端分支，因此流向振幅表现为“双峰”现象，在初始分支Ur=3时达到极值0.14D，在Ur=8时达到次峰0.075D。

图6 不同初始条件下圆柱涡激流向振幅Fig.6 In-line max amplitude in different initial conditions

2.2 无因次约化振动频率响应

图7给出了圆柱涡激振动横向约化频率随约化速度变化曲线，在图中添加两条辅助线，实线为运动频率等于自振动频率辅助线，虚线为固定圆柱卸涡频率预估线。由图7可知，采用匀加速条件时：圆柱无无因次锁定频率在1.3附近。采用匀减速和匀速条件的锁定频率在1.25附近，数值模拟的锁定频率均小于实验值的1.35。三种初始条件下的锁定区域有很大的不同。采用匀加速条件时，涡激振动锁定区范围为4.2～12.8，其中处于上端分支的为4.2～6.8，处于下端分支的为6.8～12.8；当采用匀减速条件时，涡激振动锁定区范围也为4.2～12.8，但处于上端分支的为4.2～6，处于下端分支的为6～12.8；采用匀速条件时，涡激振动锁定范围为3.8～14，其中处于上端分支的为3.8～5，处于下端分支的为5～14。在锁定区间内圆柱保持高幅值振动，在锁定区外，振动频率离开实线向虚线靠拢，表现为固定圆柱泄涡频率。采用匀加速与匀减速可以得到锁定区迟滞区为6～6.8，如图7中箭头所示，其相对于试验迟滞区7.8～8.6提前约2个约化速度。

图7 不同初始条件下圆柱横向振动约化频率对比图Fig.7 Transverse reduced frequency in three initial conditions

采用匀加速初始条件时，圆柱振动频率功率函数密度分布如图8所示，当Ur=3时，横向振动频率等于固有圆柱泄涡频率，流向振动具有两个主要频率，较大振动频率对应于两倍横向振动频率，较小频率对应单倍横向频率，观察不同约化速度下振动频率可知，当约化速度较低时(Ur≤4.2)，流向振动出现了两个主要频率，分别为单倍与双倍横向频率。当约化速度Ur=6时，圆柱振动处于锁定区内，圆柱横向振动频率集中在一个频率上即圆柱固有频率，值得注意的是，在锁定区域内，发现在除了主频外，还存在一个3倍的低频分量。Jauvtis等(2004)认为，这是由于超上端分支产生的“2T”泄涡模型从而导致产生三阶升力分量。流向振动频率为横向频率2倍，同时流向频率还存在一个2倍的低频分量。当Ur增大到12.6时，圆柱处于锁定区向解锁区跳跃的过渡区，圆柱横向振动有两个主要频率，前一个频率对应固有频率，后一频率对应于固定圆柱泄涡频率。流向主要频率为横向频率的2倍，同时出现1.5倍和2倍的低频分量。

在锁定区域内，横向振动除了主频外，还存在一个3倍的低频分量，说明在锁定区域内升力存在3倍低频分量，图9给出了Ur=6时涡激升力功率谱密度曲线。由图9可知：采用匀加速条件下，锁定区域内涡激升力除主频外，还存在3倍低频以及微弱5倍低频。在锁定区域内，横向振动受3倍低频升力影响从而产生3倍振动分量。

(a)Ur=3 (b)Ur=6 (c)Ur=12.6图8 各约化速度下横向振动功率谱密度曲线Fig.8 Transverse vibration frequency PSD in different initial conditions

图9 Ur=6时涡激升力功率谱密度曲线Fig.9 Lifting force frequency PSD in increase initial conditions

锁定区域升力3倍低频分量可由泄涡模式所解释，在锁定区域内泄涡模型为“2T”模式，如图8所示。每个“T”为一个正负交替的三重涡，一个“T”模式可以看为一对半的“2S”泄涡，如图10圆圈中为一对“2S”泄涡，认为“2T”泄涡模式产生了升力中的3倍低频。

图10 “T”模式类似于一对半的“2S”模式Fig.10 “T” mode is similar to one and a half of the “2S” mode

2.3 流体力特性响应

匀加速，匀减速以及匀速三种不同速度形式所得到的升力均方根和阻力均值随约化速度变化规律如图11所示，对于升力系数而言，采用匀加速方式时，在2

图11 三种速度形式下升力均方根与阻力均值随约化速度变化图Fig.11 The RMS lift and mean drag coefficients versus Ur in different initial conditions

对比三种初始条件下的升力系数与阻力系数结果可知，作用在圆柱表面的水动力系数也存在与圆柱振幅响应相类似的迟滞区域，其区间范围为6≤Ur≤6.8。在迟滞区间内，采用匀加速条件的升力系数较大，其对应的横向振幅也较大。反之，采用匀减速条件的升力系数较小，其对应的横向振幅也较小。

图12给出了采用匀加速初始条件时无量纲位移y/D升力系数Cl以及阻力系数Cd随时间变化曲线，当Ur=3时，此时位于“初始分支”，升力系数与阻力系数以及横向位移均表现为规则的简弦形式，升力曲线与位移曲线的波峰波谷处在相同位置，表现为典型的同相位角形式。当Ur=4时，圆柱上的升力以及阻力系数以及横向位移曲线出现了“拍”现象，此区域为初始分支向上端分支过渡阶段，说明当圆柱振动从非锁定区进入锁定区时，存在两个相近频率。当Ur=5时，圆柱振动由初始分支跳跃到上端分支，升力系数与阻力系数以及横向运动幅值均有着较大的增加，横向位移与升力系数处于同相位。当Ur=7时，圆柱振动从上端分支跳跃至下端分支了，由图12(d)可知，进入下端分之后，升力系数与阻力系数以及位移均大幅度下降，升力系数与横向振幅曲线由同相位变为反相位，出现“相位切换”现象。当Ur=9时，升力系数山腰处出现“小尖峰”现象，并且升力峰值左边的小尖峰要高于右侧的小尖峰，这是由于锁定区存在三倍频率分量导致的。当Ur=13时，圆柱脱离锁定区，升力系数与阻力系数以及横向振幅均急剧减小，阻力系数几乎减小为0，升力与振动处于反相位角的情形。

(a)Ur=3 (b)Ur=4 (c)Ur=5

(d)Ur=7 (e)Ur=9 (f)Ur=13图12 典型约化速度下升力与阻力随时间变化曲线Fig.12 The lift and drag coefficients versus time in different Ur

图13给出了采用匀减速初始条件时无量纲位移y/D升力系数Cl以及阻力系数Cd随时间变化曲线，当Ur=13时，圆柱在锁定区外，升力系数较小，阻力系数几乎减小为0，并且升力与振动处于反相位角。当Ur=8时，圆柱下端分支中，升力系数依旧在山腰处出现“小尖峰”现象，升力与振动处于反相位角。当Ur=6时，圆柱由下端分支跳跃至上端分支，升力骤然增大同时阻力小幅度增大，升力与振幅由反向变为同相，出现“相位切换”现象。当Ur=5时，圆柱处于上端分支，升力系数与阻力系数增大，升力与位移表现为同相。Ur=4时，处于“初始分支”向“上端分支”过渡区域内，升力以及阻力系数出现了“拍”现象，当Ur=3时，升力系数与阻力系数回归简谐形式，升力与位移表现为同相。

(a)Ur=13 (b)Ur=8 (c)Ur=6

(d)Ur=5 (e)Ur=4 (f)Ur=3图13 典型约化速度下升力与阻力以及横向位移随时间Fig.13 The lift and drag coefficients versus time in different Ur

2.4 从运动轨迹分析涡激振动

图14给出采用匀加速方式时圆柱中心运动轨迹，图中横坐标为流向约化位移，纵坐标为横向约化位移。当Ur=2时，主要发生流向振动，横向振动幅值很小。当Ur=3时，圆柱横向运动幅值增大，流向振动频率变为横向频率两倍，运动轨迹表现为扁平的“8”字型。当Ur=5时，圆柱进入上端分支，此时横向与流向的运动频率较为单一，圆柱表现为瘦长的8字型。当Ur=6.8时，圆柱振动轨迹表现出了超上端分支的特征。圆柱横向振动大幅增大，8字型运动轨迹被拉扯为月牙形。当Ur=8时，横向振动幅值急剧减小至0.7D，圆柱轨迹变为瘦长的8字型。当约化速度达到13.6时，此时圆柱振动已经脱离锁定区，圆柱横向振幅仅剩0.1D左右。

图15给出了三种速度方式下各约化速度时的运动轨迹，通过观察各约化速度下的运动轨迹可以清晰地认识到低质量比双自由度圆柱涡激振动的振动形式衍化过程。对比三种方式可知，采用匀加速的数值模拟方式能够很好的模拟出初始分支，上端分支，超上端分支以及下端分支的运动轨迹现象，其横向振幅呈现递增趋势直至Ur=6.8，当Ur>6.8时，其轨迹突变为瘦长8字型。采用匀减速与匀速轨迹差异主要表现在5≤Ur≤6.8之间，匀减速轨迹在Ur>6时便突变为瘦长8字型。采用匀速方式时，在Ur>5时便突变为瘦长8字型。

(a)Ur=2 (b)Ur=3 (d)Ur=5

(f)Ur=6.8 (g)Ur=8 (h)Ur=14图14 匀加速各约化速度下的轨迹图Fig.14 The VIV trajectory versus time in different Ur

图15 不同初始条件下圆柱涡激振动轨迹图Fig.15 The VIV trajectory versus time in different initial conditions

2.5 从涡脱形式分析迟滞现象

图16给出了匀加速初始条件下圆柱泄涡等值云图(左为Jauvtis等试验DPIV结果，右为CFD结果)。Ur=2时，处于初始分支中，每半个周期便会在圆柱一边卸载一个漩涡，该泄涡模式为“2S”模式。Ur=5，处于上端分支中，每半个周期便会在圆柱一边卸载一对漩涡，同时一个漩涡的涡量强度明显大于另一个漩涡，该泄涡模式为“2P”模式。Ur=6.8，处于超上端分支中，每半个周期便会在圆柱一边卸载三个漩涡，该泄涡模式为“2T”模式，当圆柱横向得到超大振幅时，便会体现出该种泄涡模式。Ur=8时，每半个周期便会在圆柱一边卸载一对漩涡，这两个漩涡强度相当，该泄涡模式为“2P”模式。采用匀加速初始条件成功模拟出“2T”模式，这得益于考虑了涡激振动“迟滞效应”。

图17给出了匀减速初始条件下圆柱泄涡等值云图，Ur=6.8时，泄涡模式表现为“2P”模式，对比图12(c)可知，当Ur=6.8时，采用匀加速与匀减速两种初始条件下，圆柱尾涡模式不同，具体为：采用匀加速时，其尾涡为“2T”模式，与其较小Ur下的模式相同，而采用匀减速时，其尾涡为“2P”模式，与其较大Ur下的模式相同。匀减速Ur=6时，此时流体力超过尾涡维持“2P”模式的最大承受能力，其泄涡模式变为“2T”模式，该泄涡模式伴随着较大的升力，从而导致振幅的突然增大。不同的尾涡模式下升力大小不同，从而导致横向振幅不同以至于产生了“迟滞现象”。图18给出匀速初始条件下圆柱泄涡等值云图，Ur=5时，泄涡模式表现为“2P”模式，Ur=6.8时，泄涡模式依旧为“2P”模式。采用匀速初始条件时，圆柱泄涡一直表现为“2P”模式，该模式难以提供较大振幅，因此匀速初始条件未能捕获上端分支。对比三种初始条件下的尾涡模式可知，圆柱的泄涡模式体现出“惯性”性质，尾涡的倾泄模式总有着维持原倾泄模式的趋势，当约化速度提供的作用力大于圆柱维持原泄涡模式的最大承受能力后，圆柱泄涡模式将出现跳跃，不同初始条件对尾涡的影响进一步体现在振幅的差异上，同Ur下的不同泄涡模式导致产生“迟滞现象”。

(i)DPIV结果 (ii)数值模拟结果 (i)DPIV结果 (ii)数值模拟结果(a)Ur=3涡量等值云图 (b)Ur=5涡量等值云图

(i)DPIV结果 (ii)数值模拟结果 (i)DPIV结果 (ii)数值模拟结果(c)Ur=6.8涡量等值云图 (d)Ur=8涡量等值云图图16 匀加速圆柱泄涡等值云图对比图Fig.16 The vorticity magnitude contours versus Ur in increasing condition

(a)Ur=6.8涡量等值云图(b)Ur=6涡量等值云图(a)Ur=5涡量等值云图(b)Ur=6.8涡量等值云图图17 匀减速数值模拟圆柱泄涡等值云图 图18 匀速数值模拟圆柱泄涡等值云图Fig.17 The vorticity magnitude contours versus Ur in Fig.18 The vorticity magnitude contours versus Ur in decreasing condition using CFD constant condition using CFD

2.6 加速度a对数值结果影响

由上节可知，不同初始条件下的涡激振动体现出很大的差异性，初始条件对圆柱泄涡模式有很大影响，当约化速度提供的作用力大于圆柱维持原泄涡模式的最大承受能力后，圆柱泄涡模式将出现跳跃，同时圆柱振动振幅将出现跳跃，以至于产生“迟滞现象”。增大速度时，圆柱泄涡将倾向于速度较小的泄涡形式。减小速度时，圆柱泄涡将倾向于速度较大的泄涡形式。由于流体具有惯性效应，因此采用匀加速和匀减速进行数值模拟研究是十分有必要的，本节将探讨采用匀加速初始条件下，加速度对结果的影响。分别采用0.012 5，0.025，0.05以及0.1的约化加速度进行数值模拟，约化速度范围为3～10，约化速度间隔为0.5，同时Ur=6.8为临界点，也加以研究。

不同加速度下圆柱横向振幅变化规律如图19所示。采用四种加速度下，均捕获到初始分支上端分支与下端分支。采用约化加速度为0.012 5，0.025的方式捕获到了超上端分支。在3≤Ur≤5.5时，四种加速度方式得到的振幅基本一致，在该分支中，加速度对振动形态影响不大。在5.5

图19 不同加速度下圆柱涡激横向振幅Fig.19 Transverse max amplitude with different acceleration

图20给出了不同加速度下升力均方根值与阻力均值变化曲线。在3≤Ur≤5.5时，四种加速度得到的升力基本一致。在5.5

(a)

(b)图20 四种加速度下升力均方根与阻力均值随约化速度关系Fig.20 The RMS lift coefficient and mean drag coefficient with different accelerations

图21给出了不同加速度下迟滞区间变化曲线。采用约化加速度为0.012 5，0.025的方式捕获到了迟滞区间，其范围为6≤Ur≤6.8。采用约化加速度为0.05，0.1的方式没有捕获到超上端分支，其迟滞区间为5.5≤Ur≤6。由此可见，加速度大小对迟滞区间有较大影响，约化加速度在0.025～0.05存在一个临界值，小于临界值与大于临界值的涡激振动迟滞区域不同。但同在大于临界值或同在小于临界值时，涡激振动得到迟滞区域相同。

图21 不同加速度下涡激振动迟滞区域Fig.21 Hysteresis region of VIV under different acceleration

图22给出Ur从6.5加速到6.8时升力阻力以及位移变化曲线，其采用约化加速度分别为0.025和0.05。可以看到，约化加速度为0.025时，升力系数与阻力系数随着时间逐渐增加，而采用约化加速度为0.05时，升力系数与阻力系数随着时间增加而陡然下降，同时对比横向位移与升力系数可知，此时发生了相位切换。这说明采用较大加速度将会提前导致“相位切换”，从而使得振幅减小。

(a)a*=0.025

(b)a*=0.05图22 Ur=6.8时下升力与阻力随时间变化曲线Fig.22 The lift coefficient and drag coefficient versus time

3 结 论

本文基于PimpleDYMFoam求解器结合SSTk-ω湍流模型，进行了双自由度低质量数值模拟研究，模拟采用了三种初始条件：匀加速，匀减速以及匀速条件。数值模拟捕获到迟滞现象并且三种条件下圆柱振动响应有显著区别，数值模拟对比得出结果如下：

(1)在低质量比圆柱涡激振动试验中，采用加速与减速两种试验条件时，上端分支与下端分支存在迟滞。在“数值水池”中采用匀加速与匀减速方式也可捕获迟滞现象。同时发现，采用匀加速与匀减速条件能捕获到超上端分支，而采用匀速条件则缺失上端分支。在Ur=6.8时，匀加速方式得到最大振幅1.4D，而采用不考虑迟滞效应的匀速条件时，圆柱最大振幅仅为0.85D。由此可知，迟滞效应的存在对振幅的大小以及突变具有显著影响。

(2)迟滞效应的机理可以通过泄涡模式来解释，圆柱泄涡模式体现出“惯性”性质，尾涡的倾泄模式总有着维持原倾泄模型的趋势，即增大速度时，圆柱泄涡将倾向于速度较小的泄涡形式。减小速度时，圆柱泄涡将倾向于速度较大的泄涡形式。当约化速度提供的作用力大于圆柱维持原泄涡模式的最大承受能力后，圆柱泄涡模式将出现跳跃。在迟滞区间内，同Ur下采用不同初始条件时的泄涡模式是不同的，同Ur下不同的泄涡模式导致了“迟滞现象”。

(3)采用匀加速匀减速以及匀速条件均捕获到“锁定”现象，但是锁定区间不同。上端分支跳跃至下端分支的约化速度大小与初始条件有关。同时发现，当Ur≤4.2时，流向振动有两个主要频率，分别为单倍与双倍横向频率。在锁定区域内，横向振动除主频外，还存在一个3倍的低频分量。

(4)作用在圆柱体上的升力与阻力也出现了迟滞现象，其迟滞区域为6≤Ur≤6.8，并在约化速度为6和6.8时，通过升力与位移关系捕获了相位切换现象。升力与横向位移刚开始为同向，随着流速增加，泄涡模式发生突变导致圆柱受到的升力发生突变，从而导致升力与横向位移变为反向，认为是泄涡模式的突变导致了相位切换。

(5)匀加速条件出现了“月牙形”运动轨迹，而匀减速与匀速只捕获了经典“8”型运动。同时发现，在迟滞区间内，相同约化速度下的运动轨迹形式是不同的，圆柱轨迹也与初始条件有一定关系。

(6)在迟滞区域内，在约化加速度0.025与0.05间存在临界值，大于该临界值，升力与位移会提前发生相位切换，小于该临界值，相位切换将会延后。数值模拟中，控制约化加速度小于等于0.025是其捕获超上端分支的一个必要条件。

[ 1 ] 白治宁, 肖龙飞, 程正顺,等. 深吃水半潜式平台涡激运动响应模型实验研究[J]. 船舶力学, 2014(4):377-384.

BAI Zhining, XIAO Longfei, CHENG Zhengshun, et al. Experimental study on vortex induced motion response of a deep draft semi-submersible platform[J]. Journal of Ship Mechanics， 2014(4): 377-384.

[ 2 ] 曹淑刚, 黄维平, 顾恩凯. 考虑流固耦合的弹性圆柱体涡激振动研究[J]. 振动与冲击, 2015, 34(1): 58-62.

CAO Shugang, HUANG Weiping, GU Enkai. Vortex-induced vibration of an elastic cylinder considering fluid-structure interaction[J]. Journal of Vibration & Shock, 2015, 34(1): 58-62.

[ 3 ] 黄维平, 刘娟, 唐世振. 考虑流固耦合的大柔性圆柱体涡激振动非线性时域模型[J]. 振动与冲击, 2012, 31(9): 140-143.

HUANG Weiping, LIU Juan, TANG Shizhen. Nonlinear model of vortex induced vibration of flexible cylinder in consideration of fluid-structure interaction[J]. Journal of Vibration & Shock, 2012, 31(9): 140-143.

[ 4 ] 范杰利, 黄维平. 细长立管两向自由度涡激振动数值研究[J]. 振动与冲击, 2012, 31(24): 65-68.

FAN Jieli, HUANG Weiping. Numerical simulation of 2-DOF vortex-induced vibration of a long riser[J]. Journal of Vibration & Shock, 2012, 31(24): 65-68.

[ 5 ] 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.

[ 6 ] 谷家扬, 杨建民, 肖龙飞. 两种典型立柱截面涡激运动的分析研究[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.

[ 7 ] GUILMINEAU E, QUEUTEY P. Numerical simulation of vortex-induced vibration of a circular cylinder with low mass-damping in a turbulent flow[J]. Journal of Fluids and Structures, 2004, 19(4): 449-466.

[ 8 ] 陈威霖, 及春宁, 徐万海. 并列双圆柱流致振动的不对称振动和对称性迟滞研究[J]. 力学学报, 2015, 47(5): 731-739.

CHEN Weilin, JI Chunning, XU Wanghai. Numerical investigation on the asymmetric vibration and symmetry hysteresis of flow-induced vibration of two side-by-side cylinders[J]. Chinese Journal of Theoretical & Applied Mechanics, 2015, 47(5): 731-739.

[ 9 ] KIM S, WILSON P A, CHEN Z M. Numerical simulation of force and wake mode of an oscillating cylinder[J]. Journal of Fluids & Structures, 2014, 44(7): 216-225.

[10] KHALAK A, WILLIAMSON C H K. Motions, forces and mode transitions in vortex-induced vibrations at low mass-damping[J]. Journal of Fluids & Structures, 1999, 13(7/8): 813-851.

[11] MINGUEZ M，LUPPI A，BERGER A. Slender buoy FSHR vortex induced rotations[C]// Proceedings of the ASME 31th International Conference on Ocean, Offshore and Arctic Engineering. Rio de Janeiro: OMAE, 2012.

[12] CHEN W, JI C, XU W, et al. Response and wake patterns of two side-by-side elastically supported circular cylinders in uniform laminar cross-flow[J]. Journal of Fluids & Structures, 2015, 55: 218-236.