振荡流作用下圆柱体结构的水动力特性研究

2021-07-01杨冲霄袁昱超薛鸿祥唐文勇

海洋工程 2021年3期

杨冲霄，袁昱超，薛鸿祥，唐文勇

(上海交通大学海洋工程国家重点实验室，上海 200240)

随着我国海洋油气资源开发逐渐向深海迈进，海洋工程结构物的安全性能越来越得到重视。海洋立管系统作为连接水面浮式装置和海底设备的导管，在复杂的载荷作用下，极易出现碰撞、振动及疲劳等破坏形式，是深海工程装备中最薄弱的环节[1]。在来流作用下，立管的交替泄涡现象产生振荡的流体力，引起结构的涡激振动[2]。早期国外学者对该问题开展了大量的试验研究[3-5]。其中Gopalkrishnan[3]对Re为10 000下的圆柱受迫振动进行了模型试验，并首次构建了受迫振动流体力系数库。在数值方法中，随着计算机技术的发展，计算流体力学方法(computational fluid dynamics，简称CFD)应运而生，为了提高计算效率，将立管简化为二维圆柱是较为常见的研究手段。Muhamad和Krish[6]使用RANS k-ω湍流模型对Re=10 000来流下的圆柱进行单自由度受迫振动数值模拟，结果表明漩涡尾迹对圆柱振动频率较为敏感，当振动频率接近Strouhal频率时，会发生漩涡的相位变换。王亚非[7]数值模拟了双自由度弹性支撑圆柱体的涡激振动问题，结果表明振荡流扩大了自激振动的锁定范围，并且使得顺流向振动幅值大大增加，甚至超过了横向所能激发的最大幅值。付博文等[8]基于切片理论，通过使用径向基函数法作为OpenFOAM中的动网格策略，模拟了长细比为1 000的柔性立管在横流向和顺流向的振动，数值模拟重现了高阶主控模态及主控模态的频繁变换等大长细比柔性立管的涡激振动特性。王凯鹏[9]通过对均匀来流和剪切来流两种来流形式下的静止圆柱、横向受迫振动圆柱绕流及圆柱涡激振动问题展开数值模拟，首次系统地分析了来流形式对圆柱绕流和圆柱涡激振动问题的影响，并指出来流形式是很重要的一个影响因素。总体来说目前的研究主要集中在定常流。

实际生产过程中，海洋浮式结构物在风、浪、流的联合作用下会带动立管在水中往复运动，此时立管的遭遇流场可等效为振荡流场[10]。相比定常流，振荡流场引起的圆柱绕流问题更为复杂，需要更加深入的研究，目前已有学者获得了相关成果。Zhao等[11]模拟了振荡流和均匀流共同作用下圆柱的涡激振动，探讨了流动比例a对圆柱的响应，并发现了在a=0.8，Vr=7时，漩涡在一个振荡流周期中经历2S、2P和2T三种泄涡模式。邓跃[12]对低雷诺数时单自由度弹性支撑圆柱在均匀流和振荡流共同作用下的受迫和自激振动进行了数值模拟，结果发现当有振荡流参与时，结构的锁定范围和振动幅值等都有明显的变化。邓迪等[13]采用OpenFOAM开源软件对在静水中做周期性振荡运动的二维刚性圆柱涡激振动进行数值模拟，发现圆柱的横向振动加剧了升力系数的变化，使得泄涡方向和圆柱表面漩涡分离点的位置发生了明显变化。

目前对于振荡流的研究主要集中在低雷诺数或小KC数工况，而已有认知表明，不同雷诺数及KC数对圆柱绕流特性影响显著[14-16]。文中重点研究高雷诺数条件并兼顾大KC数工况。首先，对均匀流下雷诺数10 000时的圆柱受迫振动进行了数值模拟，通过与试验结果进行对比，验证了基于CFD方法复现受迫振动试验并研究流体力系数的可行性。之后开展振荡流下圆柱受迫振动数值仿真，分析不同KC数下圆柱的水动力特性及漩涡形态，并归纳出大KC数和小KC数下升力系数和漩涡发放特点。

1 数值分析方法

1.1 控制方程

考虑不可压缩流场，采用雷诺时均方法(Reynolds Averaged Navier-Stokes，简称RANS)，结合剪切应力输运模型(k-ω SST)求解N-S方程，控制方程：

(1)

(2)

k-ω SST湍流模型的输运方程：

(3)

(4)

式中：μt为涡黏性；Sij为平均速度应变率张量；τtij表示雷诺应力的涡黏性模型；σk、β*、σω、σω2均为经验参数；Pω代表交错扩散项。F1为混合函数，在近壁处采用Wilcox k-ω模型，边界层边缘和自由剪切层采用k-ε模型，中间通过F1函数实现过渡。

1.2 水动力系数识别方法

在Gopalkrishnan[3]的试验中，圆柱进行受迫振动，位移函数：

y(t)=Asin(2πf0t)

(5)

式中：A为振动幅值；f0为振动频率；t为运动时间。

由此得到圆柱运动的速度函数为：

v(t)=2πf0Acos(2πf0t)

(6)

当圆柱以式(5)进行振动时，垂直来流方向的升力可表示为:

L=L0sin(2πf0t+φ0)+Lssin(2πfst+φs)

(7)

式中：L0、Ls为升力幅值；φ0、φs为相位角；fs代表Strouhal频率，其中下标为0表示与圆柱振动频率相关，下标为s表示与Strouhal频率相关。

当圆柱处于锁定状态时，Strouhal频率消失，并且圆柱并未以频率fs在振动，即该成分不参与到流体与结构的能量传递中，故在实际分析时，升力得到简化：

L=L0sin(2πf0t+φ0)

(8)

升力幅值及相位角可通过对时历曲线进行傅里叶拟合得到：

(9)

计算可得参数：

(10)

(11)

(12)

升力幅值和相位角可表示为：

(13)

(14)

升力系数由下式计算得到：

(15)

式中：l为圆柱高度，在二维计算中取单位高度；D为圆柱直径；U为遭遇流速。

升力系数还可进一步分解得到激励力系数和惯性力系数。升力系数中和圆柱运动速度同相位的部分定义为激励力系数：

CL_V0=CL0sinφ0

(16)

激励力系数为正代表能量从流体输入到圆柱结构中，此时易发生圆柱的涡激振动，激励力系数为负代表能量从结构输出到流体中。

升力系数中和圆柱运动加速度同相位的部分定义为惯性力系数：

CL_A0=CL0(-cosφ0)

(17)

无因次振幅定义为：

A*=A/D

(18)

无因次频率定义为：

f*=f0U/D

(19)

振荡流定义如下：

u(t)=2πfBsin(2πft)

(20)

式中：B为振荡流幅值；f为振荡流频率。

在均匀流中，采用雷诺数定义流体特征：

Re=UD/ν

(21)

振荡流中，引入柯莱根—卡彭特数(Keulegan-Carpenter number，KC)表征来流：

KC=Umax/fD

(22)

Umax=2πfB

(23)

式中：Umax为振荡流最大速度。

1.3 计算网格和边界条件

计算域如图1所示。X轴平行来流方向，Y轴垂直来流方向。圆柱直径D=0.025 4 m，上下边距圆柱中心8D，水平方向50D。由于圆柱周围流场变化较为剧烈，在圆柱中心3.5D范围内进行加密，如图2所示。左端为速度入口，均匀流流速0.4 m/s；上下两端为对称边界；右端为自由出流；圆柱采用无滑移固壁条件；压力—速度耦合采用SIMPLE算法；时间项采用二阶隐式积分方法；对流项采用二阶迎风离散格式。

图1 计算域Fig. 1 Computational domain

图2 局部加密Fig. 2 Partial encryption

2 均匀流工况试验验证

2.1 网格收敛性验证

首先对定常流速0.4 m/s、Re=10 000下的静止圆柱进行了数值模拟，给出了普通网格、加密网格及其他学者计算结果的对比，如表1所示。

表1 两套网格参数Tab. 1 Two sets of grid parameters

根据对比可知，两套网格的计算结果均与试验数据对应较好，在综合考虑计算精度和计算成本的前提下，选择普通网格完成振荡流流动计算。

2.2 流体力系数对比

计算了均匀流振幅比为0.3时各无因次频率下的流体力系数，并与Gopalkrishnan[3]的试验数据进行了对比，升力系数、激励力系数、惯性力系数及相位角结果如图3所示。

图3 流体力系数对比Fig. 3 Comparison of fluid force coefficient

由图3可知，随着无因次频率的增加，升力系数幅值呈现逐渐增加的趋势，文中模拟结果在趋势和数值上和试验结果吻合较好。激励力系数变化较为复杂，在无因次频率0.17处出现正峰值，在0.21处出现负峰值，模拟结果基本可复现这一现象，但数值上略有差距。惯性力系数变化与无因次频率呈现负相关关系，文中CFD模拟结果较好的吻合了试验数据。相位角在低频率处为负值，在频率0.15附近发生相位突变现象，之后随着无因次频率增加逐渐趋于0，CFD模拟结果也基本符合这些特征。但在低频率处，文中模拟结果与试验存在一定误差，这可能是因为Gopalkrishnan所进行的试验结果包含了复杂的三维效应，故二维模型暂未全面反映三维试验结果。文献[18]针对二维数值模拟与试验误差也提出了相似的分析。

综上所述，文中所采用的网格划分和CFD设置方法复现了圆柱受迫振动的试验结果，在水动力系数模拟上具有较高的准确性，可运用到振荡流中开展研究。

2.3 均匀流漩涡形态分析

均匀流工况漩涡发放形态如图4所示，选取一个脱落周期的漩涡发放结果。图4(a)时，圆柱向上方运动，尾涡在圆柱下流向产生，并逐渐向后延伸。当圆柱向下运动到图4(b)位置时，漩涡从圆柱的下尾涡末端脱落，此时上尾涡的末端也分离出了即将脱落的漩涡。图4(c)时，漩涡从上尾涡末端分离。脱落周期持续0.37 s，与圆柱此时的振动周期是吻合的。可以看到，漩涡在一个周期内呈现上下交替脱落现象，为典型的2S脱落模式。

图4 均匀流下圆柱漩涡脱落(A*=0.5，f*=0.17)Fig. 4 Vortex shedding under uniform flow(A*=0.5，f*=0.17)

3 振荡流结果与分析

3.1 工况设置

为了与前文均匀流形成对比，振荡流流速幅值定为0.4 m/s，振荡流可表示为：

v=0.4×sin(2πft)

(24)

基于第2节的研究基础，通过改变振荡流振荡周期，研究不同KC数下的二维圆柱受迫振动特性。其中，振荡流工况设置情况如表2所示。为研究圆柱振动幅度和振动频率对计算结果的影响，同一振幅比下选择8组无因次频率进行计算，无因次振幅比设置6组，共计算240个振荡流工况。

表2 振荡流工况设置Tab. 2 Setting of oscillation flow condition

3.2 流体力系数分析

提取0.50D和1.00D两组振幅比的计算结果进行分析，升力系数幅值变化情况见图5(a)。在低振动频率时，升力系数呈现较小的值，随着振动频率增加，升力系数幅值也逐渐增加。在该工况下，不同KC数对升力系数幅值影响不大。激励力系数变化情况见图5(b)。激励力系数变化较为复杂，且对KC数的敏感性较高。此时不同KC数下的激励力系数随无因次频率变化基本呈现先增加后减小的趋势，峰值区出现在0.18附近，相比均匀流的0.17略有变化，符合文献[12]的研究结论。在振动频率0.15～0.20之间，KC数对激励力系数影响较大，该区间也属于均匀流试验测得的锁定发生的关键区域。由图5(b)可知，小KC数下，激励力系数变化较大。除在峰值区，KC=31.5时激励力系数基本保持在较低水平。这可能是因为此工况下流速及流向变化较快，漩涡尚未脱落便遭遇反向流速，无法形成稳定的漩涡脱落周期。KC数增大后，激励力系数也逐渐增加，但随着KC数变大，流态逐渐趋于定常流，此时激励力系数值将保持相近。

图5 低振幅比(0.50D)流体力系数对比Fig. 5 Contrast of fluid force coefficients at low amplitude ratio (0.50D)

图6给出了振动幅值1.00D时，升力系数幅值和激励力系数随KC数的变化情况。随着振动幅值的增加，KC数对流体力系数的影响效应并不显著，这与小振幅工况时是不同的。图6(b)表明，激励力系数随着无因次频率的增加逐渐减小。大KC数工况相对较为稳定，小KC数工况KC=31.5时振荡流周期变化较快，漩涡脱落相对不稳定，规律性较弱，故在低无因次频率下呈现出了下降趋势。但整体来看，各工况下激励力系数均小于0，此时流体对结构振动起阻尼作用，能量由结构传向流体。激励力系数由升力系数幅值和相位角计算得到，由Gopalkrishnan[3]均匀流试验可知，低振幅比与高振幅比时相位角变化具有明显不同。振荡流工况下高低振幅比时相位角变化趋势也不同，1.00D时相位角均保持为负值，故激励力系数变化趋势与0.50D存在较大差别。

图6 高振幅比(1.00D)流体力系数对比Fig. 6 Contrast of fluid force coefficients at high amplitude ratio (1.00D)

综上可知，在小振动幅值时，当圆柱的无因次频率处于涡激振动锁定区间内，KC数对激励力系数影响较大，且大、小KC数下呈现不同的规律性。在大振动幅值时，KC数对该雷诺数下流体力系数的影响将逐渐减小，激励力系数保持负值。

3.3 漩涡形态分析

基于振幅比0.50D计算结果，分析KC数对漩涡发放形态的影响。前文分析可得，振动频率在0.15～0.20之间流体力系数受KC数变化较为敏感，因此选择频率0.17作为分析重点。将5组KC数分为两类，其中，小KC数工况为组1和组2(KC=31.5和KC=63.0)，大KC数工况为组3、组4和组5(KC=126.0、KC=252.0和KC=503.9)。由于大小两类工况下不同KC数的漩涡形态具有相似性，故大KC数和小KC数各选择一组典型形态进行分析。

3.3.1 大KC数工况

图7给出KC=503.9、无因次频率0.17、一个振荡流周期内漩涡脱落发展状态。随着流速的逐渐增大，尾涡在圆柱壁面上产生，并逐渐向圆柱右侧延伸。当到达图7(a)时，有漩涡从尾涡后方脱落。图7(b)对应流速达到峰值，在由(a)到(b)的过程中，速度进一步增加，圆柱后方的尾涡也逐渐伸长。此阶段内伸长的尾涡与之前脱落的漩涡连在一起，并未出现显著漩涡脱离现象。经过图7(b)后，流速逐渐下降，此时圆柱后方的漩涡开始上下交替脱落，呈现出典型的2S泄涡模式，如图7(c)所示。圆柱运动到图7(d)时，此时圆柱后方依旧有漩涡脱落，但已不再是2S模式，并且漩涡的大小和强度均开始下降。

图7 大KC数下漩涡脱落(A*=0.5，f*=0.17，KC=503.9)Fig. 7 Vortex shedding at large KC number (A*=0.5，f*=0.17，KC=503.9)

当流速继续下降时，脱落的漩涡也逐渐在圆柱后方消散。图7(e)显示此时流速很小，圆柱壁面已无法产生尾涡。之后流速反向并逐渐增加。漩涡开始在圆柱的左侧产生，但此时尾涡较短，漩涡在离圆柱较近的位置脱落，如图7(f)所示。流速进一步增加，尾涡逐渐延伸，但无法观察到脱落的漩涡。这一阶段类似图7(a)至7(b)。图7(g)时流速刚过最大值，漩涡开始从尾涡上脱落，并呈现出2S泄涡模式，与流速正向不同的是，这一阶段2S模式的持续时间要更小。图7(h)后，流速逐渐下降，圆柱的尾涡长度慢慢变短，强度逐渐减弱。

3.3.2 小KC数工况

图8为KC数为31.5时的尾涡演化。图8(a)时，流速刚经过最大值，尾涡在圆柱后侧产生延长。随着流速下降，从图8(b)可以发现有单个漩涡从圆柱壁面上脱落，但由于没有长时间的单侧流向，难以观察到漩涡一个接一个在圆柱的下流向消散的现象。图8(c)时，流速反向，脱落的漩涡被反向流速带回到圆柱附近，导致圆柱周围漩涡分布复杂化，这也是小KC数时的典型状态。之后流速增加，尾涡在圆柱的左侧产生。

图8 小KC数下漩涡脱落(A*=0.5，f*=0.17，KC=31.5)Fig. 8 Vortex shedding at low KC number (A*=0.5，f*=0.17，KC=31.5)

在f*=0.17下，振荡流随着KC数的增加，漩涡脱落的模式各有不同。小KC数下流体速度变化较快，漩涡从圆柱壁面脱落后难以保持稳定泄涡，当流速反向后，脱落的漩涡又被冲向圆柱，使得圆柱周围压力变化更加复杂。在大KC数时，漩涡以典型的2S模式从圆柱壁面脱落，且随着速度下降，漩涡强度逐渐下降，脱落的漩涡也可在流速反向前消散完成。

3.4 升力系数时历分析

为更好地分析漩涡产生发展过程和升力系数之间的关系，提取了A*=0.5，f*=0.17对应的升力系数时历曲线，并对照漩涡发放过程，探讨二者之间的内在联系。

图9展示了大KC数时的升力系数时历曲线。3种大KC数工况均观察到了振幅调制现象。图9(c)给出了KC=503.9时升力系数变化情况。该曲线可以看到两个升力系数增大区，即振荡流速度幅值附近(72 s和88 s)。当运动到64 s时，振荡流速度由0开始逐渐增大，升力系数也逐渐增加。当圆柱运动至72 s时，进入到2S发放模式，升力系数进入第一个峰值区。80 s时，流速降低为0，可以观察到此时升力系数要小于峰值区。流速反向后，速度开始增加，再次进入到2S发放模式，此时升力系数同样达到峰值区。