海管S型初始鋪設(shè)仿真算法

李　震, 張同喜, 孟凡森, 許秀軍

(哈爾濱工程大學(xué) 機(jī)電工程學(xué)院， 哈爾濱 150001)

海管S型初始鋪設(shè)仿真算法

李震, 張同喜, 孟凡森, 許秀軍

(哈爾濱工程大學(xué) 機(jī)電工程學(xué)院， 哈爾濱 150001)

摘要:為準(zhǔn)確模擬初始鋪管作業(yè)中管道和起始纜的形態(tài)，創(chuàng)建一個(gè)逼真的深水鋪管起重船鋪管作業(yè)虛擬訓(xùn)練環(huán)境. 針對(duì)S型鋪管作業(yè)中的初始鋪管作業(yè)，以Euler-Bernoulli梁理論為基礎(chǔ)，對(duì)初始鋪管管道和纜索進(jìn)行靜態(tài)分析，建立幾何非線性微分方程，且對(duì)管道與纜索微分方程邊界條件難以確定導(dǎo)致方程無(wú)法求解的問(wèn)題，提出一種基于微分求積法的迭代方法，該方法能夠準(zhǔn)確地實(shí)現(xiàn)邊界條件的確立，從而完成微分方程的求解. 仿真和實(shí)驗(yàn)分析不同作業(yè)狀態(tài)下管道與纜索的形態(tài)與內(nèi)力變化，結(jié)果證明了初始鋪管整體算法的準(zhǔn)確性. 該方法提高了微分方程求解精度且計(jì)算量少，易于程序?qū)崿F(xiàn)，可用于海上鋪管作業(yè)方案的工程預(yù)演，可行性分析和優(yōu)化作業(yè)方案等.

關(guān)鍵詞:S型鋪管；Euler-Bernoulli梁；非線性；微分求積法;迭代

對(duì)于深水海底管道鋪設(shè)，管道受力及變形分析一直是鋪管作業(yè)過(guò)程中關(guān)注的重點(diǎn).目前，對(duì)正常鋪管與收、棄管作業(yè)研究居多，同時(shí)也涌現(xiàn)出了許多關(guān)于管道分析和計(jì)算的方法. 最早由Plunkett[1]對(duì)海底立管進(jìn)行分析，并提出用奇異攝動(dòng)法求解管道微分方程； Dareing等[2]首次將有限差分法應(yīng)用管道微分方程的求解；Croll[3]首先提出對(duì)S型管道整體形態(tài)采用自然懸鏈線法進(jìn)行分析計(jì)算，對(duì)微分方程求解過(guò)程進(jìn)行簡(jiǎn)化，修正邊界條件；Dixon等[4]考慮管道剛度，提出剛性懸鏈線法，并完成管道微分方程的求解；顧永寧[5]針對(duì)不同因素影響,分別提出了求解方法，對(duì)于剛性懸鏈線法，將只能應(yīng)用于下彎段的剛性懸鏈線法拓展到上彎段，簡(jiǎn)化了計(jì)算；Kalliontzis[6]考慮了管道塑性變形對(duì)鋪管作業(yè)的影響；隨著有限元軟件的興起， Khdeir[7]提出采用ABAQUS對(duì)管道鋪設(shè)靜態(tài)進(jìn)行分析. 迄今為止，對(duì)初始鋪管作業(yè)過(guò)程研究較少，缺少對(duì)初始鋪管過(guò)程中管道與纜索的靜態(tài)與動(dòng)態(tài)分析. 初始鋪管作業(yè)是指海底管道在纜索的導(dǎo)引下從鋪管船延伸至海床，同時(shí)張緊力逐漸增加至設(shè)計(jì)值的過(guò)程. 它是整個(gè)海上鋪管作業(yè)的開始，是正常鋪設(shè)的必要前提. 為達(dá)到仿真程序設(shè)計(jì)要求，對(duì)管道和纜索進(jìn)行受力分析，針對(duì)數(shù)學(xué)模型實(shí)時(shí)解算方法進(jìn)行研究是很有必要的. 而在初始鋪管中，要求同時(shí)求出管道和纜索的形態(tài)，增加了求解的難度. 在初始鋪管作業(yè)過(guò)程中，由于管道剪切變形和橫截面的轉(zhuǎn)動(dòng)非常微小，因此其屬于一種大變形梁[8].

本文以“海洋石油201”深水鋪管起重船的初始鋪管作業(yè)為研究對(duì)象，對(duì)初始鋪管作業(yè)計(jì)算理論進(jìn)行分析，在原有的鋪管常用算法基礎(chǔ)上提出了將非線性梁理論[9]應(yīng)用于初始鋪管計(jì)算，基于Euler-Bernoulli梁理論[10]建立了管道的初始鋪設(shè)數(shù)學(xué)模型，并針對(duì)建立的微分方程特點(diǎn)，采用擬牛頓法[11-12]對(duì)方程組進(jìn)行求解，利用迭代的方法求出微分方程邊界，得出微分方程的數(shù)值解，從而解算出管道和纜索的形態(tài)和內(nèi)力. 微分求積[13]數(shù)值算法，使微分方程求解精度高,計(jì)算量少,易于程序?qū)崿F(xiàn)，更有利于視景仿真[14]系統(tǒng)對(duì)計(jì)算程序的需求. 建立起一套基于海洋石油201 船工程實(shí)踐的視景仿真系統(tǒng)，不僅可以減少訓(xùn)練的危險(xiǎn)性，縮短實(shí)船訓(xùn)練周期，還可以到達(dá)工程預(yù)演的目的，對(duì)提升作業(yè)效率與安全性,以及降低成本等有著積極重要的作用.

1初始鋪管微分方程

1.1管道力學(xué)分析

管線未變形前長(zhǎng)度為L(zhǎng),直徑為D,慣性矩為I,彈性模量為E,管線浮重度 (單位長(zhǎng)度管線在水中的重量) w，在管線端部受到水平方向載荷H,豎直方向載荷V. 管線受力見圖1.

圖1　管線微段受力

圖中θ=θ(x)為軸線的切線與x軸的夾角也是橫截面的轉(zhuǎn)角. 設(shè)弧長(zhǎng)為s0=s0(x)，由圖1可得如下幾何關(guān)系：

其中δ為軸線伸長(zhǎng)率.

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

其中：TW和M分別為軸力和彎矩，EA和EI分別為抗拉和抗彎剛度.

如圖3，管線變形后任意微分單元上有連續(xù)分布機(jī)械力q=w，利用內(nèi)力等效關(guān)系，軸向內(nèi)力TW(拉力為正)可表示為

聯(lián)立上式可以得到軸線伸長(zhǎng)率和微分方程為

圖2　初始鋪管示意

圖3　管線整體變形

1.2纜索力學(xué)分析

從纜索中取出任意微段i，微段的受力情況如圖4所示，假設(shè)微段長(zhǎng)度為si，兩端的受力分別為Ti和Ti+dT，纜索自重簡(jiǎn)化為作用在微段重心上的集中力，在拉力和自重作用下纜索產(chǎn)生一定的轉(zhuǎn)角dφci.

圖4　纜索單元受力

根據(jù)纜索微段軸向受力平衡,可以得到以下微分方程:

式中：wc為纜索的浮重度，φc為纜索軸線與水平方向的夾角.

如果忽略纜索在軸向上的伸長(zhǎng)量，將纜索分為n段，設(shè)由第i點(diǎn)至第i+1點(diǎn)的長(zhǎng)度為si，兩端的軸力分別為Ti和Ti+1，將上式寫為線性積分形式為:

2求解微分方程

2.1微分求積格式

定義一系列量綱一的系數(shù):

式中：L為管線長(zhǎng)度，I為慣性矩，E為彈性模量.

則管道大撓度控制微分方程轉(zhuǎn)化為量綱一的形式如下：

(1)

式中:

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

(2)

內(nèi)力的表達(dá)式為

其中1≤i≤N.

2.2逆Broyden法求解非線性方程組

(3)

從式(3)可以看出,有x1,x2,…xN和H、V總計(jì)n+2個(gè)參數(shù)是未知的，而只有n個(gè)方程. 假設(shè)水平拉力H和豎直拉力V已知，即對(duì)非線性方程組(4)的求解：

(4)

(5)

Step 1給定初始近似值x0∈Rn及計(jì)算精度要求的ε1與ε2；

Step 2計(jì)算初始矩陣H0∈Rnn，可取H0=[F′(x0)]-1或單位矩陣I;

Step 7k=k+1，轉(zhuǎn)Step 4；

Step 8x*=xk+1，結(jié)束.

顯然,這個(gè)過(guò)程需要水平拉力H和豎直拉力V已知,即可求解.

2.3利用邊界條件確定水平拉力H和豎直拉力V

初始鋪管作業(yè)過(guò)程中，對(duì)于管道，如果已知張緊器張力及管道下端所受力水平力H與垂直力V，則可根據(jù)管道撓度無(wú)量綱方程. 算出管道形態(tài)及管道所受彎矩M及軸向力N的. 對(duì)于纜索，如果已知一端所受的張力T 和角度φci，則可根據(jù)式(2)，從纜索已知端開始逐段計(jì)算，直至求出整條纜索的形態(tài)和內(nèi)力. 但是，在實(shí)際起始鋪管作業(yè)過(guò)程中管道下端的坐標(biāo)和纜索所受張力都是未知的，即管道下端的邊界條件是未知的. 由于管道下端和纜索相連，根據(jù)靜力平衡可知纜索和管道內(nèi)的張力大小相等而方向相反，又因?yàn)槔|索剛度極小基本可以忽略，所以管道下端所受彎矩為0.

已知的變量有張緊器張力T、管道長(zhǎng)度L、纜索長(zhǎng)度S、起始錨坐標(biāo)、船體坐標(biāo)、船舶六自由度及管道物理參數(shù).

首先，定義一系列數(shù)組和變量如下：

Ti為管道第i段拉力；Hi為管道第i段水平拉力；Vi為管道第i段豎直拉力；θi為管道第i段轉(zhuǎn)角；Mi為管道第i段彎矩；其中(i=1,2,…,N)，第N點(diǎn)為管道與纜索連接點(diǎn)(如圖2點(diǎn)B). 對(duì)于初始鋪管過(guò)程中每一時(shí)刻，定義一個(gè)數(shù)組DH(10)：用于存儲(chǔ)管道下端水平拉力；數(shù)組de(10)：用于存儲(chǔ)誤差；dt為水平拉力平均分成10份，每份大小為dt；e為計(jì)算精度.

由管線微段受力(圖1)受力分析可得

由管道下端所受彎矩為0，即由邊界條件MN=0可得

進(jìn)行量綱一化得

(6)

通過(guò)式(6)以看出，只要給出HN，就能計(jì)算出VN，這樣減少了迭代過(guò)程中的未知變量，使方程組(4)可以求解. 取起始錨位置為絕對(duì)坐標(biāo)系原點(diǎn)，船體位置坐標(biāo)(x，y，z)由衛(wèi)星定位可知.

然后判斷出10個(gè)de(i)中最小的一個(gè)，其下標(biāo)定為I，從而DH(I)對(duì)應(yīng)的誤差de(I)是最小的一個(gè)，即DH(I)是最接近精確解的.

第一輪計(jì)算完畢，進(jìn)行下一輪的逼近計(jì)算. 在DH(IP)附近取10個(gè)點(diǎn)組成新的DH(10). 根據(jù)DH(IP)值的特點(diǎn)，分3種情況計(jì)算：

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

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

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

這樣就定了新一輪計(jì)算的水平拉力DH(10). 值得注意的是，為保證每輪計(jì)算中都有10個(gè)元素，第一輪逼近中總共分割了10份，去掉為0的點(diǎn)，取上面10個(gè)點(diǎn)作為輸入. 而第二輪計(jì)算中由于第一個(gè)點(diǎn)不為0，所以分割了9份，取全部10點(diǎn)作為輸入.

如此反復(fù)進(jìn)行直到前后兩輪計(jì)算所得的DH(I)的差值滿足精度要求為止，至此求出了管道長(zhǎng)度為L(zhǎng)時(shí)管道下端水平拉力H. 至此就能求解非線性方程組(4)，就可以求出管道和纜索形態(tài)及內(nèi)力分布.

3仿真實(shí)例驗(yàn)證及其結(jié)果分析

3.1仿真作業(yè)算例

為了驗(yàn)證初始鋪管作業(yè)理論分析的準(zhǔn)確性、數(shù)值解算方法的快速性及穩(wěn)定性，鋪管作業(yè)仿真系統(tǒng)數(shù)學(xué)模型海況輸入采用第3方提供的DP數(shù)據(jù). 此次算例水深1 000 m，管道直徑323.9 mm，壁厚15.9 mm；纜繩長(zhǎng)度為1 000 m，纜繩外徑0.108 m，密度7 850 kg/m3；托管架長(zhǎng)度和曲率半徑分別為87 m和85 m. 為便于分析，提供鋪管過(guò)程中的6種狀態(tài)對(duì)應(yīng)的長(zhǎng)度及張力(見表1).

表1　起始鋪管管道長(zhǎng)度及張力

3.2仿真結(jié)果分析

將仿真結(jié)果與理論數(shù)據(jù)對(duì)比，驗(yàn)證1 000 m水深海底管道起始鋪設(shè)過(guò)程中管道和纜索的形態(tài)、軸力及彎矩分布. 仿真的管線形態(tài)如下圖5(a)～(f)，節(jié)點(diǎn)上部分為管道，下部分為纜索形態(tài).

圖5中實(shí)線代表仿真得到的曲線，虛線為商業(yè)軟件OFFPIPE計(jì)算得到的本次案例曲線，均為靜態(tài)管道鋪設(shè)形態(tài). 結(jié)果表明，隨著鋪管作業(yè)的進(jìn)行，管道長(zhǎng)度逐漸增加，管道下端高度逐漸降低. 從S1～S3可以看出，管道和纜索形態(tài)逐漸變平緩，纜索觸底點(diǎn)的位置逐步遠(yuǎn)離鋪管船，而在S3～S6中，纜索觸底點(diǎn)重新向鋪管船靠近，鋪管形態(tài)又變得陡峭.

在起始鋪管S1～S6中，管道上端的張力從180.0 kN逐漸增加到400.0 kN，鋪管作業(yè)過(guò)程中，起始鋪管軸力分布見圖6，起始鋪管彎矩分布見圖7.

由圖6、7可以看出，S1～S6過(guò)程中，纜索最大軸力逐漸變小，軸力變化趨勢(shì)逐漸變緩，S1時(shí)最陡，S6時(shí)最緩；纜索內(nèi)沒(méi)有彎矩，這是由于假設(shè)纜索的抗彎剛度為零.

(a) S1

(c) S3

(e) S5

(b) S2

(d) S4

(f) S6

圖6　起始鋪管軸力分布

圖7 起始鋪管彎矩分布

4結(jié)論

1)對(duì)初始鋪管作業(yè)過(guò)程中的管道與纜索進(jìn)行受力分析. 在原有的鋪管常用算法基礎(chǔ)上提出了將非線性彈性大撓度梁理論應(yīng)用于初始鋪管計(jì)算，建立管道數(shù)學(xué)模型，采用微分求積法求解，同時(shí)建立纜索模型，并用迭代法求解. 微分求積法極大地減少了計(jì)算量且求解精度較高，求解過(guò)程易于在計(jì)算機(jī)上程序?qū)崿F(xiàn)，提高實(shí)時(shí)仿真的真實(shí)度，有助于視景仿真系統(tǒng)的開發(fā).

2)在求解微分過(guò)程中，基于管道與纜索微分方程邊界條件難以確定導(dǎo)致方程無(wú)法求解的問(wèn)題，采用一種基于微分求積法的迭代方法，完成邊界條件的確立,然后用擬牛頓法完成方程組的求解. 利用已知的邊界條件迭代求出所需的邊界條件，在一定的誤差范圍內(nèi)結(jié)束迭代，結(jié)果表明，此方法有效解決了邊界條件問(wèn)題.

3)仿真得到的管道和纜索形態(tài)曲線與項(xiàng)目組提供的本次案例曲線的比較表明，仿真求解得到的曲線與理論曲線基本一致，表明作業(yè)理論分析的準(zhǔn)確性高,數(shù)值方法解算速度快, 穩(wěn)定性好.

參考文獻(xiàn)

[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] 顧永寧. 海底管線鋪管作業(yè)狀態(tài)分析[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] 花雷.曲梁結(jié)構(gòu)非線性大變形分析[D].揚(yáng)州:揚(yáng)州大學(xué),2011.

[10]李世榮,孫云,劉平.關(guān)于Euler-Bernoulli梁幾何非線性方程的討論[J].力學(xué)與實(shí)踐,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]王鑫偉. 微分求積法在結(jié)構(gòu)力學(xué)中的應(yīng)用[J]. 力學(xué)進(jìn)展, 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]黃玉盈, 朱達(dá)善. 海洋管線鋪設(shè)時(shí)的靜力分析[J]. 海洋工程, 1986, 4(1): 32-46.

[16]聶國(guó)雋, 仲政. 微分求積單元法在結(jié)構(gòu)工程中的應(yīng)用[J]. 力學(xué)季刊, 2005, 26(3): 423-427.

[17]馬成業(yè)，黃世華. 解非線性方程組的一個(gè)改進(jìn)牛頓法[J]. 甘肅聯(lián)合大學(xué)學(xué)報(bào)(自然科學(xué)版), 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

基金項(xiàng)目:國(guó)家科技重大專項(xiàng) (Z12SJENA0014)

作者簡(jiǎn)介:李震(1971—)，男，副教授，碩士生導(dǎo)師

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

中圖分類號(hào):TH133; TP183

文獻(xiàn)標(biāo)志碼:A

文章編號(hào):0367-6234(2016)07-0106-06