海洋热塑性增强管(RTP)涡激振动数值计算1)

2020-02-23陈东阳王国平

力学学报 2020年1期

芮雪陈东阳王国平

∗(南京理工大学发射动力学研究所,南京 210094)

†(扬州大学电气与能源动力工程学院,江苏扬州 225127)

引言

当海流绕过海洋立管时,在立管下游会产生尾流和漩涡,周期发放的漩涡对立管产生垂直于来流向的涡激升力,引起立管的涡激振动(vortex-induced vibration,VIV).立管涡激振动即立管尾流涡脱落的频率与结构频率接近时发生的自激振动,一直以来是立管结构疲劳破坏的一个主要因素[1-4],由于海流速度沿水深变化的非均匀性及流固耦合的复杂性,涡激振动长期以来是海洋工程中颇具挑战性的问题之一,涡激振动的准确预报是一个巨大的难题[5].

海洋热塑性增强管(reinforced thermoplastic pipe,RTP)是一种复合材料管道,可用作立管,具有可盘卷(近百米一卷)、海陆施工快捷、质量小、耐疲劳、保温性能好,以及高耐腐蚀等性能[6-10].使用RTP 作为立管可以减小海洋平台的体积,减少海洋平台锚张力、减少立管顶张力,显著降低开采成本.同时,使用复合材料立管操作成本低,对系统其他构件的影响相对较小,相比传统的钢制立管有更多的优越性[11-12].RTP 在海上优势突出,虽然要求高,投资大,但发展迅猛,海洋用RTP 市场可能超过陆用市场.我国海洋油气资源丰富,过去油气开采限制于100 m 深以内的浅海,近年向更深海域拓展,满足深海要求的RTP 成关键.国内对复合材料立管的研究起步较晚,关于完全非金属的复合材料RTP 立管涡激振动的研究更少.能够适用于工程上快速评估复合材料立管涡激振动性能的方法未见报道.RTP 产品的种类繁多,本文仅以如图1 所示的复合材料预制增强带RTP 立管作为研究对象,探索适合工程的复合材料立管的建模方法,以及研究其涡激振动行为,为工程上类似的复合材料立管快速建模和涡激振动快速分析提供参考.

1 海洋立管系统动力学模型

图1 RTP 结构示意图Fig.1 Structural diagram of RTP

如果水深小于500 m,一般采用传统的固定式平台或塔式平台.假设本文研究的立管长度在50∼150 m 之间,因此采用塔式平台,即固定式平台.针对密度与水接近的RTP 立管,忽略立管重力,考虑均匀来流和剪切流引起的涡激升力、立管顶部的顶张力、连接管道的刚性接头,简化后的RTP 立管模型示意如图2 所示,立管下端铰接在位于海底的万向节上,上端铰接于固定平台上.图中T为顶张力,数字1 表示均匀来流,数字2 表示剪切流.本文研究所用到的立管参数如表1 所示.忽略平台的振动,立管边界条件可视为两端简支,另外做如下假设:

(1)不考虑立管的扭转转动;

(2)立管内部流体为油,视为附加质量,不考虑内流速度对立管振动的影响;

(3)忽略立管本身的结构阻尼,仅考虑立管振动引起的水动力阻尼;

(4)立管应变很小,应力应变为线性关系;

(5)海流沿水深为均匀来流或者剪切流.

RTP 管道一般由内管、增强层和外管组成.如表1 所示,本文研究的RTP 立管的材料包括高密度聚乙烯(HDPE)、中密度聚乙烯(MDPE)和芳纶纤维.立管参数见表2.

图2 立管模型示意图Fig.2 Riser model

表1 RTP 管道材料性质Table 1 RTP material properties

表2 立管参数Table 2 Riser parameters

考虑柔性立管刚性接头等结构细节,基于多体系统传递矩阵法计算其振动特性.铰−铰约束的带有顶张力的立管的欧拉梁模型如图3 所示.本文采取将梁分段的方法来建立柔性立管的力学模型.为了便于考虑立管中刚性接头对振动特性的影响,以刚性接头的长度为参考长度,将RTP 立管分为n段,那么刚性接头的长度为L/n(本节针对150 m 长RTP 立管计算时,取n=80).这样计算中可以任意选取刚性接头的个数,只需将对应的梁元件的材料改为刚材料即可.

图3 包括状态矢量和传递方向的立管模型Fig.3 Riser model including state vector and transfer direction

2 RTP 立管动力学特性建模

文献[13-15]基于多体系统传递矩阵法[16-18],均推导了带有预应力的欧拉−伯努利梁的传递矩阵.本文的立管模型也采用该梁模型,但立管所受的定张力方向向上,与文献[13-15]所推导的传递矩阵存在正负号差别.限于篇幅,此处不作重新推导.所得的该可以考虑横向和轴向振动的预应力欧拉−伯努利梁的传递矩阵为

式中

由于柔性立管非常细长,将立管分成四部分,并使用图3 所示的传递方向,以减小同一方向传递过程中由于过多矩阵相乘引起的高阶振型的数值误差.柔性立管的传递方程为

式(3)中,状态矢量的下标up 和down 分别表示联接点处的上方和下方对应的状态矢量.另外,由图3 中的坐标系和传递方向,以及MSTMM 书中的符号约定可知,式(3)中[16]

所以柔性立管的总传递方程

式中,Upart1=Ui···U3U2U1,Upart2=Ui+1···Uj−1Uj,Upart3=Uk···Uj+2Uj+1,Upart4=Uk+1···UnUn+1,n.U1,U2···Un分别代表柔性立管每一段传递矩阵,可以充分考虑柔性立管具体的结构细节.对于铰−铰连接的柔性立管的边界条件为

将系统边界条件代入式(5),可得系统特征方程

求解方程(9)即可得到系统的固有频率ωk(k=1,2,···),然后求出系统边界点状态矢量和Zall,进而通过元件传递方程得到对应于固有频率ωk的系统全部联接点的状态矢量,得到RTP 立管系统的振型.

3 RTP 立管涡激振动动力学方程

基于MSTMM 和Van der Pol 尾流振子模型建立RTP 立管涡激振动动力学方程为

其中,M和K分别为质量矩阵和刚度矩阵.令V=[V1V2···Vn],q=[q1q2···qn]T,n表示立管系统模态叠加所用的模态阶数,本节对于立管的涡激振动计算取前10 阶模态,即n=10.v=Vq,可以将式(10)转化到模态坐标系中,则有

根据文献[1,19]可得

利用增广特征矢量的正交性,在式(12)两边同时乘以VT,即

将式(15)、式(16)写成状态空间的形式,令x1=q,x2=x3=qv,x4=,则

采用龙格库塔法求解式(17),其初始条件q==01×10,qv=0.21×80,=01×80,时间步长大小一般取参考长度与来流速度比值的0.1∼0.01 倍,此处取0.01 s.

4 模型验证及参数获取

4.1 振动特性及涡激振动模型验证

以文献[20]中的立管参数为例,采用MSTMM计算了细长立管的振动特性.图4 为采用MSTMM计算顶张力为817 N 的立管的干模态(不考虑周围流体的影响下的模态,这种模态可以称为“干模态”,即不受流体影响的模态)前8 阶圆频率计算结果,与文献[20]中的理论分析解对比结果如表3 所示,误差小于0.12%,因此验证了基于MSTMM 立管动力学特性计算模型的准确性.

图4 干模态前8 阶圆频率Fig.4 First eight circular frequencies of dry modal

在海流作用下,立管上各点由漩涡脱落引起的振动可视为立管各模态振动的叠加(一般为前几阶湿模态,湿模态即考虑周围流体影响下的模态,工程上湿模态计算一般将立管外径对应的圆柱排开的流体质量作为附加质量考虑),且一般存在主导振动模态.以文献[20]中的立管为例,基于MSTMM 和Van der Pol 尾流振子模型,通过式(17)采用模态叠加法计算带有顶张力的立管的涡激振动响应.

表3 立管干模态前8 阶频率对比Table 3 Comparison of the first eight frequencies of riser dry modal

文献[20]中,均匀来流速度为0.2 m/s,直径为0.02 m,所以Re=3988,根据S t和Re的关系图[3],S t取0.23.Van der Pol 尾流振子模型的计算精度一定程度上取决于经验系数的选取,本文选取A=12 和ε=0.3,CL0=0.25,CD=1.2[3].如图5 所示,仿真得到沿立管轴向的涡激振动最大振幅分布.通过与文献[20]的实验数据及CFD 仿真数据对比,误差较小,验证了本文方法的可行性.

图5 沿立管轴向的涡激振动振幅分布Fig.5 Amplitude of vortex induced vibration along the axis of riser

4.2 RTP 立管等效刚度计算

顶张力作用下的简支梁固有频率的理论计算公式为

式中,T为顶张力,M为单位长度的RTP 管质量,L为立管长度,E为立管的弹性模量,I为立管的截面惯性矩.

在作者前期工作中[6],采用ANSYS 软件的shell 181 单元对RTP 立管进行了复合材料建模,并验证了建模方法的准确性.本文为了获取RTP 立管的等效刚度参数,采用ANSYS 软件计算一小段两端铰接的RTP 短立管(此处取10 m 长,直径为365 mm 的短管进行有限元复合材料铺层并进行模态计算,顶张力为1 N)的第一阶频率,即可以利用式(18)反推出立管的等效刚度EI.然后采用MSTMM 就可以方便的计算出不同长度、不同顶张力的RTP 立管的振动特性.

如图6 所示,为短立管的第一阶模态分析的结果,计算出短立管的第一阶频率为7.357 9 Hz,立管的单位质量为40.081 kg/m.令式(18)中n=1,因此可以反推出

为了验证所获得的立管弯曲刚度参数的准确性,采用MSTMM 计算150 m 长,直径D=365 mm 的复合材料RTP 立管的前10 阶频率和振型,并与ANSYS有限元软件计算结果对比分析.图7 为采用MSTMM计算得到的顶张力T=60 000 N 情况下的RTP 立管前10 阶圆频率.表4 为MSTMM 与ANSYS 有限元软件计算结果的对比,两种方法计算结果相差较小,说明采用短立管获取细长复合材料立管等效刚度的方法可行.

图6 短立管模态分析Fig.6 Modal analysis of short riser

图7 基于MSTMM 的RTP 立管前10 阶圆频率Fig.7 The first ten circular frequencies of RTP riser based on MSTMM

表4 MSTMM 与ANSYS 计算结果对比Table 4 Comparison of results between MSTMM and ANSYS

将MSTMM 计算出的前10 阶振型进行质量归一化,然后与ANSYS 计算出的振型进行对比.如图8所示,可以看出MSTMM 计算出的振型与ANSYS 软件计算出的复合材料立管的振型一致.

5 仿真结果

5.1 RTP 立管涡激振动特性分析

采用本文的方法省去了每改变一次参数就要对立管进行重新建模、复合材料铺层等繁琐步骤,也无需划分大量网格单元,计算量小,且矩阵阶次低,计算效率高.直接采用MSTMM 可以快速计算出RTP立管不同参数情况下的湿模态频率.如图9(a)所示,计算了顶张力为60 000 N 情况下,不同长度RTP 立管的湿模态频率.随着立管长度的增大,立管的高阶湿模态频率增长越缓慢,且立管长度越大,频率越小.说明立管长度越大,立管整体刚度越小.立管的顶张力大小一般取立管重量的1∼2 倍之间.150 m 长立管的重量为60 000 N.如图9(b)所示,分别计算了150 m长立管顶张力为60 000,80 000,100 000 N 情况下的湿模态频率.从图中可以看出,顶张力对立管的湿模态频率影响不大,随着顶张力的增大,RTP 立管湿模态频率微弱增大.

图8 RTP 立管的1,3,5,7,9,10 阶振型Fig.8 Modes 1,3,5,7,9 and 10 of RTP riser

图9 RTP 立管湿模态频率Fig.9 RTP riser wet modal frequency

5.2 RTP 立管涡激振动特性分析

一般表层海流的水平流速从0.01 m/s 到0.3 m/s,深处的水平流速则在0.01 m/s 以下.基于MSTMM 和Van der Pol 尾流振子模型,利用式(11)∼式(17)求解立管的涡激振动响应,取Van der Pol 尾流振子模型中系数A=12 和ε=0.3,CL0=0.25,CD=1.2[3],S t取0.2.

立管参数确定以后,立管的固有频率就可确定,但涡激升力的频率会随流速的增大而增大,这就导致不同流速下立管涡激振动激发的模态不同.从图10可以看出,150 m 长,顶张力为60 000 N 的立管,流速越大,立管涡激振动激发模态越高,但沿着立管轴向的最大振幅略有减小,说明立管涡激振动位移主要受低阶模态控制.这是由于150 m 长立管相对于50 m长立管刚度大幅度减小,容易激发出高阶模态为主导的振动响应.图11(a)、图11(b)、图11(c)分别是L=150 m,T=60 000 N 的立管在U=0.1,0.2,0.3 m/s时的y/D历时云图.图中水平坐标表示时间,垂向坐标表示立管长度方向位置,颜色代表的是不同的y/D的值.从图11(a)可以看出,立管主要是以第1 阶模态为主导的振动,图11(b)可以看出立管在来流速度为0.2 m/s 时,发生第1 阶和第2 阶模态切换现象.这是由于柔性立管各阶振型对应的频率比较相近造成的.图11(c)为立管在来流速度为0.3 m/s 时,发生以第3阶模态为主导的振动历时云图.

图10 沿RTP 立管轴向的涡激振动振幅分布Fig.10 Amplitude of vortex induced vibration along the axis of RTP riser

图11 沿RTP 立管轴向的涡激振动响应历时云图Fig.11 Contour of vortex induced vibration response along the axis of RTP riser

图11 沿RTP 立管轴向的涡激振动响应历时云图(续)Fig.11 Contour of vortex induced vibration response along the axis of RTP riser(continued)

同样的,计算立管不同长度,不同顶张力情况下的涡激振动响应,并绘出沿着立管轴向的涡激振动振幅分布图.图12(a)为来流速度为0.3 m/s 情况下,不同长度立管的涡激振动振幅沿着立管轴向分布图,对照图10、图11 可以看出,立管长度为50 m和100 m 时主要以第1 阶模态为主导的振动,立管长度为150 m 时主要以第3 阶模态为主导的振动.如图12(b)所示,顶张力为80 000 N 和100 000 N 激发出以立管第2 阶模态为主导的振动;顶张力为60 000 N时,激发出以立管第3 阶模态为主导的振动.随着顶张力增大,立管振动位移略有增大.这是由于顶张力增大,立管弯曲刚度增大,立管湿模态频率升高,在同等流速下由于涡激频率不变,导致立管激起模态降低,从而使立管振动位移增大.

图12 沿RTP 立管轴向的涡激振动振幅分布Fig.12 Amplitude of vortex induced vibration along the axis of RTP riser

基于MSTMM 快速建模的优点,考虑150 m 长立管中间包括多个刚性接头对涡激振动响应的影响,只需将传递矩阵中的刚性接头部分的RTP 的参数换成刚性接头的参数(EI=9.5734×107N·m2,¯m=284.70 kg/m),即可快速求解出立管振动特性及动力响应.从图13(a)可以看出,均匀添加了3 个刚性接头对RTP 立管的频率影响较小.这是由于添加了刚性接头,同时微弱的增大了立管的质量和刚度,因此立管的湿模态频率变化微弱.从图13(b)可以看出,增加了刚性接头后,对立管的振型也只是产生微弱的影响.同理,均匀添加了2 个刚性接头对RTP 立管的频率和振型几乎没有影响.但是,在立管中均匀添加2 个或3 个刚性接头会对涡激振动的振幅分布产生一定影响.

图13 包含3 个刚性接头情况下沿RTP 的湿模态Fig.13 Wet mode along RTP with three rigid joints

从图14 可以看出,均匀添加了2 个刚性接头后,2 个刚性接头对RTP 立管的涡激振动振幅分布影响较小.从图15 中可以看出,在立管中添加3 个刚性接头后,沿着RTP 立管的涡激振动振幅却有明显变化.均匀分布的3 个刚性接头对立管的低流速激发出以第1 阶模态为主导的涡激振动振幅分布影响较小,对立管中部0.2 m/s 流速激发出以第2 阶模态为主导的涡激振动振幅响应影响稍大.当来流速度为0.3 m/s时,刚性接头对涡激振动振幅分布的影响比较显著.这说明添加了3 个刚性接头后,随着流速的增大,激发出的模态越高,对涡激振动的振幅分布影响越大.

图14 包含两个刚性接头情况下沿RTP 立管轴向的涡激振动振幅分布Fig.14 Amplitude of vortex induced vibration along the axis of RTP riser with two rigid joints

图15 包含3 个刚性接头情况下沿RTP 立管轴向的涡激振动振幅分布Fig.15 Amplitude of vortex induced vibration along the axis of RTP riser with three rigid joints

实际海洋环境中,水流条件复杂多变,立管不仅存在均匀流的作用,更多实际情况下还受到倾斜流、剪切流、阶梯流等作用.如图16 所示,来流为线性剪切流,靠近立管底部流体速度都为0,靠近海平面的三种来流速度情况分别是0.1,0.2,0.3 m/s.从图中可以看出,剪切流顶端来流速度为0.1 和0.2 m/s 主要激发出了以第1 阶模态为主导的振动.剪切流顶端来流速度为0.3 m/s 主要激发出以第2 阶模态为主导的振动.相对于均匀来流速度情况下的沿着立管轴向的涡激振动振幅分布有较大不同.剪切流顶端来流速度为0.1 和0.2 m/s 时引起的涡激振动振幅都比均匀来流情况下引起的涡激振动振幅要小,这主要是剪切流的能量相对于均匀流的能量过小的原因造成的.而当剪切流顶端来流速度为0.3 m/s 时,涡激振动引起的振幅相比0.3 m/s 的均匀来流情况下的振幅稍大.这是由于,剪切流顶端来流速度为0.3 m/s 时激起的立管模态降低,使得立管振动位移增大.