三角分布的海洋立管涡激振动数值分析

张 猛， 刘 冉， 赵桂峰， 王军雷

(1. 郑州大学 土木工程学院， 郑州 450001； 2. 郑州大学 化工与能源学院， 郑州 450001)

海洋立管是深海油气开发系统中连接海面作业平台和海底钻采设施的关键设备，如钻井隔水套管、输液立管等，是系统中薄弱易损的构件之一[1]。当海流经过立管结构时，在立管两侧会产生交替的旋涡脱落，柱体下游的流动分离会在其尾流中产生环流，根据贝努利方程，在柱体上下侧会出现压强差，由此产生横向力，随着旋涡脱落的交替出现，上下横向力也交替产生，从而在结构表面形成周期性的脉动作用力，很容易诱发海洋立管涡激振动[2]。若流体流速增大，旋涡脱落频率接近于结构的固有频率fn时，两者会发生共振，尾涡的非定常压力就会作用在结构上，引起结构加速运动，大大增加其动能和振幅。长时间的振动极易引起海洋立管结构疲劳破坏，会对结构安全运行构成很大威胁[3]。

实际工程中，海洋立管往往不是单独存在的，常以管群的方式出现[4],如图1所示。

图1 管群结构形式

对于典型的单圆柱钝体结构，已引起学者们的广泛关注，现阶段已对单圆柱绕流及涡激振动响应进行了大量试验与数值模拟研究[5]。而当流体流过多圆柱时其流动特性要比单圆柱复杂的多，当立管彼此相互靠近时，会发生流场干涉效应，这种情况下的激励诱发特性与单根立管是大不相同的。目前对管群干扰效应的研究大多数是针对双圆柱体系[6]。对于各种布置的双圆柱绕流，Zdravkovich等[7-9]通过实验方法研究发现根据圆柱间距将流体干扰可分为三类，当两个圆柱彼此足够靠近时会发生接近干扰，当下游圆柱部分或全部浸没在上游圆柱尾流时会发生尾流干扰，在第三类中，接近和尾流干扰被合并。黄钰期等[10]使用分块耦合求解方法来数值模拟单圆柱以及不同间距下的串列双圆柱绕流情况,研究分析了改变双圆柱中心间距对上下游圆柱的升阻力系数和脉动频率所产生的影响

目前国内外研究学者对于管群中的三圆柱绕流问题还少有研究，三圆柱是一种具有代表性的结构排列组合单位，由于成三角形分布的三圆柱可以看成是对多对双圆柱的组合，所以对于三圆柱绕流现象的研究具有重要的学术意义和工程价值[11]。基于此，本文拟对成三角形布置的等直径三圆柱体系在均匀来流中的涡激振动特性进行研究，模拟分析来流约化速度、圆柱间距等参数对圆柱升阻力系数、振幅及尾流旋涡脱落模式的影响，以期为海洋立管管群布置提供参考依据。

1 结构模型

XFlow是基于格子Boltzmann粒子网格技术，用介观模型来模拟流体宏观行为，具有易于使用、高效并行、边界条件处理简单、模拟精确等特点[12]。本文采用XFlow进行海洋立管涡激振动计算，计算域如图2所示。

该流场域有5个边界条件，进出口边界，上下边界和三圆柱外表面边界。把三圆柱看成一个整体，形心到进口边界和上下边界的距离相同均设置为10D，形心到出口边界的距离设置为25D。对于三个圆柱，下游两并联圆柱固定不动，解除上游圆柱的横向约束，使其可以发生振动。上游圆柱到下游圆柱的水平净距离设为L1，下游两圆柱的垂直净距离设为L2，本文中L1与L2始终保持相等，令L1=L2=L如图2所示。

图2 计算域

由于上游圆柱是非固定的，在流场中该圆柱在横向力作用下会产生横向单自由度涡激振动，我们把它简化为质量-弹簧-阻尼系统，如图3所示。

(a) 三维模型

(b) 二维模型

为了使圆柱绕流模拟结果更精确，本文采用自适应网格划分方式，分别对圆柱周围与尾流旋涡进行捕捉加密，如图4所示。

2 数学模型

本文主要研究来流约化速度及圆柱间距比对圆柱绕流特性的影响，以海水为流场介质，假设密度为常数，则该介质为定常不可压缩黏性流体。

2.1 流动控制方程

由于流体流动要受到物理守恒定律的支配，所以需要满足质量守恒定律(连续性方程)、动量守恒定律。

(a) 圆柱周围粒子网格加密

(b) 圆柱及尾流区粒子网格捕捉加密

(1)

(2)

写成矢量形式

(3)

式中：μ为流体黏性系数；fx,fy,fz为外力分量；ux,uy,uz为来流速度分量；ρf为流体密度；P为微元体上的压强。

2.2 振动方程

单自由度质量-弹簧-阻尼系统振动方程为

(4)

(5)

(6)

(7)

式中：Fl(t)、Fd(t)为圆柱所受到的升力、阻力；Cl(t)、Cd(t)为与升力所对应的升力系数、阻力系数；ξ为阻尼比。各流体动力参数如表1所示。

表1 流体动力参数

式中：fs为漩涡脱落频率, 一般是通过对升力的傅里叶变换求得；fn为圆柱固有频率；D为圆柱的直径；v为运动黏度；U∞为来流速度；m为海洋立管质量；k为系统刚度。立管和流体参数见表2。

表2 立管和海水的相关参数[13]

由于海洋立管的特点是低阻尼，低质量比，在模拟计算时根据Carmo等[14]的结果取阻尼比ξ= 0.007，将立管和海水参数取值代入质量比表达式可以求出m*= 2.5。Lee等考虑立管张力的自振频率计算公式为[15]

(8)

将参数代入可求得fn= 0.76 Hz。

3 流体求解器的验证

本文分别对单圆柱和三圆柱进行了模拟分析，为验证单圆柱流固耦合计算的正确性，其各参数取值与前文所定义的海洋立管和海水参数相同, 模拟了单圆柱在雷诺数等于200、3 900、10 000时的阻力系数与斯特劳哈尔数，并将结果与其他文献中数据进行对比。如表3所示，其结果吻合较好，说明单圆柱数值模拟的精度较高。

表3 单圆柱绕流模拟结果对比

对于三圆柱，Gu等[25]对成三角形布置的三圆柱进行了风洞试验研究，为验证本文中三圆柱流固耦合计算的正确性，其模拟参数取值与试验保持一致，圆柱直径D=48 mm，圆柱间距比L/D=1.7，来流风速为18 m/s，空气密度ρf=1.205 kg/m3，动力黏度系数u=1.81×10-5Pa·s，模拟分析了三个圆柱的平均阻力系数与平均升力系数，并与试验结果进行对比，如表4所示，其结果吻合较好。

表4 数值模拟结果与风洞试验结果对比

4 计算结果及分析

为了对比分析在流场中三圆柱间的相互干扰强度，首先模拟了相同参数的单圆柱涡激振动。在三圆柱涡激振动数值模拟中，圆柱间距分别根据小间距比，中间距比和大间距比取L/D= 0.5、2.0、4.0。考虑到来流约化速度对圆柱绕流的影响，在每一种圆柱间距工况下又分别选取了Ur= 1、1.5、2、2.5、3、3.5、4、4.5、5、5.5、6、6.5、7、7.5、8、8.5、9等工况进行模拟计算。本文中用CD1(CL1)，CD2(CL2)，CD3(CL3)分别代表上游圆柱，下游上圆柱，下游下圆柱的阻(升)力系数。

4.1 涡激振动特性

以小间距比L/D= 0.5 mm，Ur= 4为例，图5～图7中分别给出了上游圆柱的涡激振动位移时程曲线，三个圆柱的升阻力时程曲线以及各圆柱升力频谱图。

从图5可知，在上游圆柱的前几个振动周期内，振幅逐渐增大，随后趋于稳定，体现了流体与固体的耦合作用，随着上游圆柱振幅的增大，其对周围流体的扰动也越来越大，流体的变化反过来又会作用在圆柱上，抑制圆柱的振动，所以圆柱的振幅不会一直增加，当流体与固体耦合稳定时，圆柱的振幅也会趋于稳定[26]。

图5 上游圆柱位移时程曲线

图6列举了三个圆柱的升阻力系数时程曲线。从图6可知，升力系数周期是阻力系数的两倍，对于上游圆柱升力系数脉动幅值明显大于阻力系数，由此可见升力对上游圆柱的影响更大。

将图6各圆柱升力系数经过傅里叶变换可得到图7的频谱图，从图7可知，三圆柱的旋涡脱落频率均等于0.596 Hz。对比上下游圆柱的频谱图可以发现，上游圆柱只有一个对应于fs= 0.596 Hz的频率幅值突出点，但下游两并列圆柱则有多个频率幅值突出点，主要是因为小间距比时，上下游圆柱距离太近，使得上游圆柱后方来不及形成涡脱就直接作用在下游圆柱上，干扰下游圆柱的旋涡脱落。

(a) 上游圆柱

(b) 下游上圆柱

(c) 下游下圆柱

图6 升阻力系数时程图

Fig.6 Lift-drag coefficient time history curve

如图8所示，分别列出了三圆柱不同间距时，上游圆柱振幅均方根和圆截面特征长度的比值Yrms/D与来流约化速度的关系曲线，并与单圆柱进行对比。从图8可知，当来流约化速度较小时(Ur< 3)上游圆柱振幅均趋于0，激起圆柱振动不显著，当Ur达到3.5时上游圆柱振幅均有跳跃性增加，然后随来流约化速度的增加，上游圆柱均保持在一个较高幅值振动，当Ur进一步增大时，振幅又逐渐减小。由此可以看出L/D=0.5、2.0、4.0时其上游圆柱涡激振动的锁振区间分别为3.5～7.0、3.5～7.5、3.5～6.0，最大振幅分别为Yrms/D= 0.372, 0.546, 0.470，单圆柱大约在3.5～5范围内发生锁振，最大振幅为Yrms/D= 0.472。通过对比可以发现圆柱间距比的大小对上游圆柱的锁振区间和最大振幅的影响很明显，L/D= 4.0时上游圆柱的锁振区间、最大振幅与单圆柱很接近，证明对于大间距比，上游圆柱受下游圆柱的干扰很小。L/D= 2.0时上游圆柱的锁振区间和最大幅值均明显大于单圆柱，证明中间距比，下游两并列圆柱会促进上游圆柱保持较高振幅振动。而L/D= 0.5时上游圆柱锁振区间和最大幅值均明显小于单圆柱，证明小间距比，下游两并列圆柱会抑制上游圆柱的高振幅振动。所以中间距比时海洋立管振幅最大，最不安全。

(a) 上游圆柱

(b) 下游上圆柱

(c) 下游下圆柱

图7 升力系数频谱图

Fig.7 Lift coefficient spectrum diagram

图8 上游圆柱振幅随来流约化速度的变化曲线

图9列出了不同间距比时上游圆柱涡脱频率比fs/fn与来流约化速度Ur的关系曲线，由图可以看出，L/D= 0.5和L/D= 2.0关系曲线变化规律相同，L/D= 4.0与单圆柱关系曲线变化规律相同。以L/D= 2.0和L/D= 4.0为例来分别介绍：当L/D= 2.0时，在Ur= 4～7.5范围内，频率比fs/fn≈ 1，这表明上游圆柱的涡脱频率被锁定在了结构的自振频率附近，上游圆柱与流体发生了共振。当约化速度超过共振区间时，频率比随约化流速的增大而减小，其变化规律较单圆柱有较大差异，这是由于上游圆柱无明显振动，上下游圆柱间距始终较小，圆柱间的接近效应和尾流效应的影响抑制了上游圆柱的旋涡脱落。当L/D= 4.0时，在Ur= 4.5～6范围内，频率比fs/fn≈ 1，当约化速度超过共振区间时，频率比变化规律与单圆柱相同，均随约化流速的增加而增大，因为对于大间距比，上下游圆柱始终可以产生充分发展的旋涡脱落。

通过对上游圆柱锁振区与共振区对比可以发现，某些约化速度虽处于上游圆柱的锁振区间，但其涡脱频率并未锁定在柱体自振频率附近，这与本文数值模拟所采用的质量比有关，由于柱体结构与周围流场有很强的非线性关系，可能会有附加质量对质量弹簧系统固有频率和发生涡激共振时的振动频率产生影响[27]。

图9 频率比与来流约化速度的关系曲线

4.2 流动特性

如表5所示，分别给出了L/D= 0.5、2.0、4.0，来流约化速度Ur= 1.5、7、8时在同一时刻的旋涡脱落云图，所选约化速度对应着上游圆柱涡激振动锁振与非锁振阶段。

从表5云图可知，在小间距比时(L/D= 0.5)，圆柱之间的干扰较为严重，上游圆柱后方均没有出现回流区和旋涡，是由于下游两并列圆柱限制了上游圆柱自由剪切层的卷起，使其不能产生旋涡脱落，而是与下游圆柱的分离层相互作用，只在下游圆柱后方出现规则的旋涡脱落现象。其流体流动特性主要受接近效应的影响。

在中间距比时(L/D= 2.0)，当约化速度在非锁振区范围内，在上游圆柱后方产生了旋涡，但下游圆柱抑制了其旋涡的脱落，上游圆柱后方的尾流被下游两并列圆柱的尾流夹在中间，在其间隙中变窄，可见其流体流动特性受接近效应和尾流效应的共同作用。当约化速度处于锁振区时，由于上游圆柱与流体发生了共振使其横向振动加剧，振幅增大，这时上游圆柱与下游两并列圆柱的间距比成周期性变化，当间距比达到一定值时，上游圆柱后方可产生旋涡脱落，其流动特性仅受尾流效应的影响。

在大间距比时(L/D= 4.0)，上、下游圆柱后方始终会产生充分发展的旋涡脱落的现象。在三个圆柱后方形成一个接近同步的尾流模式，来自上游圆柱的自由剪切层在下游两并列圆柱的间隙中波动并与下游圆柱的尾流发生干涉，影响了下游圆柱的旋涡脱落，因而使得下游圆柱后方的涡脱形式变得不规则。而这种作用随来流约化速度的增大而减小，当Ur= 8时，上下游圆柱的尾流互不干扰，周围的流体流动模式相对于自由流方向是不对称的，且呈现出明显的偏流现象。

表5 不同圆柱间距比和约化流速下圆柱涡量图

4.3 升阻力特性

圆柱间流体的流动特性决定着圆柱的升阻力特性。为研究来流约化速度、圆柱间距比对各圆柱升阻力系数的影响，如图10～图11分别列出了L/D=0.5、2.0、4.0时的数值模拟结果。

图10给出了不同间距比时三圆柱阻力均值与来流约化速度的关系曲线，从图中可以看出：上下游圆柱的阻力均值在不同间距比时随流体约化速度的变化趋势相同。对于上游圆柱，当Ur< 3.5时，其阻力均值随流体约化速度的增加无明显变化，当Ur= 3.5时发生突变，阻力均值突然大幅增加致最大值(CD，mean=0.640, 0.593, 0.566)，随后又随流体约化速度的增加而减小。对应图8，我们可以发现在Ur= 3.5时上游圆柱阻力均值的突变与振幅的突变有关。对于下游两并列圆柱的阻力均值除了在Ur< 3.5区间内，其余均小于上游圆柱，这是由于当上游圆柱进入锁振区间后振幅增大，涡脱频率增加，其尾流更大程度上包裹了下游圆柱。

图11给出了不同间距比时三圆柱升力均方根与来流约化速度的关系曲线，从图11可知：① 上下游圆柱的升力均方根随来流约化速度的变化大体一致，Ur< 3.5范围内各圆柱的升力均方根均保持较小值，在Ur= 3.5时同阻力均值一样由于上游圆柱进入锁振区而跳跃增加致最大值(CL,rms= 0.412, 0.430, 0.382)，下游圆柱在Ur=4.0时达到最大值，随着来流约化速度的继续增加，各圆柱升力均方根逐渐减小;② 下游两并列圆柱的升力均方根普遍大于上游圆柱，是由于上游圆柱后方的尾流或旋涡作用在下游圆柱表面，会对其产生一个附加力;③ 圆柱间距比足够大时，上下游圆柱升力均方根仅在锁振区间存在较小的差值，因为对于非锁振区间，上游圆柱无明显振动，下游两并列圆柱间距较大，由表5可知，上游圆柱后方产生的旋涡直接从下游圆柱间隙中流过，不能充分作用在下游圆柱表面。对于锁振区间，由于上游圆柱发生大幅振动，其后方产生的旋涡会周期性的作用在下游圆柱表面，且由于上下游圆柱间距比较大，上游圆柱产生的旋涡需经过一段长距离才能作用于下游圆柱，其能量有较大的损耗。由此可以得出，当圆柱达到大间距比时，随圆柱间距比的增大，上下游圆柱的升力系数差值会逐渐减小;④ 在小间距比的锁振区间及中间距比的超锁振区间，下游两并列圆柱的升力均方根出现明显差异，这是由于上游圆柱后方的尾流是不对称的。

(a) L/D = 0.5

(b) L/D = 2.0

(c) L/D = 4.0

5 结 论

对阻尼比ξ= 0.007，质量比m*= 2.5，三角形分布的海洋立管在间距比为L/D= 0.5、2.0、4.0，来流约化速度为Ur= 1～9进行了涡激振动数值模拟，并与单圆柱涡激振动进行了对比，将结果归纳如下：

(1) 海洋立管间距比的大小对上游立管的锁振区间和最大振幅有很大影响，中间距比时其锁振区间和最大振幅最大，最不安全。

(a) L/D = 0.5

(b) L/D = 2.0

(c) L/D = 4.0

(2) 由三圆柱涡量云图对比可知，立管间距比较大或较小时，立管间流体的流动模式分别受尾流效应和接近效应的影响，中间距比时受两种效应的共同作用。

(3) 由于上游立管的尾流或旋涡会对下游立管产生附加力，使得上游立管的升力系数普遍低于下游立管，且当立管达到大间距比时，随立管间距比的增大，上下游立管的升力系数差值会逐渐减小。在小间距比的锁振区间和中间距比的超锁振区间，下游两并列立管表现出明显的不对称性，这是由于上游立管后方流体流动的不对称性，尾流区出现了偏流现象。