仙人掌形截面柔性圆柱体涡激振动抑制研究*

娄 敏 陈法博 时 晨

(中国石油大学(华东)石油工程学院)

0 引 言

立管是深海油气开发系统中连接海面作业平台和海底钻采设施的关键结构。当海流经过立管时，在立管的两侧产生交替脱落的漩涡，从而在结构表面形成周期性的脉动力引起立管振动，当漩涡脱落频率接近立管的某阶自振频率时，漩涡脱落频率和立管振动频率相互锁定导致立管大幅振动，这种振动称为涡激振动。涡激振动是导致海洋立管及其他圆柱形结构物疲劳破坏的主要原因之一。实际工程中为减小或者消除涡激振动的影响，延长海洋立管的使用寿命，往往需要采用涡激振动抑制措施。

目前关于抑制涡激振动的研究较多，宋振华等[1]针对3根附属杆的涡激振动抑制效果进行分析研究；娄敏等[2]对减振器、分隔板和控制杆的抑制效果进行了对比；刘志慧等[3]通过波流水槽试验对盘球的抑制效果进行了研究。工程上常用的涡激振动抑制装置有螺旋列板和整流罩等。螺旋列板会增大立管受到的拖曳力；整流罩体积大、不易储存，而且如果因为水生生物附着等导致整流罩不能自由转动，整流罩不但无法抑制涡激振动，而且还会引起很大的阻力。鉴于螺旋列板和整流罩等传统抑振装置的不足，本文根据仿生学原理，设计了仙人掌形截面抑振装置。

当前关于仙人掌形截面研究主要集中在试验和静止圆柱数值模拟方面，S.TALLEY等[4]利用风洞试验对不同凹槽深度的仙人掌形截面圆柱和圆形截面圆柱在雷诺数(2～20)×104的流场中进行了研究，研究结果显示，随着凹槽深度的增加，仙人掌形截面的拖曳力系数逐渐减小。Y.YAMAGISHI等[5]研究了高雷诺数下具有沟槽的圆柱体涡激振动响应特性，结果表明带有沟槽的圆柱体的拖曳力比相同雷诺数下的光滑圆柱体减小了约15%。WANG J.等[6]利用ANSYS Fluent模拟了低雷诺数情况下固定仙人掌形截面圆柱绕流，研究结果表明，仙人掌形截面可以降低拖曳力、升力和斯特劳哈尔数。余俊等[7-8]利用激光例子图像速度场PIV(Particle Image Velocimetry)测试技术对仙人掌形截面圆柱进行了试验研究，分析了仙人掌形截面圆柱绕流流场的特性。钱权等[9-10]对仙人掌形截面和圆形截面圆柱绕流进行了三维流场数值模拟，研究表明仙人掌形截面具有较好的减阻效果。本研究通过ANSYS Fluent动网格技术并结合波流水槽试验对相同脊高、不同脊数的仙人掌形立管涡激振动响应进行了研究，同时分析了不同脊数的仙人掌形截面的抑振效果。

1 数值方法

1.1 动网格模型

数值模拟圆柱涡激振动的过程简述如下：首先利用ANSYS Fluent对流场N-S方程进行求解，然后通过自定义函数(UDF)提取作用在圆柱上的流体力，并将力代入圆柱运动的控制方程中，用四阶龙格库塔法求解方程得到圆柱体运动的速度和位移，再通过UDF中的宏函数实现圆柱的运动和网格的更新。Fluent中动网格有3种更新方法：弹簧光顺法、动态层法和局部重构法。本文采用弹簧光顺法和局部重构法相结合来实现网格变形及重构。

弹簧光顺模型将网格的边看作连接节点的弹簧，初始网格间距相当于弹簧系统的平衡状态。网格边界节点发生相对位移后，产生与位移成正比例的力，力的大小由胡克定律计算，在外力作用下，弹簧系统经过调整将达到新的平衡状态。

当运动边界位移超过网格尺寸太多时，应用弹簧光顺模型会导致网格质量降低，甚至会出现网格线交叉重叠，从而导致计算错误，或者引起网格畸变导致迭代不收敛。因此，本文采用弹簧光顺法和局部重构法相结合的模型，该模型可以将弹簧光顺模型生成的畸变过大或尺寸变化过大的网格进行局部重新划分。网格尺寸和网格畸变率是判断网格是否需要重新划分的两个指标，在计算过程中如果网格尺寸不满足以下任一条件，网格就会标记并重新划分：①网格尺寸小于设定的最小网格尺寸；②网格尺寸大于设定的最大网格尺寸；③网格畸变率大于设定的畸变率。

1.2 涡激振动微分方程

本研究将三维立管模型简化为二维弹簧-质量-阻尼系统，建立的圆柱振动模型如图1所示。通过ANSYS Fluent研究其涡激振动响应，振动控制方程为：

(1)

(2)

图1 圆柱振动模型Fig.1 Cylindrical vibration model

2 数值模型

2.1 模型参数

根据T.K.PRASANTH等[11]对网格尺寸的研究，当满足尾流区长度大于25D、整体高度区域大于20D时(其中D为圆柱的外径)，圆柱体双自由度涡激振动的响应将不再受流体区域边界的影响。

流场区域计算模型如图2所示。选取圆柱直径D=0.017 m，模拟的均匀流场为28D×44D的矩形，立管的尾流区长度为30D，整个流场划分为“刚性运动区域(区域③)+动网格区域(区域②)+静止网格区域(区域①)”3个区域， 壁面设置15层四边形结构化网格，向外过渡采用非结构化网格直至与动网格区域交界处。为了保证圆形截面附近的网格具有初始的高质量状态，将刚性运动区域设置成直径为3D的圆形区域；动网格区域设置成直径为16.5D的圆形区域，采用Smoothing + Remeshing 的方式更新网格，此区域为网格变形区域，采用非结构化网格；静止区域采用结构化网格。网格划分结果如图3所示。

图2 流场区域计算模型Fig.2 Calculation model of flow field area

图3 网格划分Fig.3 Mesh generation

数值模拟参数取值如下：系统质量为0.338 kg，系统阻尼比为0.117 6，系统刚度为96.221 6 N/m，固有频率为2.686 0 Hz。圆形截面和仙人掌形截面圆柱所取参数相同。圆柱的直径D为0.017 m，仙人掌形截面最大外径为0.018 4 m，脊高度为0.04D。R8-4表示脊数为8、脊高为0.04D的仙人掌形截面圆柱，如图4所示；R20-4表示脊数为20、脊高为0.04D的仙人掌形截面圆柱；R24-4表示脊数为24、脊高为0.04D的仙人掌形截面圆柱。求解设置时，采用SSTk-ω湍流模型，对流项采用二阶迎风离散格式，时间项采用全隐式格式；控制方程中速度和压力的耦合采用PISO算法。左侧边界采用速度入口(inlet)；右侧边界采用压力出口条件(pressure outlet)；上下边界采用对称边界条件(symmetry)；立管表面采用无滑移壁面条件(no-slip wall)。

2.2 模型验证

验证上述计算模型与网格划分的可靠性，在雷诺数Re=3 900时进行单圆柱绕流数值模拟，并将结果与其他文献中的升力、阻力系数和斯特劳哈尔数对比。利用MATLAB进行频域分析，使用其内置函数“pwelch”计算升力系数的功率谱密度(PSD)得到漩涡脱落频率fs，再通过fs计算出斯特劳哈尔数St。其表达式为：

(3)

式中：u为来流速度，m/s。

图4 脊数为8、脊高为0.04D的仙人掌形截面图Fig.4 A cactus-shaped cross-section with 8 ridges and a ridge height of 0.04D

当流体流经圆柱时，圆柱两侧产生漩涡泄放，在横流向会产生升力Fl使得圆柱发生横向运动，并在顺流向产生拖曳力Fd使得圆柱能够在顺流向运动，升力和拖曳力分别表示为：

(4)

(5)

式中：Cl为升力系数；Cd为拖曳力系数；ρ为流体密度，kg/m3。

随着流场域网格数量的增加，拖曳力系数和升力系数逐渐趋于稳定，考虑到计算时间随着网格数量的增加而延长，在保证计算结果准确的前提下，本研究将流场域的网格数量取35万。

表1 雷诺数3 900时流场网格数敏感性分析Table 1 Sensitivity analysis of flow field grid number under Reynolds number 3 900

表2 雷诺数3 900时单圆柱绕流模拟结果与文献对比Table 2 Comparison of the simulation results of the flow around a single cylinder under the Reynolds number 3900 with the results in literatures

2.3 计算结果

多组离散数据y1，y2，……，yn的均方值可以表示为：

(6)

均方值开平方后取正值可得均方根值，该值能够反映出圆柱在每个时刻偏离平衡位置的程度，表示为：

(7)

4种不同截面单圆柱无量纲位移均方根值对比情况如图5所示。

由图5可以看出：在约化速度4～8范围内，R20-4和R24-4仙人掌形截面横流向无量纲位移均方根值均大于圆形截面，说明R20-4和R24-4没有抑制涡激振动的效果，且增大了振动幅值；R8-4仙人掌形截面横流向无量纲位移均方根值小于圆形截面，说明R8-4仙人掌形截面具有较好抑制涡激振动的效果。R8-4仙人掌形截面顺流向无量纲位移均方根值小于圆形截面，说明R8-4不仅可以抑制涡激振动，还可以降低立管在顺流向的振动幅值。

图5 4种不同截面单圆柱无量纲位移均方根值对比Fig.5 Comparison of the dimensionless displacement root-mean-square value of a single cylinder with different cross-sections

3 试验模型

3.1 模型建立

2019年12月，在中国石油大学(华东)海工实验室的波流水槽完成试验。波流水槽的横截面宽0.8 m，高1.5 m，造流机可造稳定流速0.05～0.40 m/s。立管模型如图6a所示。立管模型采用有机玻璃管，将其部分浸入水中，管子下端封闭且管内充满水，两端采用万向节铰接固定在铝合金框架上， 铝合金框架的刚度足够大，其固有频率远高于立管模型。试验水深0.7 m，试验流速0.20～0.36 m/s，其对应的雷诺数范围3 400～6 120。立管模型参数如下：立管总长度L=1.5 m，浸入水中部分L1=0.7 m，外径d=0.014 m，壁厚t=0.002 m，弹性模量E=3.2 GPa，单位长度质量mp=0.089 kg/m，单位长度排开水质量md=0.216 kg/m；管套的最大外径D=0.017 m，单位长度质量ms=0.069 kg/m；立管等效质量密度mequ=0.291 kg/m。

在立管中部布置了4个应变传感器用于测量立管的弯曲微应变，传感器的布置型式如图6b所示。通过3D打印技术制作了不同截面型式的仙人掌形管套和圆形管套，保证不同截面形状管套的横截面积相同，然后将管套布置在立管浸入水中的部分。为了减小管套对立管刚度的影响，将管套分成每节0.017 m的小段，且相邻的管套之间留有0.003 m的空隙。

图6 立管模型和应变传感器布置示意图Fig.6 Schematic diagram of standpipe model and strain sensor layout

3.2 单独管水中自振试验

3.2.1 单独管水中自振频率测量

对装有圆形截面管套的立管模型进行水中自振试验，振动时程曲线和频谱图如图7所示。

图7 圆形截面立管水中自由振动时程曲线和频谱图Fig.7 Free vibration time history curve and frequency spectrum of a circular cross-section riser in water

圆形截面立管模型静水中的自振频率为2.686 Hz。同样，测得不同脊数的仙人掌形截面静水中自振频率为2.686 Hz。

3.2.2 立管水中自振频率理论计算

对于两端铰支固定的均匀梁，可由Euler Bernoulli梁理论计算出单管立管模型在水中的自振频率，计算公式如下:

(8)

式中：m为立管模型单位长度质量，kg/m；I为截面惯性矩，m4。

立管模型只有部分浸入水中且水下部分套有管套，管内充满水，则单位长度质量m表示为：

(9)

式中：z为距离立管底端的高度，m；Ca为附加质量系数，圆柱体一般取1。

S.LENCI等[16]提出了一种计算非均匀梁自振频率的方法，采用等效单位长度质量mequ代替公式(8)中的m，其表示为:

(10)

通过理论计算得到立管模型水中自振频率理论值fw=2.740 Hz。水槽试验测得立管模型静水中自振频率为2.686 Hz，理论值和试验值较为接近。

3.3 单独管波流水槽试验

本研究进行了3种截面单独管的波流水槽试验，3种截面分别为：圆形截面立管、R8-4型截面立管和R20-4型截面立管。计算了立管在不同流速下的峰值频率、微应变均方根值以及年疲劳损伤率。

3.3.1 峰值频率对比

利用MATLAB进行频率分析，使用其内置函数“pwelch”计算立管微应变的功率谱密度，得到峰值频率，如图8所示。

从图8可以看出，3种立管模型的峰值频率都随着约化速度的增大而增大，且都存在“倍频”现象。R8-4仙人掌形截面立管的频率比其他2种立管要大，可能是其尾流区宽度变窄造成的。

3.3.2 微应变均方根值对比

3种不同截面单独管微应变均方根值对比情况如图9所示。从图9可以看出，在约化速度4～8范围内，圆形截面立管和R20-4型立管微应变均方根值随着约化速度增大而增大，而R8-4型立管呈现先减小后增大趋势，每个流速下横流向微应变的均方根值均低于圆形立管，且其横流向微应变的均方根值最大可降低44%。

图8 3种不同截面单独管的峰值频率对比Fig.8 Comparison of peak frequencies of three individual tubes with different cross-sections

图9 3种不同截面单独管微应变均方根值对比Fig.9 Comparison of micro-strain root-mean-square values of three individual tubes with different cross-sections

3.3.3 年疲劳损伤率对比

立管的年疲劳损伤率可使用雨流循环计算法和Miner准则计算。本文仅出于说明目的，借用钢材的S-N曲线[17]来计算有机玻璃管的疲劳损伤，其表达式如下：

N=aS-b

(11)

式中：N为应力范围内达到失效时的周期数，S为应力范围，a=4.266×1011，b=3.0。

从图10可以看出，R8-4型立管横流向和顺流向的疲劳损伤率相对于圆形截面立管和R20-4型立管大幅度降低，说明脊数为8、脊高0.04D的仙人掌形截面管套具有良好的振动抑制效果。

图10 3种不同截面单独管疲劳损伤率对比Fig.10 Comparison of fatigue damage rate of three individual tubes with different cross-sections

4 数值模拟和水槽试验对比

本研究选择圆形截面和1组仙人掌形截面的数值模拟数据与水槽试验进行对比，利用张建侨等[18]模态分解的方法将试验采集的微应变值转化为该截面的位移值。将数值模拟得到的无量纲位移均方根值与水槽试验值进行对比，结果如图11和图12所示。由图11和图12可以看出，数值模拟值和水槽试验值具有相同的趋势，各个流速下的无量纲位移均方根值误差均在15%以内，数值模拟和水槽试验符合性较好。由于数值模拟没有考虑立管三维效应的影响，使得数值模拟和波流水槽试验结果存在一定的误差。

图11 圆形截面数值模拟结果和水槽试验结果对比Fig.11 Comparison of numerical simulation results of circular cross-section and channel experiment results

图12 R8-4仙人掌形截面数值模拟结果和水槽试验结果对比Fig.12 Comparison of numerical simulation results of R8-4 cactus-shaped cross-section and channel experiment results

5 结 论

(1)数值模拟结果显示：在约化速度4～8范围内，R8-4型仙人掌形截面不仅能够降低立管横流向的振动幅值，具有抑制涡激振动的效果，还能大幅降低立管在顺流向的振动幅值；而R20-4和R24-4型仙人掌形截面没有抑制涡激振动的效果。

(2)水槽试验结果显示：在约化速度4～8范围内，R8-4型仙人掌形截面立管各流速下的频率值大于圆形立管，其横流向微应变均方根值和年疲劳损伤率均小于圆形立管，且横流向微应变均方根值相比于圆形立管最大可降低44%；R20-4型仙人掌形截面立管在各流速下的频率值小于圆形立管，但其横流向微应变均方根值和年疲劳损伤均大于圆形立管，说明R8-4型仙人掌形截面具有较好的振动抑制效果，而R20-4型仙人掌形截面不具有振动抑制效果。

(3)数值模拟和波流水槽试验结果符合性较好，说明数值模拟具有一定的可靠性，当水槽试验不能满足研究要求时，可以通过数值模拟来实现。