海管S型初始铺设仿真算法

李　震, 张同喜, 孟凡森, 许秀军

(哈尔滨工程大学 机电工程学院， 哈尔滨 150001)

海管S型初始铺设仿真算法

李震, 张同喜, 孟凡森, 许秀军

(哈尔滨工程大学 机电工程学院， 哈尔滨 150001)

摘要:为准确模拟初始铺管作业中管道和起始缆的形态，创建一个逼真的深水铺管起重船铺管作业虚拟训练环境. 针对S型铺管作业中的初始铺管作业，以Euler-Bernoulli梁理论为基础，对初始铺管管道和缆索进行静态分析，建立几何非线性微分方程，且对管道与缆索微分方程边界条件难以确定导致方程无法求解的问题，提出一种基于微分求积法的迭代方法，该方法能够准确地实现边界条件的确立，从而完成微分方程的求解. 仿真和实验分析不同作业状态下管道与缆索的形态与内力变化，结果证明了初始铺管整体算法的准确性. 该方法提高了微分方程求解精度且计算量少，易于程序实现，可用于海上铺管作业方案的工程预演，可行性分析和优化作业方案等.

关键词:S型铺管；Euler-Bernoulli梁；非线性；微分求积法;迭代

对于深水海底管道铺设，管道受力及变形分析一直是铺管作业过程中关注的重点.目前，对正常铺管与收、弃管作业研究居多，同时也涌现出了许多关于管道分析和计算的方法. 最早由Plunkett[1]对海底立管进行分析，并提出用奇异摄动法求解管道微分方程； Dareing等[2]首次将有限差分法应用管道微分方程的求解；Croll[3]首先提出对S型管道整体形态采用自然悬链线法进行分析计算，对微分方程求解过程进行简化，修正边界条件；Dixon等[4]考虑管道刚度，提出刚性悬链线法，并完成管道微分方程的求解；顾永宁[5]针对不同因素影响,分别提出了求解方法，对于刚性悬链线法，将只能应用于下弯段的刚性悬链线法拓展到上弯段，简化了计算；Kalliontzis[6]考虑了管道塑性变形对铺管作业的影响；随着有限元软件的兴起， Khdeir[7]提出采用ABAQUS对管道铺设静态进行分析. 迄今为止，对初始铺管作业过程研究较少，缺少对初始铺管过程中管道与缆索的静态与动态分析. 初始铺管作业是指海底管道在缆索的导引下从铺管船延伸至海床，同时张紧力逐渐增加至设计值的过程. 它是整个海上铺管作业的开始，是正常铺设的必要前提. 为达到仿真程序设计要求，对管道和缆索进行受力分析，针对数学模型实时解算方法进行研究是很有必要的. 而在初始铺管中，要求同时求出管道和缆索的形态，增加了求解的难度. 在初始铺管作业过程中，由于管道剪切变形和横截面的转动非常微小，因此其属于一种大变形梁[8].

本文以“海洋石油201”深水铺管起重船的初始铺管作业为研究对象，对初始铺管作业计算理论进行分析，在原有的铺管常用算法基础上提出了将非线性梁理论[9]应用于初始铺管计算，基于Euler-Bernoulli梁理论[10]建立了管道的初始铺设数学模型，并针对建立的微分方程特点，采用拟牛顿法[11-12]对方程组进行求解，利用迭代的方法求出微分方程边界，得出微分方程的数值解，从而解算出管道和缆索的形态和内力. 微分求积[13]数值算法，使微分方程求解精度高,计算量少,易于程序实现，更有利于视景仿真[14]系统对计算程序的需求. 建立起一套基于海洋石油201 船工程实践的视景仿真系统，不仅可以减少训练的危险性，缩短实船训练周期，还可以到达工程预演的目的，对提升作业效率与安全性,以及降低成本等有着积极重要的作用.

1初始铺管微分方程

1.1管道力学分析

管线未变形前长度为L,直径为D,惯性矩为I,弹性模量为E,管线浮重度 (单位长度管线在水中的重量) w，在管线端部受到水平方向载荷H,竖直方向载荷V. 管线受力见图1.

图1　管线微段受力

图中θ=θ(x)为轴线的切线与x轴的夹角也是横截面的转角. 设弧长为s0=s0(x)，由图1可得如下几何关系：

其中δ为轴线伸长率.

由梁理论和胡克定律[15]可得

其中：TW和M分别为轴力和弯矩，EA和EI分别为抗拉和抗弯刚度.

如图3，管线变形后任意微分单元上有连续分布机械力q=w，利用内力等效关系，轴向内力TW(拉力为正)可表示为

联立上式可以得到轴线伸长率和微分方程为

图2　初始铺管示意

图3　管线整体变形

1.2缆索力学分析

从缆索中取出任意微段i，微段的受力情况如图4所示，假设微段长度为si，两端的受力分别为Ti和Ti+dT，缆索自重简化为作用在微段重心上的集中力，在拉力和自重作用下缆索产生一定的转角dφci.

图4　缆索单元受力

根据缆索微段轴向受力平衡,可以得到以下微分方程:

式中：wc为缆索的浮重度，φc为缆索轴线与水平方向的夹角.

如果忽略缆索在轴向上的伸长量，将缆索分为n段，设由第i点至第i+1点的长度为si，两端的轴力分别为Ti和Ti+1，将上式写为线性积分形式为:

2求解微分方程

2.1微分求积格式

定义一系列量纲一的系数:

式中：L为管线长度，I为惯性矩，E为弹性模量.

则管道大挠度控制微分方程转化为量纲一的形式如下：

(1)

式中:

由微分求积法可得控制微分方程(1)的微分求积格式为[16]

(2)

内力的表达式为

其中1≤i≤N.

2.2逆Broyden法求解非线性方程组

(3)

从式(3)可以看出,有x1,x2,…xN和H、V总计n+2个参数是未知的，而只有n个方程. 假设水平拉力H和竖直拉力V已知，即对非线性方程组(4)的求解：

(4)

(5)

Step 1给定初始近似值x0∈Rn及计算精度要求的ε1与ε2；

Step 2计算初始矩阵H0∈Rnn，可取H0=[F′(x0)]-1或单位矩阵I;

Step 7k=k+1，转Step 4；

Step 8x*=xk+1，结束.

显然,这个过程需要水平拉力H和竖直拉力V已知,即可求解.

2.3利用边界条件确定水平拉力H和竖直拉力V

初始铺管作业过程中，对于管道，如果已知张紧器张力及管道下端所受力水平力H与垂直力V，则可根据管道挠度无量纲方程. 算出管道形态及管道所受弯矩M及轴向力N的. 对于缆索，如果已知一端所受的张力T 和角度φci，则可根据式(2)，从缆索已知端开始逐段计算，直至求出整条缆索的形态和内力. 但是，在实际起始铺管作业过程中管道下端的坐标和缆索所受张力都是未知的，即管道下端的边界条件是未知的. 由于管道下端和缆索相连，根据静力平衡可知缆索和管道内的张力大小相等而方向相反，又因为缆索刚度极小基本可以忽略，所以管道下端所受弯矩为0.

已知的变量有张紧器张力T、管道长度L、缆索长度S、起始锚坐标、船体坐标、船舶六自由度及管道物理参数.

首先，定义一系列数组和变量如下：

Ti为管道第i段拉力；Hi为管道第i段水平拉力；Vi为管道第i段竖直拉力；θi为管道第i段转角；Mi为管道第i段弯矩；其中(i=1,2,…,N)，第N点为管道与缆索连接点(如图2点B). 对于初始铺管过程中每一时刻，定义一个数组DH(10)：用于存储管道下端水平拉力；数组de(10)：用于存储误差；dt为水平拉力平均分成10份，每份大小为dt；e为计算精度.

由管线微段受力(图1)受力分析可得

由管道下端所受弯矩为0，即由边界条件MN=0可得

进行量纲一化得

(6)

通过式(6)以看出，只要给出HN，就能计算出VN，这样减少了迭代过程中的未知变量，使方程组(4)可以求解. 取起始锚位置为绝对坐标系原点，船体位置坐标(x，y，z)由卫星定位可知.

然后判断出10个de(i)中最小的一个，其下标定为I，从而DH(I)对应的误差de(I)是最小的一个，即DH(I)是最接近精确解的.

第一轮计算完毕，进行下一轮的逼近计算. 在DH(IP)附近取10个点组成新的DH(10). 根据DH(IP)值的特点，分3种情况计算：

1)如果DH(I)是10个中最小的，则取DH(1)=DH(I), dt=dt/9，DH(i)=DH(1)+(i-1)dt，i=2,…10.

2)如果DH(I)是10个中最大的，则取DH(1)=DH(I-1), dt=dt/9,DH(i)=DH(1)+(i-1) dt，i=2,…10.

3)如果DH(I)在10个中不是最大也不是最小的，则取DH(1)=DH(I-1)，dt=2dt/9，DH(i)=DH(1)+(i-1) dt，i=2,…10.

这样就定了新一轮计算的水平拉力DH(10). 值得注意的是，为保证每轮计算中都有10个元素，第一轮逼近中总共分割了10份，去掉为0的点，取上面10个点作为输入. 而第二轮计算中由于第一个点不为0，所以分割了9份，取全部10点作为输入.

如此反复进行直到前后两轮计算所得的DH(I)的差值满足精度要求为止，至此求出了管道长度为L时管道下端水平拉力H. 至此就能求解非线性方程组(4)，就可以求出管道和缆索形态及内力分布.

3仿真实例验证及其结果分析

3.1仿真作业算例

为了验证初始铺管作业理论分析的准确性、数值解算方法的快速性及稳定性，铺管作业仿真系统数学模型海况输入采用第3方提供的DP数据. 此次算例水深1 000 m，管道直径323.9 mm，壁厚15.9 mm；缆绳长度为1 000 m，缆绳外径0.108 m，密度7 850 kg/m3；托管架长度和曲率半径分别为87 m和85 m. 为便于分析，提供铺管过程中的6种状态对应的长度及张力(见表1).

表1　起始铺管管道长度及张力

3.2仿真结果分析

将仿真结果与理论数据对比，验证1 000 m水深海底管道起始铺设过程中管道和缆索的形态、轴力及弯矩分布. 仿真的管线形态如下图5(a)～(f)，节点上部分为管道，下部分为缆索形态.

图5中实线代表仿真得到的曲线，虚线为商业软件OFFPIPE计算得到的本次案例曲线，均为静态管道铺设形态. 结果表明，随着铺管作业的进行，管道长度逐渐增加，管道下端高度逐渐降低. 从S1～S3可以看出，管道和缆索形态逐渐变平缓，缆索触底点的位置逐步远离铺管船，而在S3～S6中，缆索触底点重新向铺管船靠近，铺管形态又变得陡峭.

在起始铺管S1～S6中，管道上端的张力从180.0 kN逐渐增加到400.0 kN，铺管作业过程中，起始铺管轴力分布见图6，起始铺管弯矩分布见图7.

由图6、7可以看出，S1～S6过程中，缆索最大轴力逐渐变小，轴力变化趋势逐渐变缓，S1时最陡，S6时最缓；缆索内没有弯矩，这是由于假设缆索的抗弯刚度为零.

(a) S1

(c) S3

(e) S5

(b) S2

(d) S4

(f) S6

图6　起始铺管轴力分布

图7 起始铺管弯矩分布

4结论

1)对初始铺管作业过程中的管道与缆索进行受力分析. 在原有的铺管常用算法基础上提出了将非线性弹性大挠度梁理论应用于初始铺管计算，建立管道数学模型，采用微分求积法求解，同时建立缆索模型，并用迭代法求解. 微分求积法极大地减少了计算量且求解精度较高，求解过程易于在计算机上程序实现，提高实时仿真的真实度，有助于视景仿真系统的开发.

2)在求解微分过程中，基于管道与缆索微分方程边界条件难以确定导致方程无法求解的问题，采用一种基于微分求积法的迭代方法，完成边界条件的确立,然后用拟牛顿法完成方程组的求解. 利用已知的边界条件迭代求出所需的边界条件，在一定的误差范围内结束迭代，结果表明，此方法有效解决了边界条件问题.

3)仿真得到的管道和缆索形态曲线与项目组提供的本次案例曲线的比较表明，仿真求解得到的曲线与理论曲线基本一致，表明作业理论分析的准确性高,数值方法解算速度快, 稳定性好.

参考文献

[1] PLUNKETT R. Static bending stresses in catenaries and drill strings [J]. Journal of Manufacturing Science and Engineering, 1967, 89(1): 31-36.

[2] DAREING D W, NEATHERY R F. Marine pipeline analysis based on Newton’s method with an arctic application[J]. Journal of Manufacturing Science and Engineering, 1970, 92(4): 827-833.

[3] CROLL J G A. Bending boundary layers in tensioned cables and rods[J]. Applied ocean research, 2000, 22(4): 241-253.

[4] DIXON D A, RUTLEDGE D R. Stiffened catenary calculations in pipeline laying problem[J]. Journal of Manufacturing Science and Engineering, 1968, 90(1): 153-160.

[5] 顾永宁. 海底管线铺管作业状态分析[J]. 海洋工程, 1988, 6(2): 11-23.

[6] KALLIONTZIS C. Non-linear finite element simulations of highly curved submarine pipelines[J]. Communications in numerical methods in engineering, 1998, 14(11): 1067-1088.

[7]KHDEIR A A. Thermal buckling of cross-ply laminated composite beams[J]. Acta Mechanica, 2001, 149(1/2/3/4):201-213.

[8] DABROWSKI R. Curved thin-walled girders: Theory and Analysis[R]. Workingham: Transport and Road Research Laboratory, 1900.

[9] 花雷.曲梁结构非线性大变形分析[D].扬州:扬州大学,2011.

[10]李世荣,孙云,刘平.关于Euler-Bernoulli梁几何非线性方程的讨论[J].力学与实践,2013,35(2):77-80.

[12]ZIANI M, GUYOMARC’H F. An autoadaptative limited memory Broyden’s method to solve systems of nonlinear equations[J].Applied Mathematics and Computation, 2008, 205(1): 202-211.

[13]王鑫伟. 微分求积法在结构力学中的应用[J]. 力学进展, 1995, 25(2): 232-240.

[14]BURDEA G, COIFFET P. Virtual reality technology[J]. Presence: Teleoperators and virtual environments, 2003, 12(6): 663-664.

[15]黄玉盈, 朱达善. 海洋管线铺设时的静力分析[J]. 海洋工程, 1986, 4(1): 32-46.

[16]聂国隽, 仲政. 微分求积单元法在结构工程中的应用[J]. 力学季刊, 2005, 26(3): 423-427.

[17]马成业，黄世华. 解非线性方程组的一个改进牛顿法[J]. 甘肃联合大学学报(自然科学版), 2009, 23(1): 116-117.

(编辑杨波)

Simulation algorithm for initial S-lay

LI Zhen, ZHANG Tongxi, MENG Fansen, XU Xiujun

(College of Mechanic and Electrical Engineering, Harbin Engineering University, Harbin 150001, China)

Abstract:To simulate the shape of the pipe and cable in the initial pipe laying operation accurately, a realistic virtual training environment of pipe-laying operation was created. For the initial S-laying and based on the Euler-Bernoulli beam theory, the geometric nonlinear differential equations is established by analyzing the static force of initial pipe laying and cable. To solve the problem of the differential equations boundary conditions, an iterative method based on differential quadrature method is proposed. By simulation and experiments, the shape and variation of the internal force of the pipe and cable under different operating conditions are analyzed, and the accuracy of the algorithm for initial S-laying is proved. The practical example shows that the algorithm can be applied to the initial S-laying effectively. The accuracy of the differential equation is improved and it is easy to be programmed by using this method, which contributes to the preview of offshore pipe laying operation , feasibility analysis and optimization of operation plan.

Keywords:S-lay; Euler-Bernoulli beam; nonlinear; differential quadrature method; iterative

doi：10.11918/j.issn.0367-6234.2016.07.017

收稿日期:2015-05-04

基金项目:国家科技重大专项 (Z12SJENA0014)

作者简介:李震(1971—)，男，副教授，硕士生导师

通信作者:张同喜,yanyulou126@126.com

中图分类号:TH133; TP183

文献标志码:A

文章编号:0367-6234(2016)07-0106-06