表面粗糙度对圆柱体涡激振动响应特性影响数值研究

高 云, 杨家栋, 邹 丽, 宗 智

(1. 西南石油大学 油气藏地质及开发工程国家重点实验室,成都 610500;2.东京大学 机械工程学院,东京 113-8656;3.大连理工大学 船舶工程学院,大连 116024)

圆柱体在一定速度的来流下，会在其两侧形成交替脱落的漩涡，周期性的漩涡脱落会在圆柱体上产生周期性的横向(Cross-Flow, CF)升力以及流向(In-Line, IL)拖曳力。若圆柱体为弹性支撑，周期性的升力以及拖曳力会引起圆柱体在横向以及流向发生振动，称之为涡激振动(Vortex-Induced Vibration, VIV)[1]。圆柱体涡激振动又会发过来作用于其尾部流场，从而影响到作用在圆柱体上的水动力载荷。因此，圆柱体涡激振动属于典型的非线性流固耦合问题[2]，受周围流场特性的影响显著，而圆柱体的表面粗糙度是影响圆柱体周围流场的一个重要参数。

关于表面粗糙度对圆柱体周围流场以及圆柱体VIV响应的影响，已有较多的学者进行了相关研究。根据加以研究的圆柱体其自身边界条件，可分为静止圆柱体(完全固定)研究[3-8]以及振荡圆柱体(弹性支撑)研究[9-14]。对于静止圆柱体，主要是针对圆柱体周围流场特性进行研究，包括漩涡泄放频率、漩涡泄放模式以及由漩涡泄放引起的升力以及拖曳力等参数。对于振荡圆柱体，除了针对圆柱体周围流场特性进行研究外，还需要对结构振动响应特性进行研究，包括圆柱体的VIV响应幅值以及响应频率等参数。

如表1所示，早期学者主要针对不同粗糙度下空气中圆柱体的尾部流场特性进行了研究。通过研究发现: 当流体从层流变成湍流时，作用在圆柱体上的拖曳力会极速下降，这种现象称为拖曳力危机(drag crisis)。且随着表面粗糙度的上升，出现拖曳力危机的临界Re数逐渐变小，当表面粗糙度达到一定值(3.0×10-3)时，临界区域、超临界区域以及超高临界区域会合并到一个窄带区域。近些年，由于海洋工程的快速发展，关于表面粗糙度对水中圆柱体VIV响应特性的研究则得到了关注。通过研究表明: 随着表面粗糙度的上升，圆柱体的最大位移响应以及最大拖曳力均值逐渐降低，且趋于恒定值。与光滑立管相比，粗糙立管的斯脱哈尔数呈上升趋势。

表1 不同表面粗糙度下的圆柱体研究

从目前的研究现状来看,针对水中粗糙圆柱体的尾部流场特性以及结构响应特性的研究还是很缺乏，尚存在很多问题需要进行更深入的研究。比如当圆柱体尾部流场漩涡泄放频率接近结构固有频率时，便会发生锁定现象(lock-in phenomenon)，锁定区间内，圆柱体会发生大幅的、危险的涡激振动。随着折合速度的增加，具有不同粗糙度圆柱体的锁定区间以及锁定频率会发生怎样的变化？基于此，本文对表面粗糙度对圆柱体涡激振动响应特性影响进行了数值研究。对不同粗糙度下圆柱体CF以及IL方向的涡激振动响应幅值、漩涡泄放频率、结构振动频率、结构振动轨迹、锁定区间进行了系统地分析和讨论。

1 问题描述

本文研究流场如图1所示，在均匀流场中放置一个直径为D的圆柱，圆柱直径D取为0.038 1 m。流场长度取为32D，流场宽度取为22D。因此阻塞率D/22D=0.045。当阻塞率小于0.05时，流场宽度对圆柱体响应的影响可以忽略[15-16]。入流边界为一定速度的不可压缩均匀流，流体为标准状态下的液态水，从入口边界以恒定速度U∞流向圆柱。

如图1所示，入口边界条件为u=U∞以及v=0，入口速度U∞与雷诺数(Re数)相互对应；出口边界条件为∂u/∂x=0和∂v/∂x=0；上下边界条件满足∂u/∂y=0以及v=0。圆柱表面为无滑移壁面，在无滑移壁面上满足∂u/∂y=0和v=0。

图1 流场区域以及边界条件

这个系统可以进行流向和横向的二维振动。圆柱质量比(m*)取为2.6(m*=m/mw，其中m和mw分别是单位长度的圆柱体质量以及单位长度圆柱体所排开的流体质量，m=ρcπD2/4，mw=ρπD2/4，ρc和ρ分别是圆柱体和流体的密度)。本文Re数保持不变，Re=5 000，也就意味着流场的入流速度U∞保持不变(U∞=Re·ν/D，其中ν为流体的运动黏性系数)。通过改变固有频率的大小来改变折合速度，在过去的很多数值方法中，都采用了通过改变固有频率的方法来改变折合速度[17-19]。固有频率fn与折合速度Vr之间的关系为Vr=U∞/(fn×D)。计算中取了4种不同粗糙度的圆柱体，粗糙度系数Ks/D分别取为0、5×10-3,1×10-2以及2×10-2，对于每种圆柱体，计算了分布在1～14的21种不同折合速度。因此，一共计算了84种工况。

2 控制方程以及数值方法

2.1 控制方程

圆柱体涡激振动响应系统可以看作一个二自由度质量-弹簧-阻尼模型，二自由度振动圆柱体的动力学方程为[20-21]

(1)

(2)

式中:c,k以及m分别为系统的阻尼系数、刚度以及单位长度圆柱体的质量，这3个参数的取值来源于Jauvtis等[22]的实验数据。圆柱体边界条件取为

(3)

圆柱体初始条件取为

(4)

2.2 数值方法

由于本文将圆柱体视为在xoy平面内刚性运动的质点，因此圆柱动力学方程可以分别依照单自由度质点非线性振动独立求解。采用四阶龙格库塔法对动力学方程进行离散。以横向振动为例，方程可以表达为

(5)

使用经典四阶龙格库塔法，式(5)可以离散为

(6)

式中：K1,K2,K3,K4,L1,L2,L3以及L4分别为

(7)

如图2所示，本文使用等效砂粒粗糙度模型来模拟圆柱体的表面粗糙度。假设粗糙度的存在会导致一个阻塞效应，大概占粗糙高度的50%。粗糙度的存在使圆柱体的直径增加了一个粗糙高度。即D′=D+Ks, 其中D′为校正后的圆柱体直径。

图2 等效砂粒粗糙度模拟

3 分析与讨论

3.1 网格独立性以及数值模型验证

本文的网格由网格划分软件ICEM CFD生成，如图3所示。

计算域被分成了两个块：随动区域和变形区域。随动区域是用四边形网格划分的，变形区域是用三角形网格划分的。随动区域的网格随着圆柱的振动同步变动；变形区域的网格可以根据随动区域的移动进行调整和适应。为了验证网格的无关性，对光滑圆柱体进行了验证，验证中选取的参数为：Vr=8.0,m*=2.6,CA=1,ζ=0.003 6。对4种不同的网格进行了计算，计算得到的Ay/D和fs/fn如表2所示。Ay/D为圆柱体在CF方向的无量纲最大位移，fs/fn为结构振动频率与固有频率的比值。

(a) 整个区域网格

(b) 随动区域网格

表2网格独立性验证(Vr=8.0,m*=2.6,CA=1,ζ=0.0036)

Tab.2Meshindependencystudy(Vr=8.0,m*=2.6,CA=1,ζ=0.0036)

网格单元数Ay/Dfs/fnM1155700.6021.24M2197400.613(1.83%)1.274(2.7%)M3237180.620(1.30%)1.281(0.55%)M4263230.621(0.16%)1.282(0.08%)

从表2可以看出：网格从M1加密至M2时，Ay/D和fs/fn的变化幅值分别为1.83%和2.7%；从M2加密至M3时，Ay/D和fs/fn的变化幅值缩小为1.30%和0.55%；当网格从M3加密至M4，Ay/D和fs/fn的变化幅值进一步缩小至0.16%和0.08%，即M4计算得到的Ay/D和fs/fn与M3非常接近。在Intel Core i5-4590平台上，使用M4完成计算用时为12 h，而M3完成计算用时为9 h，M4占用计算资源明显高于M3，在同样满足计算精度的条件下，使用M3计算更加节约计算成本，因此本文中选择M3作为计算网格。

为了验证本文数值模型的可靠性，这里对本文的数值模型进行了验证。验证时选取的参数与Jauvtis等(2004)的实验参数相同(m*=2.6,CA=1,ζ=0.003 6)。图4和图5依次给出了本文数值模型计算得到的最大无量纲位移Ay/D以及频率比fs/fn与Jauvtis等的实验结果以及Zhao等[23]的数值结果(Zhao等的参数为：m*=2.6,CA=1,ζ=0.003 6)的对比。由图4的位移响应可以看出位移响应存在明显的4个区间，分别为初始分支、上分支、低分支以及解锁区间。值得注意的是：本文的数值结果以及Zhao等的数值结果均没有得到超上端分支，原因是由采用改变固有频率的方法来得到不同的折合速度所导致。由图5的频率响应可以很明显地看出锁定区间和非锁定区间，在非锁定区间里，振动频率与计算得到的斯脱哈尔漩涡泄放频率吻合；而在锁定区间里，振动频率基本会出现在固有频率附近，偏移斯脱哈尔漩涡泄放频率。总体上来说，本文的数值结果与Jauvtis等的实验结果以及Zhao等的数值结果吻合良好。

图4 本文数值模型光滑圆柱体涡激振动响应幅值与Jauvtis等以及Zhao等的结果对比

Fig.4 Comparison of response amplitudes of smooth cylinder with Jauvtis’s and Zhao’s

图5 本文数值模型光滑圆柱体涡激振动响应频率与Jauvtis等以及Zhao等的结果对比

Fig.5 Comparison of response frequencies of smooth cylinder with Jauvtis’s and Zhao’s

3.2 涡激振动响应特性分析

图6到图9分别给出了光滑圆柱体以及带有三种不同粗糙度的圆柱体的CF以及IL方向的无量纲最大位移、尾部流场漩涡泄放频率、结构振动频率以及锁定区间。尾部流场漩涡泄放频率与升力相关，对于稳态流中的圆柱体涡激振动，漩涡泄放频率与升力频率是相等的。漩涡泄放频率以及结构固有频率可分别对升力时间历程以及结构位移响应求快速傅里叶(Fast Fourier Transform, FFT)变换后取主导频率得到。这里简要地介绍一下如何判断锁定区间。对于刚性圆柱体，传统上认为当结构振动频率(fs)、尾部流场漩涡泄放频率(fv)与静止流体中的结构固有频率(fn)一致时圆柱体处于锁定状态。这种判定方法非常适合高质量比m*=O(100)(如空气)下，此时频率比f*=fs/fn非常接近1[24-25]；但当结构在低质量比m*=O(10)(如水中)流体中发生锁定时，f*要很明显偏移1[26-28]。

Sarpkaya[29]对这种偏移现象给出了很好的解释：通常我们对振动频率进行无量纲化所采用的结构固有频率是处在静止流体中的固有频率。但当结构发生振动时，结构的真实固有频率会逐渐偏移静止流体中的结构固有频率，这会导致与结构真实固有频率吻合的结构振动频率逐渐偏移静止流体中的结构固有频率。这个偏移量会随着质量比的减小而呈现上升趋势。当结构处在高质量比流体中，这个偏移量可以忽略，这就会出现结构的振动频率非常接近结构在静水中的固有频率；但当结构处在低质量比流体中，这个偏移量较为明显，最大偏移值可以达到1.5。综合以上因素，这里我们给出了判断刚性圆柱体锁定区间的方法，锁定区间的涡激振动响应特性需满足以下几个条件：① 结构的振动频率与流场的漩涡泄放频率吻合；② 结构振动频率以及流场漩涡泄放频率与结构在静水中的固有频率比较接近；③ 结构的涡激振动响应轨迹呈现规则的8字形状。

图6 光滑圆柱体无量纲涡激振动响应幅值、结构振动频率、漩涡泄放频率以及锁定区间

Fig.6 Nondimensional displacements, structural vibration frequencies, vortex shedding frequencies and lock-in regions versus reduced velocity for the smooth cylinder

图7 粗糙圆柱体(Ks/D=5×10-3)无量纲涡激振动响应幅值、结构振动频率、漩涡泄放频率以及锁定区间

Fig.7 Nondimensional displacements, structural vibration frequencies, vortex shedding frequencies and lock-in regions for the rough cylinder withKs/D=5×10-3

图8 粗糙圆柱体(Ks/D=1×10-2)无量纲涡激振动响应幅值、结构振动频率、漩涡泄放频率以及锁定区间

Fig.8 Nondimensional displacements, structural vibration frequencies, vortex shedding frequencies and lock-in regions for the rough cylinder withKs/D=1×10-2

图9 粗糙圆柱体(Ks/D=2×10-2)无量纲涡激振动响应幅值、结构振动频率、漩涡泄放频率以及锁定区间

Fig.9 Nondimensional displacements, structural vibration frequencies, vortex shedding frequencies and lock-in regions for the rough cylinder withKs/D=2×10-2

基于以上判定方法，由图6～图9中的响应频率以及振动轨迹，我们可以将整个区间分成锁定区间以及非锁定区间。锁定区间可以按照幅值变化趋势细分为两个不连续的区间，记为区间II、区间III。区间II里，响应幅值随着折合速度的上升呈增加趋势，当响应幅值呈现下降趋势时则进入区间III，区间II和区间III对应涡激振动响应的“上端分支”以及“低分支”。进入锁定区间前的非锁定区间记为区间I，对应涡激振动响应的“初始分支”。锁定区间结束后的非锁定区间记为区间IV，对应涡激振动响应的“解锁区间”。

表3依次给出了光滑圆柱体以及3个不同粗糙度圆柱体对应的4个区间的折合速度分布范围。由表1可以看出：初始分支分布区间随着粗糙度的上升基本不发生变化，只有当粗糙度从1×10-2增加到2×10-2时，区间结束点对应的折合速度减小了0.5；上端分支分布区间随着粗糙度的上升呈现明显的减小趋势；低分支分布区间同样随着粗糙度的上升呈现明显的降低趋势；随着粗糙的上升，由上端分支和低分支联合组成的整个锁定区间锁定开始点呈现缓慢提前的趋势，而锁定结束点呈现明显的提前趋势，导致整个锁定区间的宽度随着粗糙的上升而逐渐减小。

表3 折合速度划分区间

由图6可以看出：当圆柱体为光滑圆柱体时，当Vr较小处在区间I时，此时结构振动频率与漩涡泄放频率吻合，且二者与斯脱哈尔漩涡泄放频率fst=St×V/D,St=0.18非常接近，但远离结构固有频率，此时圆柱体没有发生锁定，在大多数情况下结构振动轨迹呈现非常混乱的非8字形状。值得注意的是当Vr=3时，振动轨迹呈现规则的8字形状，这是由IL方向产生的Pure lock in所导致[30]。随着Vr的上升，由区间I进入区间II时，振动轨迹会出现由非8字混乱形状变换到规则8字形状的现象。当折合速度处在区间II时，此时结构振动频率与漩涡泄放频率吻合，且二者逐渐偏移斯脱哈尔漩涡泄放频率fst，但接近结构的固有频率，圆柱体发生锁定，锁定区间内振动幅值呈现明显的上升趋势且振动轨迹呈现明显的规则8字形状。在上分支区间，频率比f*在0.9附近；在低分支区间频率比f*在1.3附近。当Vr=6时，CF以及IL方向的响应幅值依次达到最大值1.08D以及0.2D。随着Vr的进一步上升，由区间II进入区间III，此时CF以及IL方向的振动响应幅值呈现明显的突降趋势，此时结构的振动频率会出现一个明显的跳跃增加现象，这是由漩涡泄放模式发生转变所导致。当Vr处在区间III时，随着Vr的上升，CF方向的振动响应幅值基本处于稳定状态，稳定在0.6D附近。当Vr进入区间IV时，此时结构振动频率脱离漩涡泄放频率，进入解锁区间，振动轨迹则再次出现混乱形状，且振动响应幅值呈现下降趋势。

对比图7、图8以及图9与图6进行对比可以看出：随着粗糙度的上升，CF方向的最大无量纲位移呈现递减趋势，由光滑圆柱体的1.08D逐渐减小到具有最大粗糙度的圆柱体的0.78D。对于光滑圆柱体以及粗糙度为5.0×10-3以及1.0×10-2的粗糙圆柱体，CF方向的最大位移响应均出现在区间II(上端分支)，但对于粗糙度为2.0×10-2的粗糙圆柱体，CF方向的最大位移响应则出现在区间III(低分支)。当圆柱体为光滑圆柱体或粗糙度为5.0×10-3的小粗糙圆柱体时，当Vr从区间II进入区间III时，振动响应幅值会出现明显的跳跃下降现象，而振动响应频率则会出现明显的跳跃增加现象，由0.9fn直接跳跃到1.25fn；但当圆柱体为粗糙度为1.0×10-2以及1.0×10-2的大粗糙圆柱体时，当Vr从区间II进入区间III时，振动响应幅值呈缓慢下降趋势，此时振动响应频率呈现缓慢上升趋势，最后基本稳定在1.20fn。

4 结 论

本文针对粗糙度对圆柱体涡激振动响应特性的影响进行了数值分析。圆柱的质量比、阻尼比以及Re数分别取为2.6,0.003 6以及5 000。对光滑圆柱体以及3种具有不同粗糙度的粗糙圆柱体的涡激振动响应特性进行了研究，研究参数包括：涡激振动位移响应、结构振动频率、漩涡泄放频率、锁定区间等。通过以上研究，可得到如下结论：

(1) 根据涡激振动响应幅值、响应频率以及振动轨迹，整个折合速度区间可以大致分为4个区间(区间I、区间II、区间III以及区间IV)，这4个区间依次对应涡激振动响应的初始分支、上分支、低分支以及解锁区域。

(2) 区间II和区间III同属锁定区间，锁定区间内，CF方向的涡激振动响应要明显大于IL方向，且该区间响应轨迹呈现明显的竖8字形状。区间I和区间IV属于非锁定区间。当Vr处于区间I内，涡激振动响应会出现Pure-lock-in现象，此时圆柱体IL方向的涡激振动响应与CF方向的涡激振动响应非常接近。

(3) 随着粗糙度的上升，圆柱体涡激振动响应最大值呈下降趋势；随着粗糙度的上升，锁定区域开始点对应的Vr呈缓慢提前的趋势，而锁定区域结束点对应的Vr呈快速提前的趋势，因此整个锁定区域宽度会随着粗糙度的上升而逐渐变窄。

参 考 文 献

[1] 高云, 付世晓, 熊友明,等. 剪切来流下柔性圆柱体涡激振动响应试验研究[J]. 振动与冲击, 2016, 35(20): 142-148.

GAO Yun, FU Shixiao, XIONG Youming, et al. Experimental study on vortex induced vibration responses of a flexible cylinder in sheared current[J]. Journal of Vibration and Shock, 2016, 35(20): 142-148.

[2] 杨林. 非线性流固耦合问题的数值模拟方法研究[D]. 青岛: 中国海洋大学, 2011.

[3] ACHENBACH E. Influence of surface roughness on the cross-flow around a circular cylinder[J]. Journal of Fluid Mechanics, 1971,46(2): 321-335.

[4] ACHENBACH E, HEINECKE A. On vortex shedding from smooth and rough cylinders in the range of Reyonlds numbers 6×103to 5×106[J]. Journal of Fluid Mechanics,1981,109:239-251.

[5] NAKAMURA Y, TOMONARI Y. The effect of surface roughness on the flow past circular cylinders at high Reynolds numbers[J]. Journal of Fluid Mechanics, 1982,123: 363-378.

[6] RIBEIRO L J D. Effects of surface roughness on the two-dimensional flow past circular cylinders I: mean forces and pressures[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1991,37(3): 299-309.

[7] RIBEIRO L J D. Effects of surface roughness on the two-dimensional flow past circular cylinders II: fluctuating forces and pressures[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1991,37(3): 311-326.

[8] BEARMAN P W, HARVEY J K. Control of circular cylinder flow by the use of dimples[J]. AIAA Journal, 1993,31(10): 1753-1756.

[9] OKAJIMA A, NAGAMORI T, MATSUNAGA F, et al. Some experiments on flow-induced vibration of a circular cylinder with surface roughness[J]. Journal of Fluids and Structures, 1999,13(7): 853-864.

[10] ALLEN D W, HENNING D L. Surface roughness effects on vortex-induced vibration of cylindrical structures at critical and supercritical Reynolds numbers[C]∥Proceedings of the Offshore Technology Conference. Houston: OTC, 2001.

[11] BERNITSAS M M, RAGHAVAN K, DUCHENE G. Induced separation and vorticity using roughness in VIV of circular cylinders at 8×103

[12] BERNITSAS M M, RAGHAVAN K. Reduction/suppression of VIV of circular cylinders through roughness distribution at 8×103

[13] KIU K Y, STAPPENBELT B, THIAGARAJAN K P. Effects of uniform surface roughness on vortex-induced vibration of towed vertical cylinders[J]. Journal of Sound and Vibration, 2011,330(20): 4753-4763.

[14] GAO Y, FU S, WANG J, et al. Experimental study of the effects of surface roughness on the vortex-induced vibration response of a flexible cylinder[J]. Ocean Engineering, 2015, 103: 40-54.

[15] ZHAO M, TONG F, CHENG L. Numerical simulation of two-degree-of-freedom vortex-induced vibration of a circular cylinder between two lateral plane walls in steady currents[J]. Journal of Fluids Engineering, 2012, 134(10): 377-390.

[16] NAVROSE, MITTAL S. Free vibrations of a cylinder: 3D computations atRe=1 000[J]. Journal of Fluids and Structures, 2013, 41: 109-118.

[17] LEONTINI J, THOMPSON M, HOURIGAN K. The beginning of branching behavior of vortex-induced vibration during two-dimensional flow[J]. Journal of Fluids and Structures, 2006, 22(6/7): 857-864.

[18] BOURGUET R, JACONO D. Flow-induced vibrations of a rotating cylinder[J]. Journal of Fluid Mechanics, 2014, 740: 342-380.

[19] ZHAO M, CUI Z, KWOK K, et al. Wake-induced vibration of a small cylinder in the wake of a large cylinder[J]. Ocean Engineering, 2016, 113: 75-89.

[20] PRASANTH T, MITTAL S. Vortex-induced vibration of two circular cylinders at low Reynolds[J]. Journal of Fluids and Structures, 2009, 25(4): 731-741.

[21] BAO Y, HUANG C, ZHOU D, et al Two-degree-of-freedom flow induced vibrations of isolated and tandem cylinders with varying natural frequencies[J]. Journal of Fluids and Structures, 2012, 35: 50-75.

[22] 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, 59: 23-62.

[23] ZHAO M, KAJA K, XIANG Y, et al. Vortex-induced vibration (VIV) of a circular cylinder in combined steady and oscillatory flow[J]. Ocean Engineering, 2013, 73: 83-95.

[24] FENG C. The measurement of Vortex-induced effects in flow past a stationary and oscillating circular and D-section cylinders[D]. Vancouver: University of British Columbia, 1968.

[25] SUMER B M, FREDSØE J. Hydrodynamics around cylindrical structures[J]. Singapore: World Scientific,1997:26.

[26] MOE G, WU Z J. The lift force on a cylinder vibrating in a current[J]. ASME Journal of Offshore Mechanics and Arctic Engineering,1990,112(4):297-303.

[27] KHALAK A, WILLIAMSON C H K. Fluid forces and dynamics of a hydroelastic structure with very low mass and damping[J]. Journal of Fluids and Structures,1997,11(8): 973-982.

[28] KHALAK A, WILLIAMSON C H K. Motions, forces and motion transitions in vortex-induced vibration at low mass-damping[J]. Journal of Fluids and Structures,1999, 13: 813-851.

[29] SARPKAYA T. Fluid forces on oscillating cylinders[J]. Nasa Sti/recon Technical Report A,1979,104(3): 275-290.

[30] ARONSEN K H. An experimental investigation of in-line and combined in-line and cross-flow vortex induced vibrations[D]. Trondheim: NTNU, 2007.