基于向量式有限元的深水顶张力立管动力响应分析

李效民， 张　林， 牛建杰， 韩永庆， 郭海燕

(中国海洋大学 工程学院，山东　青岛　266100)

基于向量式有限元的深水顶张力立管动力响应分析

李效民， 张林， 牛建杰， 韩永庆， 郭海燕

(中国海洋大学 工程学院，山东青岛266100)

向量式有限元法亦称有限质点法是以向量式力学与数值计算为基础的一种结构行为分析的创新性方法。该方法将结构离散为一系列质点的集合，质点间采用单元连接，单元没有质量，只承受内力，质点的运动则直接采用牛顿第二定律来表达，采用中央差分的显式积分法求解质点群各时刻运动方程。首次将该方法应用于顶张力立管的动力行为分析中，采用平面弯曲杆件元模拟质点间的相互作用力，推导了立管质点运动公式，编制相应Matlab求解程序，计算了立管在海流、波浪作用下的动力响应特性并与基于传统有限元法的结果进行了对比分析，表明该方法在立管行为分析中的正确和有效性，并指出了该方法在处理其它立管形式行为分析中的优势。

向量式有限元(VFIFE)；有限元法(FEM)；向量式力学；顶张力立管(TTR)

顶张力立管(TTR)是各种形式的海洋立管中既可用于钻井、完井及修井又可以用于注水、生产等功能的立管，是油气勘探开发中必不可少的立管形式。近几十年来受到了国内外科技界和工业界的广泛关注，对其展开了深入研究。总体来说，研究分析主要集中在涡激振动、参数激励(稳定性)、顺流向振动响应研究等几方面，并采用理论分析、试验研究与数值模拟相结合的方式展开相关研究。

Morkookaza等[1]在时域和频域内对比分析了顶张力立管在不同波浪力和海洋环境条件下的动力响应。Chatjigeorgiou等[2]研究了垂直海洋立管在其顶端承受参数激励下的非线性响应特性，详细分析了水动力阻尼对系统动力不稳定区域的影响。Kuiper等[3]研究了深水立管参数激励的稳定性，表明浮式结构的升沉运动能够导致立管出现不稳定性。Chaplin等[4]为增加立管涡激振动物理特性的理解并得到详实的试验数据来验证涡激振动预测模型，在Delta 水槽中设计了一精巧的顶张力立管涡激振动实验研究方案。高云等[5]比较系统的研究了不同鳍高、不同螺距的Strake下立管状态下的响应特性。王腾等[6]模拟了一钻井隔水管在波、流和平台运动联合作用下的非线性动力响应，研究了平台的静偏移和动态运动对隔水管最大弯矩分布的影响。畅元江等[7]考虑波流联合作用下隔水管与导向架之间的间隙接触条件，用ANSYS研究了非线性边界条件下隔水管的动态响应。杨和振等[8]采用Floquet理论分析了深水立管参数稳定性，指出了参数共振区域，分析了阻尼对参数激励不稳定性的影响。唐友刚等[9]研究了顶张力立管在剪切流中的涡激振动响应和参激-涡激耦合振动响应特性，分析参数激励对立管横向涡激振动的影响。陈伟民等[10]采用有限元法对浮体升沉与立管涡激振动的耦合动力响应进行了分析。唐国强等[11]建立了一种立管涡激振动的经验模型，研究了立管顶部张力、水流速度、立管外径以及立管管壁厚度的变化对于立管动力响应的影响。

综上所述，当前人们对顶张力立管的研究大都是采用基于小变形小位移理论的梁单元的传统有限元方法进行模拟。在立管的运动响应计算过程中刚度矩阵不随时间变化，也即始终以其初始垂直位形作为立管的计算刚度矩阵。对于深水立管来说，当其相对位移较小时该方法是足够精确的，然而当立管在海流和顶端较大浮体偏移下的运动响应计算误差则相对偏大，因为严格意义来说，此时立管运动行为属于小变形大变位的几何非线性问题。对于立管的大变位问题可以采用在应变计算中包含高阶项、推导平衡方程时考虑变形的影响、Lagrange(T.L.)法、修正的Lagrange(U.L.)法和动坐标法/共旋坐标法等诸多方法来求解，不过对其进行建模求解需要特殊技巧和方案，同时计算中难以区分立管单元的刚体运动和纯变形，此时对立管动态响应求解将变得及其复杂。

向量式有限元法(VFIFE)亦被称作有限质点法(FPM)是由丁承先等[12]基于向量力学与数值计算相结合提出的一种新型结构行为分析方法，是求解结构大变形、大变位、弹塑性、碰撞等非线性或不连续性力学行为的新方法。与传统有限元需要建立完整的数学理论描述不同，向量式有限元法以结构的物理模式为基础，将结构离散为质点群，质点间采用单元连接，单元没有质量，只承受内力，利用牛顿第二定律描述这些质点的运动，采用途径单元来描述质点的运动，采用虚拟的逆向运动计算结构的纯变形，可以同时计算结构的运动和变形，而不采用传统有限元方法使用的连体控制方程描述结构。倪秋斌等[13]基于向量式有限元法考虑拉索几何和材料非线性对斜拉桥振动控制进行了数值模拟。喻莹等[14-15]基于有限质点法对结构的倒塌破坏过程进行了模拟和分析。罗尧治等[16]阐述了有限质点法这种新的数值分析方法在空间结构复杂行为研究领域的优势及应用，并对该方法的发展趋势作出了展望。班自愿等[17]对悬臂梁分别承受集中载荷和弯矩下的大变形进行算例分析。袁峰[18]采用弹性海床和塑性海床模拟管土相互作用，将向量式有限元法初步应用于海洋管土动力循环中，取得了一些有益的分析结果。

本文首次将向量式有限元法应用于深水顶张力立管的动力响应研究中。将立管离散为一群质点的集合，采用平面弯曲杆件元模拟质点间的相互作用力，建立了立管质点运动计算公式，编制相应Matlab求解程序，计算了深水立管在海流、波浪作用下的动力响应特性并与基于传统有限元法(FEM)的结果进行了对比分析。通过对比分析，验证向量式有限元法在深水顶张力立管动力行为分析中的有效性及其优势，并为进一步将这种新型数值计算方法在深水各种形式的海洋立管平面内外复杂环境载荷条件下的行为分析中的应用前景作出展望。

1立管质点运动方程的建立及载荷

1.1立管质点的点值描述

立管整体坐标系及其振动变形及坐标示意图如图1所示。假设顶张力立管在初始位置时完全垂直。设定一组域坐标(x，y，z)，基底为(i，j，k)，如图1(a)所示。坐标系定义如下：以立管未发生变形时垂直轴为y轴，向上为正，并取立管底部为原点坐标；x轴与来流方向平行，z轴垂直于x和y轴。 将立管分为N段共N+1(从底向上编号依次为1, 2，… ，I,J, …,N,N+1)个空间点来进行描述。等效节点质量用tn-1时立管的密度、横截面积、跨度以及内部流体的密度、断面面积、计算跨度确定。初始时刻各点沿x,y向平移速度、转角和旋转速度均为0。

图1　立管振动变形及坐标示意图Fig.1 Diagram of the riser vibration amd coordinate

(1)

1.2节点变形

(2)

(3)

(4)

θI=βI-ΔθθJ=βJ-Δθ

(5)

1.3内力计算

然后将I、J节点内力向量分别转换成域坐标分量，节点弯矩依然不需要转换。

(6)

最后对每一个运动质点作内力的集成。

1.4质点控制方程

考虑结构阻尼，各质点运动遵循牛顿第二定律，控制方程为：

(7)

对于立管平面振动问题，其各质点有三个位移分量包括x，y轴向的平移分量和绕z轴的转动。上式具体展开表达式为：

(8)

1.5水动力载荷

实际作用于立管上的波浪、海流载荷是一种分布载荷，可以采用传统有限元的处理方法将其等效到立管各质点上。对于单位长度立管所受到的作用力可以采用Morison方程进行计算。

(9)

2数值求解

采用中央差分的显式时间积分法进行求解，可以避免隐式解法带来的迭代和收敛问题。但是时间步长必须进行控制，否则也会导致计算错误。对于每一个质点运动公式的计算采用第n-1步的运动公式，代入差分公式计算第n步的点位置值xn。然后通过计算出的第n步点位置值xn以及第n步的内力公式来计算第n步的节点内力和等效力，如此循环往复直到模拟计算时刻最后一步。

假设时间步长为h，那么质点的速度和加速度分别可以表示为

(10)

将其代入质点运动公式，得到其对应的差分计算公式为：

(11)

式中C1=1/(1+ζh/2)，C2=C1(1-ζh/2)。

3程序编制及算例分析

基于上述理论和计算方法，编制相应的Matlab求解程序，对深水立管动态响应进行数值模拟，并将其结果与以小变形理论为基础的传统有限元模拟结果进行对比分析。传统有限元模拟中选取梁单元模拟立管，对立管运动方程采用Galerkin法进行有限元离散和Newmark-β法进行数值求解。本文选取两种具有不同长度和顶张力的立管进行模拟。立管一：水深300 m，长度300 m，顶端张力500 kN，立管二：水深1 000 m，长度1 000 m，顶端张力1 500 kN，其他主要参数相同：外径0.26 m，内径0.2 m，海水密度为1 025 kg/m3, 内部流体密度为800 kg/m3，立管管材密度7 850 kg/m3，弹性模量2.07×11N/m4，曳力系数0.7，惯性系数2。分析中两种立管上下两端均假设为铰接，划分单元数N都为60。

图2和图3绘制了两种立管在海流流速分别为0.1 m/s，0.2 m/s，0.3 m/s和0.4 m/s四种工况下的最大位移图。从两图中可以看出当流速较小时基于向量式有限元方法(VFIFE)计算的结果与基于传统有限元法(FEM)计算的结果相当吻合，而当流速增大时，两者之间的差异逐渐变大，并且前者始终小于后者。其主要原因在于向量式有限元方法在变形计算中考虑了立管张力的变化，随着海流流速的增加，立管位移增大，管内轴力也随之增加，四种工况下立管一的顶端张力分别为500.3 kN、504 kN、518 kN和544 kN，立管二顶端张力分别为1 501 kN、1 510 kN、1 541 kN和1 598 kN，轴力的增大使得立管位移有所减小。

图2　不同流速下立管一最大位移图Fig.2 The maximum displacement envelope of the first model riser under different currentvelocity

图3　不同流速下立管二最大位移图Fig.3 The maximum displacement envelope of the second model riser under different current velocity

图4和图5分别绘制了立管一和立管二在波高3 m周期10 s线性波浪激励下的最大位移包络图。从图4中可以看出立管一基于向量式有限元方法计算的结果与基于传统有限元法计算的结果偏小，而从图5中可以看出立管二基于向量式有限元方法计算的结果与基于传统有限元法计算的结果差异可以忽略。

图4　波浪作用下立管一最大位移包络图Fig.4 The maximum displacement envelope of the first model riser under wave excited

图5　波浪作用下立管二最大位移包络图Fig.5 The maximum displacement envelope of the second model riser under wave excited

图6　立管一100 m和200 m处位移时程及其频谱Fig.6 Displacement history and frequency spectrum of the first model riser at 100 m and 300 m

图7　立管一底端和顶端张力时程及其频谱图 Fig.7 Bottom and top end tension history and frequency spectrum of the first model riser

图6～9分别绘制了立管一和立管二两个节点处位移时程曲线及其频谱图和立管底端和顶端张力时程曲线及其频谱图。从图6和图8中可以看出立管两个节点处的水平振动频率和横向激励频率保持一致，从图7和图9中可以看出轴力是随着时间周期变化的，频率为0.2 Hz是其横向激励频率的两倍，这与有限元分析中轴力保持不变是不同的，显然轴力的变化与立管实际受力状态更加符合。立管一初始状态下立管底端轴力是85.4 kN 顶端是500 kN，振动过程中立管底端和顶端轴力最大分别都增加12 kN左右，分别相对增大14%和2.4%，正是轴力的较大变化导致基于向量式有限元计算得到的振动包络图比基于传统有限元的图4结果偏小。立管二初始状态下底端轴力是117.5 kN顶端是1 500 kN，振动过程中立管底端和顶端轴力最大分别都增加5.5 kN左右，分别相对增大4.7%和0.37%，此时立管轴力变化较小从而基于向量式有限元计算得到的振动包络图与基于传统有限元的图5结果基本一致。

图8　立管二300 m和700 m处位移时程及其频谱图Fig.8 Displacement history and frequency spectrum of the second model riser at 300 m and 700 m

图9　立管二底端和顶端张力时程及其频谱图Fig.9 Bottom and top end tension history and frequency spectrum of the second model riser

4结论

对于深水顶张力立管的动力特性及动力响应分析方面，与以小变形理论为基础的传统有限元分析相比，本文基于向量式有限元方法分析有以下结论：

(1) 向量式有限元方法分析中考虑了立管在变形过程中轴向运动及其引起的轴向张力的变化，其振动频率为横向激励频率的两倍，这与立管实际受力状况更加符合。而有限元分析中一般不考虑立管变形过程中轴向张力的变化，也就是始终等于立管静止状态下初始有效张力。

(2) 当海流流速较小时，立管整体变形较小，轴力改变很小，向量式有限元法与有限元结果基本一致。随着流速增大，立管整体变形增大，轴力也随之增大，导致向量式有限元法结果比有限元结果偏小。

(3) 对于波浪激励下的立管，当立管相对较短时，轴力波动幅值相对有效张力较大，向量式有限元方法比有限元结果偏小。而当立管相对较长时，轴力波动幅值相对有效张力较小，向量式有限元法与有限元结果基本一致。

通过本文对深水顶张力立管的向量式有限元理论分析和数值模拟，可以看出向量式有限元法能够真实准确地模拟立管的振动行为，在立管动态行为分析中是正确和有效性，并且具有其独特的优势，点运动公式的描述方式使得控制方程大大简化，无需求解复杂的非线性联立方程组，采用中央差分的显式时间积分求解避免了迭代和收敛问题，能够同时计算立管的内力(应力应变)和变形。

总之向量式有限元法在处理其它各种形式的悬链线立管和柔性立管中将同样具有显著优势，可以采用统一的架构处理不同的立管结构形式及其复杂的力学特征(小变形大变位、大变形大变位、接触碰撞等)，为进一步对各种形式的海洋立管复杂行为分析提供了一种新的可行的求解方案。

[ 1 ] Morkookaza C K, Coelho F M, Shiguemoto D A. Dynamic behavior of a top tensioned riser in frequency and time domain [C]//Proceedings of the Sixteenth International Offshore and Polar Engineering Conference, ISOPE, San Francisco, 2006, 2:31-36.

[ 2 ] Chatjigeorgiou I K, Mavrakos S A. Bounded and unbounded coupled transverse response of parametrically excited vertical marine risers and tensioned cable legs for marine applications [J].Applied Ocean Research, 2002, 24(6):341-354.

[ 3 ] Kuiper G L, Brugmans J, Metrikine A V. Destabilization of deep-water risers by a heaving platform[J]. Journal of Sound and Vibration, 2008, 310(3):541-557.

[ 4 ] Chaplin J R, Bearman P W, Huera Huarte F J,et al. Laboratory measurements of vortex-induced vibrations of a vertical tension riser in a stepped current[J]. Journal of Fluids and Structures, 2005, 21(1):3-24.

[ 5 ] 高云,付世晓,宋磊建.柔性立管涡激振动抑制装置试验研究[J].振动与冲,2014,33(14):77-83.GAO Yun,FU Shi-xiao,SONG Lei-jian.Experimental investigation on the suppression device of VIV of a flexible riser[J].Journal of Vibration and Shock,2014,33(14):77-83.[ 6 ] 王腾，张修占，朱为全. 平台运动下深水钻井隔水管非线性动力响应研究[J]. 海洋工程，2008, 26(3):21-26.WANG Teng, ZHANG Xiu-zhan, ZHU Wei-quan. Study on nonlinear dynamic response of deepwater drilling riser with vessel motion[J].The Ocean Engineering,2008, 26(3):21-26.

[ 7 ] 畅元江，陈国明，许亮斌. 导向架隔水管在波流联合作用下的非线性动力响应[J].中国石油大学学报,2006, 30(5):74-77.

CHANG Yuan-jiang, CHEN Guo-ming, XU Liang-bin. Nonlinear dynamic response induced by wave-current for marine risers with guide-frames[J].Journal of China University of Petroleum, 2006, 30(5):74-77.

[ 8 ] Yang He-zhen, Li Hua-jun. Instability Assessment of deep-sea risers under parametric excitation[J].China Ocean Engineering, 2009, 23(4):603-612.

[ 9 ] 唐友刚，邵卫东，张杰，等.深海顶张力立管参激-涡激耦合振动响应分析[J].工程力学，2013,30(5):282-286.

TANG You-gang, SHAO Wei-dong, ZHANG Jie,et al. Dynamic response analysis for coupled parametric vibration and vortex-induced vibration of top-tensioned riser in deep-sea[J]. Engineering Mechanics,2013,30(5):282-286.

[10] Chen Wei-min, Li Min, Guo Shuang-xi,et al. Dynamic analysis of coupling between floating top-end heave and riser’s vortex induced vibration by using finite element simulations[J]. Applied Ocean Research, 2014,48:1-9.

[11] 唐国强，滕斌，吕林，等.深海柔性立管涡激振动经验模型建立及应用[J].中国海洋平台, 2010，25(3):12-16.

TANG Guo-qiang, TENG Bin, LÜ Lin,et al. Development and application of an empirical model for vortex-induced vibration of flxible deep water riser[J].China Offshore Platform, 2010, 25(3):12-16.

[12] 丁承先，段元锋，吴东岳.向量式结构力学[M].北京:科学出版社，2012.

[13] 倪秋斌，段元锋，高博青.采用向量式有限元的斜拉索振动控制仿真[J].振动工程学报，2014, 27(2):238-245.

NI Qiu-bin，DUAN Yuan-feng，GAO Bo-qing.Vector form intrinsic finite element (VFIFE) based simulation on vibration control of stay cables[J]. Journal of Vibration Engineering,2014, 27(2):238-245.

[14] 喻莹，罗尧治.基于有限质点法的结构倒塌破坏研究Ⅰ：基本方法[J].建筑结构学报，2011,32(11):17-26.

YU Ying, LUO Yao-zhi.Collapse analysis based on finite particle method I: basic approach[J].Journal of Building Structures, 2011,32(11):17-26.

[15] 喻莹，罗尧治.基于有限质点法的结构倒塌破坏研究Ⅱ：关键问题与数值算例[J].建筑结构学报，2011,32(11):27-35.

YU Ying, LUO Yao-zhi.Collapse analysis based on finite particle methodmethodⅡ:key problems and numerical examples[J]. Journal of Building Structures,2011,32(11):27-35.

[16] 罗尧治，郑延丰，杨超，等.结构复杂行为分析的有限质点法研究综述[J].工程力学，2014，31(8):1-7.

LUO Yao-zhi, ZHENG Yan-feng, YANG Chao, et al. Review of the finite partical method for complex behavior of structures[J]. Engineering Mechanics,2014,31(8):1-7.

[17] 班自愿，张若京. 悬臂梁大变形的向量式有限元分析[J].计算机辅助工程，2010,19(4): 25-28.

BAN Zi-yuan,ZHANG Ruo-jing.Vector form intrinsic finite element analysis on large deformation of cantilever beam[J]. Coputer Aided Engineering, 2010,19(4): 25-28.

[18] 袁峰.深海管道铺设及在位稳定性分析[D].杭州：浙江大学，2013.

Dynamic response of a deep-sea top tensioned riser based on vector form intrinsic finite element

LI Xiao-min, ZHANG Lin, NIU Jian-jie, HAN Yong-qing, GUO Hai-yan

(College of Engineering, Ocean University of China, Qingdao 266100, China)

The vector form intrinsic finite element (VFIFE) method also called the finite particle method (FPM) is an innovative method in the structural behavior analysis based on the vector mechanics and numerical calculations. With this method, a structure was discretized into a series of mass particles connected by massless elements only bearing internal forces. The motion of mass paticles obeyed Newton’s second law. Their motion equations were solved using the explicit time domain integration and the central difference method. VFIFE was adoped for the first time to analyze the dynamic behavior of a top tensioned riser (TTR) and the plane beam element was adoped to simulate interaction forces between mass particles. The particle motion equation of the riser was derived and the corresponding solving program was compiled with Matlab. The dynamic behavior of the riser under sea current and wave excitations was calculated and the results were compared with those with the finite element method (FEM). The results showed that the VFIFE method is correct and effective to analyze the dynamic behavior analysis of a riser and it has more advantages than traditional methods do in riser complex behavior analyses.

vector form intrinsic finite element (VFIFE) method; finite element method (FEM); vector mechanics; top tensioned riser (TTR)

10.13465/j.cnki.jvs.2016.11.035

国家自然科学基金(51279187);山东省优秀中青年科学家科研奖励基金(BS2013HZ014)

2015-04-23修改稿收到日期：2015-06-23

李效民 男，博士，讲师，1982年3月生

TU311.3

A