质量比对近壁面两向自由度圆柱涡激振动的影响

刘旭菲， 陈威霖， 及春宁

(1. 浙江水利水电学院 水利与环境工程学院，杭州 310018; 2. 天津大学 水利工程仿真与安全国家重点实验室，天津 300072)

过去几十年对圆柱绕流和涡激振动(vortex-induced vibrations，VIV)的研究多限定在无限区域内的自由来流条件[1-7]。然而，实际工程多涉及有壁面存在的情况，如海床上的输油管道、电缆以及近壁的热交换管等。壁面的存在使得近壁面一侧剪切层的发展受到抑制，会改变圆柱的受力和振动响应。因此，对近壁面圆柱涡激振动的研究具有重要的工程价值和科学意义。

对近壁面圆柱绕流而言，无量纲壁面边界层厚度(δ/D，D为圆柱直径)、间隙比(G/D，G为圆柱下表面与壁面之间的距离)以及雷诺数(Re)等对圆柱周围流场以及脱落旋涡与壁面边界层的相互作用有显著影响，并因此影响圆柱受到的流体力[8-9]。根据G/D的不同，Price等[10]将近壁面圆柱绕流分成了4个不同的区域：当G/D=0.125时，旋涡脱落被抑制，仅上侧剪切层在较远下游出现小幅摆动；当G/D=0.250～0.375时，圆柱下侧剪切层开始与上侧剪切层发生相互作用，使得圆柱上侧发生旋涡脱落；当0.501.0时，壁面的影响几乎可以忽略，近壁圆柱旋涡脱落与孤立圆柱的情况相同。类似地，Lei等[11]发现当G/D=0.2～0.3时，圆柱旋涡脱落被抑制，而Buresti等[12]试验研究发现当G/D<0.4时，圆柱旋涡脱落被完全抑制。此外，随δ/D的增加，旋涡脱落被抑制的临界G/D增加。Zdravkovich[13]发现当圆柱位于壁面边界层之外时，G/D对圆柱阻力的影响较弱，而当圆柱位于壁面边界层之内时，阻力随G/δ的增加而减小。

当圆柱可以自由振动时，圆柱与壁面剪切层发生复杂的相互作用，近壁面的影响变得更加显著。Wang等[14]研究发现，壁面效应会显著改变圆柱的振幅和频率，使得圆柱与壁面边界层之间发生显著的非线性相互作用。Zhao等[15]数值模拟研究发现，即使间隙比很小(G/D=0.002)，圆柱的涡激振动仍可以被激发。根据圆柱与壁面作用的不同，他们将脱涡模式分成了单涡模式、碰撞反弹前脱涡模式以及碰撞反弹后脱涡模式。Yang等[16]试验研究发现，当G/D>0.66时，圆柱的振动响应明显不同于G/D<0.3的情况。Li等[17]发现壁面抑制了圆柱下侧的剪切层，使得流向和横向振动具有相同的振动频率。Tham等[18]数值模拟研究发现，当Re=100和m*=10时，随G/D的减小，锁定区域变窄。当G/D<1.0时，交替脱落旋涡模式被明显抑制；当G/D≤0.6时，在初始和下端分支之间出现了第三分支。Li等[19]通过二维和三维的数值模拟研究了Re=100,300的近壁面圆柱涡激振动，其中间隙比固定为G/D=0.9。研究发现，随δ/D的增加，圆柱的振动开始于更大的折合流速上，共出现了4种不同的尾流模式：W2S(A) ——旋涡从圆柱两侧脱落，且下侧旋涡可以向下游移动；W2S(B)——与W2S(A)相似，但下侧旋涡更弱，且很快消失；1S ——仅圆柱上侧脱落旋涡；NS——无旋涡从圆柱脱落。Barbosa等[20]研究发现，当G/D>2.0时，壁面对锁定区域内的振动影响可以忽略；当G/D=0.75～2.00时，壁面边界层使得圆柱振幅减小，但流场基本保持对称；当G/D<0.75时，振动的圆柱在某些情况下与壁面发生碰撞。对近壁面圆柱涡激振动轨迹的研究发现，近壁面圆柱涡激振动的轨迹呈椭圆形，与孤立圆柱的8字型轨迹形成了鲜明的对比[21-23]。陈蓥等[24]通过试验研究了均匀流下强迫振动的近壁面圆柱的水动力特性，其中雷诺数为Re=2×105。他们发现圆柱受到的平均阻力随G/D的减小而下降；近壁面的存在对圆柱的能量传递有重要的影响。杨兵等[25]通过水槽试验研究了固壁条件下单圆柱的流向振动，发现随着来流流速的增加圆柱的流向振动经历发生、发展及最后消失的过程。刘俊等[26]通过试验研究了近壁面涡激振动触发和停振的临界流速，发现触发的临界流速随G/D的减小而增加。

以上的研究表明，近壁面圆柱涡激振动受雷诺数和间隙比等因素的影响显著。然而，孤立圆柱涡激振动的研究表明，质量比也是影响涡激振动的重要因素[27-31]。但已有的研究较少关注质量比对近壁面圆柱涡激振动的影响。为此，本文选取3个质量比(m*=2，10和20)，对质量比对近壁面圆柱涡激振动的影响进行深入研究，以期对近壁面圆柱涡激振动机理进行更深入的探讨。

1 数值方法及其验证

1.1 控制方程

流固耦合问题的数值模拟采用浸入边界法[32]，控制方程如下

(1)

∇·u=0

(2)

式中：u为速度；t为时间;p为压强;ρ为流体密度;υ为运动黏滞系数；∇为梯度算子；f为附加体积力矢量，代表流固耦合边界条件。

对以上控制方程采用二阶的Admas-Bashforth时间格式进行离散，可得守恒形式如下

(3)

∇·un+1=0

(4)

式中，h=∇·(-uu+ν(∇u+∇uT))由对流项与扩散项组成，上标T为矩阵转置，附加体积力表示为

(5)

式中：I(φ,Xi)和D(Φ,x)为引入的插值函数，φ为定义在流体网格节点x上的物理量，如流体的速度u、压强p等,Φ为定义在物面边界点Xi上的物理量，如固体速度Vn+1等;上标n+1，n+1/2，n，n-1为时间步;Δt为时间步长。

针对传统浸入边界法施加边界条件精度不高的情况，Ji等提出了基于嵌入式迭代的浸入边界法，将浸入边界法嵌入到压强泊松方程的迭代求解中，利用压强的中间解比初始值更接近真实值的特点，迭代修正附加体积力，在不显著增加额外计算耗时的前提下，提高整个算法的求解精度。有关浸入边界法的细节，请参考Ji等的研究，此处不再赘述。

对做两自由度运动的刚性圆柱，其运动方程可以用下述方程来描述

(6)

1.2 模型设置

在近壁面圆柱涡激振动的数值模拟中，计算域的流向和横向宽度分别为Lx和H。如图1所示，Lx=Lu+Ld，其中Lu=9D和Ld=80D；H=H1+H2，其中H1=1.1D和H2=97D。此外，为保证模拟结果的准确性，在圆柱周围设置一个大小为8D×8D的加密区域，加密区域内流向和横向网格大小均为1/64D。加密区域外的网格大小以一定的比例增加[33-36]。计算域边界条件设置如下。入口为Dirichlet边界条件，出口为Neumann边界条件，上边界为自由可滑移边界，下边界为不可滑移边界。无量纲时间步长为ΔtU∞/D=0.005。需要说明的是本文所采用网格精度和时间步长均与Chen等研究中的一致，参数收敛性得到了充分验证，读者可自行查阅，此处不再赘述。

图1 近壁面圆柱涡激振动的计算域和边界条件设置Fig.1 Computational domain and boundary conditions for VIV of a near-wall cylinder

1.3 数值模拟验证

以近壁面圆柱涡激振动为例，验证本文数值方法的准确性。为保证结果的可比性，这里采用与验证算例相等的计算域：圆柱中心距离入口为29D，距离出口为45D，圆柱中心距离上边界为10D。间隙比G/D=0.6，折合流速Ur=6.0，质量比m*=10，阻尼比为零。圆柱可两向自由度振动。如表1所示，本文的模拟结果与Tham等吻合很好，最大的误差仅为2.4%，说明了本文数值方法的准确性。表1中，横向无量纲最大振幅定义为Ymax/D=(ymax-ymin)/2D, 其中:ymax为最大位移和ymin为最小位移；流向无量纲均方根振幅Xrms/D，其中Xrms为流向位移的均方根值；Clrms和Cdrms分别为升力系数和阻力系数均方根值。

表1 近壁面单圆柱涡激振动的数值模拟结果与验证算例[18]的对比Tab.1 Comparison of VIV of a near-wall cylinder with Ref [18]

2 结果和讨论

2.1 质量比对振动响应的影响

图2给出了不同质量比下圆柱横向和流向振幅随Ur的变化情况。为与孤立圆柱(无壁面)进行对比，图2给出了m*=2和Re=100条件下两向自由振动孤立圆柱的结果。本节旨在研究质量比对近壁面圆柱涡激振动的影响，故未给出m*=10,20条件下孤立圆柱涡激振动的结果。总体上来看，随着m*的增加，振幅呈现下降的趋势，且圆柱振动也开始于更大的折合流速。当m*=2时，振动从Ur=3.0开始，而当m*=10和20时，振动开始于Ur=3.6，这与孤立圆柱的情况类似[38-39]。m*=10和20工况下横向和流向振幅差别不大。两种工况下，圆柱的横向最大振幅分别为Ymax/D=0.59和Ymax/D=0.56。圆柱的流向振动表现为两个峰值，且第一个峰值明显低于第二个峰值，这与低雷诺数条件下孤立圆柱的单峰响应不同，而与高雷诺数条件下孤立圆柱流向振动的双峰响应[40]类似。近壁圆柱的顺流向和横流向振动均出现迟滞，随着m*的增加，迟滞区域宽度先减小后增大，并在m*=10时达到最小。迟滞区域后，m*=2工况下圆柱的振幅缓慢下降，最终趋于零；而m*=10和20工况下，圆柱的振幅维持在零上，圆柱不再振动。需要说明的是：尽管个别工况下圆柱的横流向振幅大于圆柱与壁面之间的间隙，但由于圆柱振动的平衡位置发生了上移，圆柱并未与壁面发生碰撞。

图3给出了不同质量比下圆柱运动轨迹图。如图3(a)所示，当m*=2时，轨迹均为雨滴形。此时，圆柱的横向和流向振动同频。随着Ur的增加，轨迹的上部逐渐趋向于下游，由于旋涡脱落时向下游振动的圆柱与旋涡之间的作用时间增加，旋涡产生的更大激励作用，可激发较大振幅的振动[41]。而迟滞发生之后，振动轨迹上部趋向上游，旋涡脱落时圆柱与旋涡之间的作用时间减少，旋涡产生的激励作用也降低，圆柱振幅下降。当m*=10和20时，轨迹基本上保持与m*=2工况相同的趋势，但是也出现了两个明显的不同：一是由于流向振动2倍频成分占比增加，在Ur较低时，轨迹表现出类8字型；二是当振幅较大时，轨迹更加趋向于椭圆形。需要说明的时：本文所有工况的圆柱振动轨迹在上部均为顺时针方向，即圆柱由上游向下游通过上顶点。出于简洁，图3仅标记了m*=4,Ur=4工况的轨迹方向，其余工况的轨迹方向与之相同。

图2 不同质量比条件下圆柱横向和流向振幅随折合流速的变化(VIV为m*=2的孤立圆柱)Fig.2 Variation of the transverse and in-line vibration amplitudes with the reduced velocity at different mass ratios (VIV represents the results of an isolated cylinderwith m*=2)

图3 不同质量比条件下近壁圆柱运动轨迹随Ur的变化Fig.3 Variation of the vibration trajectoryat different mass ratios

2.2 质量比对脱涡频率的影响

图4给出了不同质量比下无量纲脱涡频率(St)随Ur的变化情况。当圆柱的振动开始后，St随着Ur的增加逐渐下降，直到迟滞时发生跳跃。m*=2工况下，尾流中旋涡脱落一直存在，但由于脱涡模式的变化，迟滞后的St比迟滞前的明显较大，且稳定在St=0.15附近。m*=10和20工况下，迟滞后尾流不脱落旋涡。此外，对比不同质量比下的St值可以发现，在相同的折合流速下，m*较大时，圆柱的脱涡频率较高。当脱涡频率(圆柱振动频率)与圆柱自然频率相近时，锁定现象出现。通过与自然频率对比，在m*=2工况下，锁定出现在Ur=3.0～6.5，而m*=10和20工况下，锁定分别出现在Ur= 3.5～7.3和Ur=3.5～7.6上。可见，随着m*的增加，圆柱涡激振动的锁定区间宽度增加。

图4 不同质量比条件下圆柱脱涡频率随Ur的变化(VIV为m*=2的孤立圆柱)Fig.4 Variation of the vortex shedding frequency with the reduced velocity at different mass ratios (VIV represents the results of an isolated cylinder with m*=2)

2.3 质量比对尾流模式的影响

图5给出了不同质量比下圆柱尾流随Ur的变化情况。由于m*=10和20工况下圆柱尾流类似，这里仅给出m*=2和10的情况。当Ur较低时，振幅较小，圆柱下侧的剪切层出现明显摆动，但并没有旋涡脱落，此时的尾流为1S模式，如图5(ai)和图5(bi)所示。随着Ur增加，振幅增加，圆柱下侧剪切层与壁面边界层发生相互作用，下侧开始发生旋涡脱落，此时圆柱下侧脱落的旋涡迅速分裂为两个，形成为 C+S(C指两个同相旋转的旋涡)模式[42]，如图5(aii)～图5(aiii)和图5(bii)～图5(biv)所示。但是由于下侧剪切层受到壁面的抑制，脱落的旋涡仅存在于较短距离的下游。随着Ur继续增大，振幅进一步增加，剪切层的摆动进一步增强，圆柱下侧有一个旋涡脱落，此时的尾流为2S模式[43]，如图5(aiv)～图5(aviii)和图5(bv)～图5(bvii)所示。当m*=2且Ur=5.8(Ur减小)和6.4时，由于边界层完全由间隙间通过，从下侧剪切层脱落的漩涡被显著拉长。迟滞发生后，m*=2工况下，振幅随Ur的增加而下降，剪切层的摆动逐渐减弱，圆柱下侧剪切层也不再脱落旋涡，尾流呈1S模式，如图5(aix)所示；m*=10工况下，圆柱振幅迅速下降到零，圆柱的上下剪切层变不摆动，形成稳定尾流，如图5(bix)所示。

图5 不同质量比和折合流速条件下圆柱尾涡的变化Fig.5 Wake patterns at different mass ratios and reduced velocities

2.4 迟滞现象的机理解释

在研究的质量比范围内，近壁面圆柱涡激振动均出现了迟滞现象。本节尝试对迟滞现象给出机理解释。如图6所示，增折合流速时，顺流向平衡位置偏移逐渐增大，在迟滞发生时出现小幅跳跃。另一方面，迟滞发生前，横流向平衡位置偏移随着折合流速的增加逐渐增加，然而迟滞发生时，平衡位置跳跃到距离壁面更近的位置上，圆柱与壁面之间的间隙明显变小。而后，随着折合流速的增大，平衡位置偏移逐渐减小。而减折合流速时，圆柱顺流向和横流向平衡位置偏移呈现出与上述相反的趋势。

图6 当m*=10时圆柱横向和流向平衡位置偏移随Ur的变化Fig.6 Variation of the transverse and in-line shifts of the balanced position with the reduced velocity at m*=10

横流向振动平衡位置随折合流速偏移的这种变化，直接导致了壁面边界层与圆柱及其剪切层相互作用模式的不同。图7给出了m*=10和Ur=7.2工况下增、减折合流速时尾流在一个振动周期内的变化情况。如图7(a)所示，增折合流速时，圆柱可以达到距离壁面更远的地方，当圆柱从最高处向壁面移动时，圆柱下侧剪切层可以充分发展。由于壁面之间的距离较大，间隙流也相对较强，壁面边界层通过圆柱下方。在壁面边界层的推动下，圆柱下侧剪切层包裹圆柱的下侧以及底部，使得圆柱下侧的压强更低，进而促进了圆柱向壁面运动。当圆柱距离壁面较近时，壁面边界层重附着于圆柱的上侧，与圆柱上侧剪切层相融合，使得圆柱上侧脱落的旋涡更强。由于壁面边界层给圆柱上侧剪切层提供涡量，使得圆柱上侧的压强更低。因此，当圆柱远离壁面时，圆柱受到了更强的激励作用。由于以上两种激励作用的影响，圆柱发生较大幅度的振动。如图7(b)所示，减折合流速时，壁面边界层不再与圆柱发生直接的相互作用，仅从圆柱下侧通过，抑制了圆柱下侧剪切层的发展，圆柱下侧剪切层和壁面边界层稳定不摆动，而上侧剪切层仅在较远下游出现摆动，这就解释了为什么圆柱的振动较小。

综合以上分析可知，增折合流速时，壁面边界层会在圆柱向下运动时通过圆柱与壁面之间的间隙，而在圆柱向上运动时重复着于圆柱上侧表面上，这种不稳定的相互作用模式对圆柱振动起到促进作用。而减折合流速时，壁面边界层与振动的圆柱之间没有直接作用，并且限制了圆柱下侧剪切层的发展，对圆柱振动起到抑制作用。总结来说，圆柱振动的迟滞现象与壁面边界层在增、减折合流速条件下的双稳态有关。需要说明的是，该迟滞现象仅出现在圆柱距离壁面较近的工况下。Chen等指出，当圆柱距离固壁较远时，壁面边界层与运动圆柱之间的双稳态现象消失，相应地，该迟滞现象也不再出现。进一步增大圆柱与壁面的距离(如壁面无限远的孤立圆柱)，在特定条件下圆柱振动的迟滞现象又重新出现，但孤立圆柱振动迟滞的成因与尾涡模式切换等有关，这与近壁圆柱振动迟滞的物理机制不同。

图7 m*=10和Ur=7.2条件下一个振动周期内圆柱尾涡的变化Fig.7 Vortex shedding process in a vibration cycle at m*=10 and Ur=7.2

3 结 论

本文研究了质量比对近壁面两向自由度圆柱涡激振动响应、脱涡频率以及尾流模式的影响，解释了近壁面涡激振动迟滞现象的机理。研究发现，质量比对近壁面圆柱涡激振动有着显著的影响。随着质量比的增加，圆柱的振动开始于更高的折合流速上，而且圆柱的振幅也更小一些。受壁面边界层影响，圆柱流向和横向的振动频率相等。因此，在绝大多数情况下，圆柱振动轨迹为雨滴形。但当m*=10和20时，流向振动的两倍频成分占比增加，在折合流速较小时，呈现出类8字形轨迹。随着质量比的增加，运动轨迹更加趋向于椭圆型。此外，还发现当振动轨迹上部偏向下游时，振动圆柱与脱落旋涡的作用时间更长，旋涡对圆柱振动的促进作用更大，而轨迹上部偏向上游时则恰好相反。

当振幅较小时，各质量比条件下，圆柱下侧剪切层出现明显摆动，但并无旋涡形成。对应地，尾流出现了1S模式。随着振幅的增加，圆柱下侧剪切层脱落两个旋涡，形成了C+S模式。振幅更大时，圆柱下层剪切层与壁面产生强烈的相互作用，圆柱下侧剪切层仅脱落一个旋涡，对应的尾流为2S模式。在迟滞区域之后，m*=2工况下，尾流为1S模式，而m*=10和20工况下，尾流为稳定模式。对增、减折合流速条件下的尾流模式分析发现，壁面边界层重复着的双稳态性是导致近壁面圆柱涡激振动出现迟滞的原因。当增折合流速时，壁面边界层周期性地重复着于圆柱上侧，对圆柱上侧剪切层的发展起到了显著的促进作用，并进而激励圆柱产生更大的振幅；当减折合流速时，壁面边界层从圆柱的下侧通过，抑制了下侧剪切层的发展和圆柱振动。

Vol．41 No．12 2022