刘旭菲， 花 阳， 卫昱含， 及春宁

(1. 浙江水利水电学院 水利与环境工程学院， 杭州 310018； 2. 上海勘测设计研究院有限公司， 上海 200335；3. 天津大学 水利工程仿真与安全国家重点实验室， 天津 300072)

海上油气的集疏运主要通过管道运输实现，其破坏后，修复难，经济成本和生态代价高[1]。造成其破坏的原因为悬空管道在近底水流作用下的涡激振动。当水流流经管道时，从管道的上下两侧交替泄放旋涡，产生周期性或拟周期性的升、阻力，导致管道发生弯曲振动[2-4]。在一定的流速范围内，结构的振动将对旋涡脱落产生控制作用，使得旋涡的脱落频率偏离Strauhal频率，与管道振动的固有频率相接近，发生“锁定”，导致管道大幅振动，进而诱发疲劳破坏[5-6]。因此，对于近壁面海底管道的涡激振动现象及其诱发机理进行研究，具有十分重要的工程价值和理论意义。

近壁面圆柱由于受壁面和壁面边界层流动的影响，其涡激振动响应与孤立圆柱明显不同。Taneda[7]开展了试验研究，探讨了低雷诺数(Re=140)条件下，圆柱与壁面之间的间隙G对圆柱尾涡的影响。此后，Bearman等[8-9]试验研究了高雷诺数(Re=2×104～2×105)内近壁圆柱的振动响应，当间隙比G/D小于临界值0.3时，振动响应被明显抑制；大于临界值后，振动响应逐渐恢复。然而，Jacobsen等[10]的试验研究结果表明，小间隙比工况下，即便圆柱泄涡被抑制，但圆柱仍发生微小幅值的振动。同样，Wang等[11-12]发现即使间隙比很小，圆柱下部泄涡被抑制，但上部仍泄放单列涡街，圆柱发生微幅振动。Tsahalis[13]的试验研究发现：在近壁影响下，顺流向振幅接近横流向振幅，圆柱的振动轨迹为椭圆形，这与无壁面圆柱涡激振动的“8”字形轨迹明显不同。Fu等[14]发现在近壁面的影响下，圆柱涡激振动的振动频率相较孤立圆柱工况有所增加，与速度同相位的升力值出现增大。Yang等[15-17]通过系列试验发现，随着间隙比的减小，涡激振动锁定区间的宽度不断增大，圆柱横流向振幅也逐渐增加。Nielsen等[18]考察了圆柱初始下垂对涡激振动的影响，并探讨了大长细比圆柱涡激振动的多模态特征。Li等[19-20]试验研究了水流和波浪作用下弹性圆柱的涡激振动，同样发现了高折合流速条件下的二阶振动模态。娄敏等[21-22]进一步分析了管道内流流速、管内压强等参数对涡激振动特性的影响。

在数值模拟研究方面，Zhao等[23]开展了间隙比为0.002和0.3条件下近壁面圆柱的二维数值模拟研究，重点讨论了壁面碰撞反弹对振动响应的影响。Jiang等[24]发现近壁圆柱尾涡的稳定性随间隙比减小而增强。Chern等[25-26]分别从壁面剪切力和圆柱升阻力等方面对近壁面圆柱涡激振动进行了数值模拟研究。

综合以上研究可以发现，由于近壁圆柱涡激振动问题相对复杂，与无壁面工况相比，涉及间隙比和边界层厚度的影响。另一方面，已有的近壁圆柱涡激振动数值模拟研究多针对弹性支撑的二维刚性圆柱，近壁弹性圆柱涡激振动的三维数值模拟研究开展得较少。本文对不同间隙比和边界层厚度下的近壁面弹性圆柱涡激振动开展了三维数值模拟研究，以期对以下两个问题作出解答。① 不同间隙比和边界层厚度下流场旋涡结构有何变化。② 圆柱振动振幅、振动频率、振动同步性与振动轨迹随圆柱和边界层相对位置改变发生怎样的变化。

1 数值模拟方法

流固耦合的数值模拟方法采用浸入边界法(immersed boundary method, IBM)。在IBM方法框架下，流体域通过固定的正交网格进行空间离散，而可移动的固体采用离散的拉格朗日点(浸入边界点)表示。流固之间的耦合作用通过在动量守恒方程中引入额外的附加体积力项来实现。由于通常情况下正交网格与可移动物面并不重合，故需要应用插值函数进行流固边界上的信息交换。基于IBM方法的流固耦合控制方程表示如下

(1)

∇·u=0

(2)

式中：u为速度；t为时间；p为压强；ρ为密度；ν为运动黏滞系数；∇为梯度算子；f为基于浸入边界法的附加体积力矢量，表示流固耦合边界条件，具体可参考文献[27]。

时间离散后的控制方程守恒形式如下

(3)

∇·un+1=0

(4)

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

(5)

针对传统浸入边界法施加边界条件精度不高的问题，Ji等提出了的嵌入式迭代的浸入边界法，将浸入边界法嵌入到压强泊松方程的迭代求解中，利用压强的中间解比初始值更接近真实值的特点，迭代修正附加体积力，在不显著增加计算耗时的前提下，显著提高了整个算法的求解精度。详细求解过程参见Ji等的研究。

本文采用三维二结点梁单元模拟管道弯曲振动。每个结点有6个自由度，即3个位移分量和3个转动分量。将管道进行空间离散，得到矩阵形式有限元控制方程为

(6)

式中：M，C和K分别为质量、阻尼和刚度矩阵；x(t)和F(t)分别为位移和流体力向量。阻尼矩阵采用瑞利阻尼形式，表示为

C=αM+βK

(7)

式中，α和β分别为与质量矩阵和刚度矩阵成比例的阻尼参数。

2 数值模型

2.1 问题描述

本文对近壁面圆柱涡激振动进行了三维直接数值模拟研究，弹性圆柱两端为铰接。数值模拟相关参数设置如下：圆柱直径为D；根据实际工程中海底管道的最大悬跨长度，设置圆柱展向长度L=25D；圆柱质量比m*=m/mf=6，m为圆柱的质量，mf为圆柱排开水的质量；圆柱阻尼系数分别为α=2.4×10-3和β=2.7×10-4，泊松比μ=0.3；间隙比G/D=0.6和1.0；雷诺数Re=U∞D/v=350，U∞为来流流速。考虑到亚临界雷诺数范围内(300.0

2.2 模型设置

本文的计算域与边界条件设置，如图1所示。来流为均匀来流，壁面边界层自由发展形成。通过调整圆柱距上游入口的距离，来生成不同厚度的边界层(δ/D=0，0.7和1.6)。在无边界层(δ/D=0)和薄边界层(δ/D=0.7)工况中，计算域大小为51D(X)×40D(Y)×25D(Z)，圆柱距速度入口7D(X1)；在厚边界层(δ/D=1.6)工况中，计算域大小为80D×40D×25D，圆柱距速度入口36D。圆柱周围16D(X)×6D(Y)的范围内采用正交均匀网格，网格尺寸为Δx=Δy=D/32，加密区以外采用正交渐变网格。计算域在展向上采用均匀网格，对于δ/D=0和0.7的工况，计算域沿展向共划分为256层网格，对于δ/D=1.6的工况，计算域沿展向共划分为384层网格。

图1 计算域大小与边界条件设置Fig.1 The computational domain size and boundary conditions

计算域的速度入口采用狄利克雷型边界条件，u=U，v=w=0。速度出口采用诺依曼型边界条件，∂u/∂x=∂v/∂x=∂w/∂x=0。上边界取自由滑移边界条件，∂u/∂y=0，v=0，∂w/∂y=0。对于无边界层工况，壁面边界设置为自由滑移边界条件，∂u/∂y=0，v=0，∂w/∂y=0；对于有边界层工况，壁面边界设置为不可滑移边界条件，u=v=w=0。展向边界为周期边界，圆柱表面设置为不可滑移边界条件。速度出口压强为零，上、下边界和速度入口采用诺依曼型边界条件，即压强的法向梯度为零。

2.3 网格无关性和计算精度验证

首先对网格无关性进行验证。针对G/D=0.6，δ/D=0和Re=350的工况，选取两种网格，即Δx=Δy=D/32(粗网格)和Δx=Δy=D/48(细网格)，进行数值模拟，结果如表1所示。可以看到两组数值结果之间的差异可以忽略不计，最大误差仅为0.62%。故后续模拟采用计算效率更高的粗网格。

表1 网格无关性验证Tab.1 Grid independence verification

进一步验证本文数值模拟方法的准确性。对Re=350的孤立圆柱绕流进行三维数值模拟，并于已有结果进行对比，如表2所示。可以看出本文的计算结果与Williamson等[28]的试验结果和Jiang等[29]的数值结果吻合良好，证明了本文数值模拟方法的准确性。

表2 Re=350条件下孤立圆柱绕流水动力系数对比Tab.2 Comparison of hydrodynamic coefficients for flow around an isolated cylinder at Re = 350

3 结果和讨论

为了全面考察圆柱和边界层相对位置对流场特性和振动特性的影响，本文精心设计了5种工况，分别采用两种间隙比(G/D=0.6和1.0)和3种边界层厚度(δ/D=0，0.7和1.6)。本文采用无量纲参数δ/G定量表示圆柱和边界层相对位置，其值越大，表明圆柱越浸没入边界层中。5种工况对应的δ/G，如表3所示。G/D=1.0，δ/D=0和G/D=1.0，δ/D=0.7工况的圆柱完全处于边界层外；G/D=0.6，δ/D=0.7和G/D=1.0，δ/D=1.6工况的圆柱部分处于边界层中；G/D=0.6，δ/D=1.6工况的圆柱完全处于边界层中。

表3 圆柱和边界层相对位置Tab.3 The relative postion of the cylinder in the boundary layer

3.1 流场特性

瞬时流场和旋涡结构的空间分布，如图2所示。旋涡结构采用λ2等值面进行可视化。λ2为张量S2+Ω2的第二特征根，其中，S和Ω分别为速度梯度张量∇u的对称部分和不对称部分。整体来看，流场中尾涡沿圆柱展向的分布表现出三维特性，即圆柱两端泄涡窄，跨中泄涡宽，展向旋涡成“弓”形分布，这与圆柱两端振幅小，跨中振幅大有关。这一点尤其在大间隙比和薄边界层工况下(见图2(c)和2(d))更为明显。与之相应，流向涡也表现出类似的特征，并且，随着边界层厚度的增大，流向涡的数量减少，强度降低。

(a) G/D=0.6, δ/D=0.7

(b) G/D=0.6, δ/D=1.6

(c) G/D=1.0, δ/D=0

(d) G/D=1.0, δ/D=0.7

(e) G/D=1.0, δ/D=1.6图2 不同间隙比和边界层厚度下的流场涡核分布图Fig.2 The vortex-shedding patterns with different gap ratio and boundary layer thickness

跨中位置(Z/D=12.5)截面的瞬时展向涡量图，如图3所示。整体来看，圆柱上下侧泄涡不对称，尾涡出现较明显的向上偏斜。

(a) G/D=0.6, δ/D=0.7

(b) G/D=0.6, δ/D=1.6

(c) G/D=1.0, δ/D=0

(d) G/D=1.0, δ/D=0.7

(e) G/D=1.0, δ/D=1.6图3 不同间隙比和边界层厚度下的流场瞬时涡量图Fig.3 The contours of instantaneous vorticity with different gap ratio and boundary layer thickness

当G/D=0.6时，壁面边界层与圆柱上部剪切层的耦合作用明显，壁面边界层向上倾斜，直接与从圆柱上部脱落的旋涡合并。而下部旋涡的强度和尺寸均较小，泄涡受到抑制。在薄边界层工况(δ/D=0.7)中，圆柱上下侧交替泄涡，并略微向上偏斜，形成2S泄涡模式。而在厚边界层工况(δ/D=1.6)中，圆柱下部剪切层脱落较小的的逆时针旋涡，并在向下游移动的过程中很快消散，圆柱上部泄涡与壁面平行，圆柱下游的壁面附近涡量较低。

当G/D=1.0时，圆柱泄涡宽度比G/D=0.6工况的更大，这与圆柱更大的振幅有关。在无边界层工况(δ/D=0)中，一个周期内圆柱上下侧各泄放一对旋涡，形成2P泄涡模式，上侧的涡对倾斜向上，而下侧的涡对由于受到壁面的限制与壁面平行。由于下侧涡对中的顺时针旋转旋涡(深灰色)强度较低，并很快耗散，因此泄涡模式逐渐变为P+S模式。随着边界层厚度的增加，圆柱的下侧泄涡受到抑制，圆柱上部泄涡与壁面边界层相互融合，形成倾斜向上的2S泄涡模式。

3.2 振动响应特性

分别作圆柱顺流向和横流向振动位移时程图，如图4所示。由于本文弹性圆柱的长细比较小，各工况下，近壁圆柱的振动响应均表现为一阶振型。圆柱跨中位置处的振幅最大，受圆柱端部支撑条件的限制，圆柱两端的振幅为零。此外，由于近壁面的影响，顺流向振动频率等于横流向振动频率，这与无壁面孤立圆柱涡激振动的“顺流向振动频率是横流向振动频率两倍”的结论不同。

(a) G/D=0.6, δ/D=0.7

(b) G/D=0.6, δ/D=1.6

(c) G/D=1.0, δ/D=0

(d) G/D=1.0, δ/D=0.7

(e) G/D=1.0, δ/D=1.6

图4 不同间隙比和边界层厚度下圆柱振动位移时程图Fig.4 The time history of the vibration amplitude of the cylinder with different gap ratio and boundary layer thickness

跨中截面振幅的最大值和均方根值，如表4～表7所示。当近壁圆柱处于边界层外时(δ/G<1.0)，在间隙比不变的条件下，增加边界层厚度，可使圆柱顺流向振幅增大，而横流向振幅减小。例如：当G/D=1.0时，随着δ/D从0增加到0.7，Ax,max由0.057增大到0.138，Ay,max由0.621减小到0.6。升阻力的均方根值也观察到相似的变化。其原因为：一方面随着边界层厚度的增加，圆柱下部的泄涡受到抑制，因此圆柱所受的升力减小，横流向振幅减小；另一方面，由于边界层中顺流向速度存在梯度，振动圆柱在较高位置处阻力较大，而较低位置处阻力较小。这个随位置变化的阻力叠加上圆柱泄涡产生的周期变化的阻力，导致了圆柱阻力变幅增大，顺流向振幅增大。Chen等[30]的二维近壁面两自由度圆柱涡激振动结果中也观察到了类似的现象。

表4 圆柱跨中位置顺流向最大振幅Tab.4 The maximum in-line vibration amplitude of the cylinder at the mid-span position

表5 圆柱跨中位置横流向最大振幅Tab.5 The maximum transverse vibration amplitude of the cylinder at the mid-span position

表6 圆柱跨中位置顺流向均方根振幅Tab.6 The root mean square of in-line vibration amplitude of the cylinder at the mid-span position

表7 圆柱跨中位置横流向均方根振幅Tab.7 The root mean square of transverse vibration amplitude of the cylinder at the mid-span position

然而，当圆柱由处于边界层外过渡到部分处于边界层中时(1.0<δ/G<2.0)，在间隙比不变的条件下，增加边界层的厚度，圆柱的顺流向振幅进一步增大，而横流向振幅则小幅增加。例如：当G/D=1.0时，随着δ/D从0.7增加到1.6，Ax,max由0.138增大到0.146，Ay,max由0.6小幅增大到0.606。产生这样结果的原因与顺流向振动对横流向振动的促进作用有关。Williamson等[31]指出：当圆柱发生横流向涡激振动时，无量纲最大振幅约为1.0；而圆柱发生两自由度涡激振动时，由于顺流向振动的促进作用，圆柱的横流向无量纲振幅可达1.5。

最后，当圆柱由部分没入边界层到完全没入边界层时，在间隙比不变的条件下，圆柱的顺流向和横流向振幅随着边界层厚度增加均出现减小的趋势。例如：当G/D=0.6时，随着δ/D从0.7增加到1.6，Ax,max由0.152减小到0.046，Ay,max由0.539减小到0.396。这是由于随着近壁圆柱完全浸没入边界层中(G/D=0.6和δ/D=1.6工况)，圆柱所在位置处的来流速度明显降低，导致了圆柱所受升、阻力大幅减小。

可见，近壁圆柱的振动幅值与圆柱和边界层之间的相对位置密切相关。当圆柱处于边界层之外时，随着边界层厚度的增大，顺流向振幅增大，而横流向振幅减小；当圆柱由在边界层外到部分没入边界层时，顺流向振幅增大，而横流向振幅小幅增大；当圆柱由部分处于边界层中到全部淹没在边界层中时，圆柱的顺流向和横流向振幅均大幅减小。

圆柱顺流向和横流向振动能量谱，如图5所示。整体来看，两个方向的振动能量均集中跨中部分，且具有相同的频率，表现出明显的单频特性。需要说明的是，G/D=1.0和δ/D=0工况的顺流向振动在跨中部分出现了微幅的二阶振动分量，这与此时圆柱的2P泄涡模式有关(见图3)。

(a) G/D=0.6, δ/D=0.7

(b) G/D=0.6, δ/D=1.6

(c) G/D=1.0, δ/D=0

(d) G/D=1.0, δ/D=0.7

(e) G/D=1.0, δ/D=1.6图5 不同间隙比和边界层厚度下圆柱两向振动频谱图Fig.5 The vibration frequency spectrum of cylinder with different gap ratio and boundary layer thickness

圆柱的顺流向和横流向振动频率，如表8、表9所示。可见，两个方向的振动频率均保持一致。当圆柱处于边界层外，或者部分处于边界层中的工况中，圆柱两向振动频率保持不变，均为St=0.222。然而，当圆柱完全浸没于边界层内时(G/D=0.6和δ/D=1.6工况)，圆柱两向振动频率出现明显下降，这与降低的来流流速有关。

表8 顺流向振动频率Tab.8 The in-line vibration frequency

表9 横流向振动频率Tab.9 The transverse vibration frequency

通过Hilbert变换求得圆柱顺流向和横流向的瞬时相位(即φx和φy)，进而通过式求得圆柱两向振动的相位差φx,y。

φx,y=[pφx-qφy,mod 360°]

(8)

式中，p和q为整数，并满足fx/fy=p/q，fx和fy分别为圆柱顺流向和横流向的振动频率。在近壁面条件下，圆柱的两向振动频率基本保持一致，故取p=1，q=1。

不同间隙比和边界层厚度下圆柱两向振动相位差均值，如表10所示。在两个方向同频振动的条件下，圆柱振动轨迹的旋转方向可以通过两向振动的相位差表征，而振动轨迹的偏斜程度可以通过两向振动的幅值和相位差共同表征。当δ/D=0，0.7时，圆柱两向振动相位差小于180°，圆柱振动表现为逆时针轨迹。而当δ/D=1.6时，圆柱顺流向和横流向振动相位差超过了180°，圆柱振动表现为顺时针轨迹。可以看出，随着间隙比的减小和边界层厚度的增加，即随着圆柱逐渐没入壁面边界层中，两向振动的相位差均值逐渐增大。这与壁面边界层中的剪切流动有关。当圆柱靠近壁面时，圆柱周围的流速较低；反之，当圆柱远离壁面时，圆柱周围的流速较高。这就造成了部分或全部浸没的圆柱在远离壁面过程中，同时向下游运动，即圆柱的振动轨迹为顺时针方向，两向振动的相位差较大。并且，圆柱没入边界层越深，两向振动的相位差越大。

表10 圆柱顺流向和横流向振动相位差均值Tab.10 The mean value of the vibration phase difference between in-line and transverse vibration (°)

沿圆柱展向等间距截取10个截面，以振动平衡位置为中心，作不同间隙比和边界层厚度下圆柱振动轨迹。圆柱的振动轨迹为雨滴形，这与无壁面孤立圆柱涡激振动的“8”字形轨迹明显不同，如图6所示。

图6 不同间隙比和边界层厚度下圆柱各截面振动轨迹图Fig.6 The XY-trajectory of each section of cylinder with different gap ratio and boundary layer thickness

当G/D=0.6，δ/D=0.7时，由于两向振动的相位差均值(171°)接近180°，振动轨迹表现为瘦长雨滴形，且下部较为尖锐，上部向来流方向偏斜。当G/D=0.6，δ/D=1.6时，相位差均值(290°)远离180°，圆柱逆时针振动，且轨迹较为圆润，上部略偏向下游。其他3个工况具有类似的结果，此处不再赘述。

4 结 论

本文应用嵌入式迭代浸入边界法对不同间隙比和边界层厚度下的近壁面弹性圆柱涡激振动开展了三维数值模拟研究。分析了各工况下的流场特性和圆柱振动响应特性。研究结果表明：

(1) 圆柱发生涡激振动时，跨中部分泄涡强度较高，尾涡较宽，与壁面边界层的相互作用较强；两端的情况则相反。

(2) 近壁圆柱的振动幅值与圆柱和边界层之间的相对位置密切相关。当圆柱处于边界层之外时，在间隙比不变条件下，随着边界层厚度的增大，顺流向振幅增大，而横流向振幅减小；当圆柱部分没入边界层时，在间隙比不变条件下，随着边界层厚度的增大，顺流向振幅增大，而横流向振幅小幅增大；当圆柱全部淹没在边界层中时，在间隙比不变条件下，随着边界层厚度的增大，圆柱的顺流向和横流向振幅均大幅减小。

(3) 圆柱顺流向和横流向振动同频，均呈现一阶振型。当圆柱位于边界层外或部分处于边界层中时，圆柱振动频率几乎不受边界层厚度和间隙比影响，频率均集中于St=0.222附近。当圆柱完全位于边界层中时，频率明显减小，集中于St=0.199附近。

(4) 在近壁面条件下，圆柱振动轨迹呈现水滴形或椭圆形，随着近壁圆柱逐渐没入边界层中，圆柱顺流向和横流向振动相位差逐渐增大，在无边界层和薄边界层工况中，圆柱振动表现为逆时针轨迹，在厚边界层工况中，圆柱振动表现出顺时针轨迹。