海管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