剪切流下海洋立管涡激振动的三维数值模拟

罗冬冬，朱仁庆

(江苏科技大学 舰舶与海洋工程学院，江苏 镇江212003)

0 引 言

近年来，由于地球人口急剧增加和能源的不断消耗，各国能源危机越发明显。同时伴随陆上石油资源的日益枯竭，使得世界油气的开发重点逐步向海洋转移，尤其是深海区域。海洋立管在深海海洋平台中应用最为广泛，是连接海底的资源与海上作业平台的关键设备，但也是薄弱易损的构件之一。一旦采油立管在涡激振动作用下发生断裂，带来的不仅仅是经济上的巨大损失，严重的是将会造成巨大的环境污染。总而言之，目前对海洋立管进行VIV 问题的研究已成为国内外的热点之一。

经过数十年来对涡激振动的研究，学者们虽然还未完全把握涡激振动的机理，但取得了许多阶段性的成果，构成了当今涡激振动研究的基石。无论是实验研究还是数值模拟，大多数的研究主要针对在均匀来流情况下的涡激振动，如张建侨［7］研究质量比对柔性细长立管涡激振动的影响。在同一流速下，质量比大的立管模型响应位移小，其中空管的涡激振动响应一直处于大振幅的锁定状态下。在相同流速下，质量比大的立管模型所激起的模态更高。王嘉松等［4］对1 根长502.92 m，直径为0.534 m 的细长立管的进行模拟，指出了当立管发生大的、非均匀弯曲变形时，多模态的振动和线性动力响应不应该被忽视。同时表明预张力的增加，能有效的减小响应振幅，从而有效地减小疲劳损伤。但真实海洋来流情况往往比较复杂，不可能出现这种均匀来流的情况，相反更类似于具有一定斜率的剪切流的情况。因此剪切流下的立管涡激振动的模拟更加贴近海洋中的真实情况，符合实际的工程需要。目前，对于剪切流下的立管涡激振动的研究还十分少，尤其在国内无论是实验研究还是数值模拟都十分稀少。在国外，2003 年ExxonMobil［1］公司在MARINTEK 的深水水池中利用旋转的架子进行了均匀流和剪切流的一系列的立管涡激振动系列实验，说明了只有在低阶模态参与的响应中，轴向疲劳损伤十分重要。近年来，Remi Bourguet［3］研究了立管在剪切流下的锁定情况。结果表明在所有模拟条件中，在涡的发放与垂向振动局部同步时，在高流速区域并且长度至少超过圆柱长度的30%，才会发生锁定现象。在锁定情况下，振幅与在相似雷偌数下受迫振动的刚性圆柱的局部频率有关。梁勇［12］基于Matlab 编写了预测涡激振动响应的计算程序，来模拟剪切流下的立管。结果表明约束条件和立管密度对立管的振动响应有着十分重要的影响。

本文利用Workbench 平台，建立三维流固耦合模型，自编了进口剪切流的程序，采用了k － ω 模型，结合基于薄壳理论的有限元方法，并通过一种新的方法System Coupling 实现流－固耦合交界面的数据交换，从而实现剪切流下三维流固耦合模拟。与相关文献进行对比，研究立管的振动变形响应及尾涡形态特征等问题。

1 数值模拟方法概述

1.1 两场控制方程

1.1.1 流体场控制方程

流场海水可看做粘性不可压缩的流体，流场满足以下微分方程:

1)连续方程

连续方程是质量守恒定律在流体力学中的表现形式，该方程可以表述为:单位时间内流入该微元的净质量等于单位时间流出该微元的净质量。对于三维流动，其平面直角坐标形式可表示如下:

式中u，v，w 分别为流体在x，y，z 方向上的速度，m/s。

2)动量方程

动量方程是动量守恒定律在流体力学中的表现形式，其本质是满足牛顿第二定律。该方程可以表述为:对于一给定的流体微元，其具有的动量对时间的导数等于作用在该微元上的各种力之和。其数学表达形式为Navier －Stokes 方程(简称N －S)方程，对于不可压缩流体，写成Cartesian 坐标形式为:

式中:f 为流体在x，y，z 各方向上所受质量力，m/s2;p 为流体压强，Pa;ρ 为密度，kg/m3;ν 为运动粘性系数，m2/s。

1.1.2 固体场控制方程

对海洋立管的动力响应采用基于三维实体单元的有限元方法进行模拟。对运动方程进行有限元离散后可得到有限元方程:

2 算例分析

2.1 模型的建立

本文模型参数如表1 所示。

表1 立管模型参数Tab.1 Parameters of riser model

进口速度为剪切流v=(0.03 × Z +0.2)m/s，Z为立管长度，立管轴线与流场速度入口边界距离为10 D，与出口边界距离为30 D，与两侧对称边界距离为10 D。考虑的现有的计算机能力，使用较精密的平面网格体系，如图1 所示。

图1 二维网格拓扑图Fig.1 Two dimensional grid

流场采用k － ω 模型，迭代方法选用具有高稳定性的SIMPLEC 算法。选取立管的圆柱面为流固耦合面，受到动网格的限制，时间的离散方式选用了一阶隐式(1nd－Order－Implicit)。时间步长为t=0.002 s，计算总时长为6 s。

图2 结构有限元网格图Fig.2 Finite element of structure

对于结构模型的网格划分如图2 所示，同样结构中，设定立管的圆柱面为流固耦合面。同时，为了保持迭代耦合的一致性，时间步长与时长同流场相一致。

2.2 海洋立管的运动分析

2.2.1 顺流向的运动分析

在顺流向方向，立管因为受到拖曳力的影响，发生弯曲变形，如图3 所示，在t = 2.57 s 时，弯曲达到最大，立管的最大弯曲变形处发生在约z/L=0.6 处，最大变形为1.2 D，此后，流场趋于稳定，立管变形很小，只是在最大位置处产生小幅振动，振动幅值为0.1 D。如图4 所示。与文献［11］对比发现:其运动趋势与均匀流下立管的运动趋势相同，并都多出现了顺流向的不可忽略的振动，只有最大弯曲变形位置向上移动。该结果说明由于剪切流流速递增的特性，立管受到水流的拖曳力的中心向上移动，使得立管的最大弯曲变形的位置上移，但顺流向上对立管的运动趋势的影响与均匀流下相同。同时说明了System Coupling 实现流固耦合具有一定的可行性。

图3 顺流向运动效果图Fig.3 In-line motion of riser

图4 最大位移处点的位移随时间的变化图Fig.4 Scatter diagram of displacement varying with time of point with largest displacement

2.2.2 横流向的运动分析

从图5 可看出，在横流向上，立管的振动在1阶、2 阶振动模态之间往复，立管的弯曲变形呈非对称性，当流体趋于稳定时，最终锁定在2 阶振动模态，同时出现了3 阶振动模态的情况。

而根据文献［11］得到立管自振频率公式为:

式中:md为单位长度隔水管排开海水的质量;CA为附加质量系数，一般取1.0;I 为截面惯性矩;T 为预紧力。同时根据“锁定”的定义:

图5 立管不同时刻横流向的振动效果图Fig.5 Cross-flow vibration shape of riser at different moment

式中:f 为隔水管实际振动频率;fn为隔水管理论自振频率;fst为涡脱落频率。

由此得到前3 阶振动模态的激发流速:1 阶模态的激发流速0.188 m/s;2 阶和3 阶模态的激发流速为0.389 m/s 和0.607 m/s。而本文的最高流速为0.488 9 m/s，最低流速为0.2 m/s，并没有达到3 阶模态的激发流速。本文激发2 阶振动模态的来流约占立管长度的40%，符合上文Remi Bourguet 的指出:在涡的发放与垂向振动局部同步时，在高流速区域并且长度至少超过圆柱长度的30%，锁定才会发生。所以立管的振动最终锁定在2 阶振动模态。同时激发1阶振动模态的来流约占立管长度的60%，对立管的影响也很大，因此可能导致了3 阶振动模态的出现。该现象表明:剪切流增加了立管振动的复杂性，立管出现了比剪切流最大流速下的振动模态更高1 阶的3 阶模态振动，但最终锁定模态处于2 阶振动模态。

2.3 海洋立管的流场分析

图6 为立管在前期和后期不同截面的尾流区涡的三维特性图。图中5 个截面分别对应Z/L=0.1，0.3，0.5，0.7，0.9 的涡量图。在t = 1.2 时，由于剪切流的关系，在Z/L=0.7 和0.9 的截面上都形成了涡的发放，而其他3 个界面涡还没有形成。当最终整个立管的VIV 稳定后，涡的发放形式主要是“2S”形的涡。从图中可以看出:由于立管的运动和剪切流的作用，涡在立管不同的位置有明显的区别，在Z/L=0.1，0.5，0.9 处，立管在横流向的运动范围比较小;两涡在横流向上的距离明显小于Z/L=0.3，0.7 处。

在立管涡的发放稳定前，还伴随其他不同形式的涡，如图7 所示。图7 为t = 2.4 s 时，Z/L=0.3、0.5 的截面的涡量图。图7(a)为Z/L=0.3 截面涡量图，此时涡为“P +S”的形式。而到图7(b)时即Z/L =0.5 截面，涡为“2P”的形式。

图6 立管不同截面不同时刻的涡量图Fig.6 Vorticity of riser at different sections for different moment

图7 t=2.4 s 时涡量图Fig.7 Vorticity of riser at t=2.4 s

3 结 语

本文基于计算流体力学(CFD)和计算结构力学(CSD)方法，借助商业软件Ansys 14.0，参考长径比为482 的均匀流的立管模拟，实现了VIV 的三维数值模拟。结果表明，数值模拟与均匀流有部分相似;大长径比柔性隔水管发生VIV 时流场受到明显扰动，涡的脱落呈现如2P、2S、P+S 等模态;立管最大位移处向上移动，立管在横流向上出现了大的弯曲变形，出现了1 阶、2 阶模态之间的转变，最终锁定在2 阶模态。

［1］TOGNARALLI M A．VIV response of a long fexible cylinder in uniform and linearly sheared currents［C］//Offshore Technology Conference，2004．

［2］FARSHIDIANFAR A，ZANGANEH H． A modified wake oscillator model for vortex-induced vibration of circular cylinders for a wide range of mass － damping ratio［J］．Journal of Fluids and Structures，2010．

［3］BAHMANI M H，AKBARI M H． Effects of mass and damping ratios on VIV of a circular cylinder［J］． Ocean Engineering，2010．

［4］REMI B，MICHAEL S T． Lock-in of the vortex-induced vibrations of a long tensioned beam in shear flow［J］．Journal of Fluids and Structures，2011．

［5］WANG Jia-song，ZHAN Lu-lu． Three dimensional numerical simulation of vortex induced vibration for a 500 － m － long marine riser［C］//ISOPE，2012．

［6］LIN Li-ming，ZHONG Xing-fu，WU Ying-xiang． Vortexinduced vibrations of a circular cylinder with different geometric disturbances［C］//ISOPE，2012．

［7］HUANG Zhi-yong，CUI Wei-cheng．Vortex-induced vibration of flexible riser in time domain［D］． Shanghai:Shanghai Jiao Tong University，2008．

［8］ZHANG Jian-qiao，SONG Ji-ning，LV Lin，et al．Analysis on the vortexoinduced vibrations experimental results of flexible risers［J］． China Offshore Platform，2009．

［9］TANG Guo-qiang，TENG Bin，LV Lin，et al．Development and application of an empirical model for vortex-induced vibration of flexible deep water riser［J］．China Offshore Platform，2010．

［10］CAI Jie． Investigation on numerical simulation of the fluidsolidInteraction characteristics for oceanic riser［D］．Shanghai:Shanghai Jiao Tong University，2010．

［11］CAI Jie，YOU Yun-xiang，LI Wei，et al． The VIV characteristics of deep-sea risers with high aspect ratio in a uniform current profile［J］．Chinese Journal of Hydrodynamics，2010．

［12］WANG Cheng-guan，TIAN Zhong-xu，QIAO Xin-qi． Study on three dimensional numerical simulation of VIV characteristics for marine riser［J］． Chinese Journal of Hydrodynamics，2011．

［13］LIANG Yong．Response analysis of flexible riserssubjeted to Vortex-induced vibration in sheared flow［D］． Harbin:Harbin Institute of Technology，2011．