海洋管道單點(diǎn)提放的若干數(shù)學(xué)模型研究

曾霞光，段夢(mèng)蘭，陳景皓

(中國(guó)石油大學(xué)(北京)機(jī)械與儲(chǔ)運(yùn)工程學(xué)院，北京 102249)

隨著中國(guó)深海油氣開發(fā)的進(jìn)程，中國(guó)必須掌握深水管道鋪設(shè)和安裝的分析設(shè)計(jì)技術(shù)。深水海洋油氣開發(fā)中，主要的鋪管安裝方式有S型鋪設(shè)，J型鋪設(shè)等。鋪管和安裝作業(yè)中，回收和棄置管道操作非常普遍。在這些作業(yè)過(guò)程中，管道的軸力，彎矩，應(yīng)力，應(yīng)變等需要控制在一定的范圍內(nèi)以保證管道完好，因此海洋管道鋪設(shè)和收棄作業(yè)方案，尤其是船舶位置移動(dòng)設(shè)計(jì)和纜繩的收放速度非常重要。為了進(jìn)行設(shè)計(jì)分析，需要建立海洋管道的收棄作業(yè)過(guò)程分析模型和方法。

鋪管作業(yè)可以看作棄管作業(yè)的前期工作或者收管作業(yè)的后期工作，因此用于分析的力學(xué)模型與收棄管分析模型相同。棄置時(shí)海洋管道通常將脫離托管架或其他支撐結(jié)構(gòu)，A＆R絞車的鋼繩通過(guò)海洋管道封頭將管道從海面緩慢下放到海底，回收時(shí)A＆R絞車的鋼繩通過(guò)海洋管道封頭將管道從海底緩慢拉升到海面上，此過(guò)程被稱為海洋管道的單點(diǎn)提放。在這個(gè)往復(fù)過(guò)程中，管道經(jīng)歷了一個(gè)大變形和小變形相互轉(zhuǎn)化的變形過(guò)程。

1 控制微分方程

海洋管道從海底提升到海面或者從海面下放到海底，這兩個(gè)物理過(guò)程是互逆的，可以用同一個(gè)數(shù)學(xué)模型進(jìn)行描述。在管道提放作業(yè)時(shí)通常要求海洋環(huán)境平穩(wěn)，管道提升或者下放速度緩慢，因此可把問(wèn)題處理成靜態(tài)問(wèn)題。另外可認(rèn)為海底為剛性水平面［1］，海洋管道由各向同性材料組成，且提升和下放過(guò)程中始終處于彈性狀態(tài)。在上述條件下，不考慮海水流動(dòng)效應(yīng)，取懸空管段進(jìn)行研究，如圖1所示。以懸空管道與海床的接觸點(diǎn)(觸地點(diǎn))為原點(diǎn)o，建立笛卡爾坐標(biāo)系。再沿管線建立自然坐標(biāo)系，顯然該管線上的物理量均為管線弧長(zhǎng)的單變量函數(shù)。令管道與水平方向的夾角為θ，長(zhǎng)度為L(zhǎng)，材料彈性模量為E，截面慣性矩為I，水中的單位線重為w，上端點(diǎn)受拉力T0作用，與水平方向成夾角θ0，管線原點(diǎn)處橫截面的水平力為H0，鉛垂力為V0。其中鉛垂力V0的可能來(lái)源有且僅有兩個(gè)，一個(gè)是左邊管道對(duì)其剪力，另一個(gè)是剛性接觸海床對(duì)其向上的支撐力。由于通常認(rèn)為左邊管道沒(méi)有彎曲變形，其彎矩為零，因此其相互剪力也為零，否則左邊管道必然彎曲。由此可知原點(diǎn)處的鉛垂力就是剛性海床對(duì)其向上的支撐力。

根據(jù)文獻(xiàn)［2-5］，如圖1所示的海管提放問(wèn)題可以用如下微分方程描述:

圖1 管道單點(diǎn)提放力學(xué)模型Fig.1 Mechanical model of pipe lifting or lowering by one point

根據(jù)水平和垂直方向上力的平衡有:

將式(2)、(3)代入方程(1)得:

根據(jù)文獻(xiàn)［6］，該問(wèn)題也可以用以下控制微分方程組描述:

用這兩種微分方程對(duì)管道收放問(wèn)題進(jìn)行分析，結(jié)果差異何在?這兩種控制微分方程是否都能很好地描述海洋管道的提放過(guò)程?哪種微分方程更好地描述海洋管道的提放過(guò)程?為了回答上述及相關(guān)問(wèn)題，將以兩種微分方程為核心建立海洋管道提放過(guò)程的三個(gè)分析模型，并借助數(shù)值計(jì)算方法分別進(jìn)行求解，然后對(duì)計(jì)算結(jié)果進(jìn)行分析。

2 邊界條件

分析海洋管道提放過(guò)程中邊界處的形態(tài)和受力情況，通常認(rèn)為有如下邊界條件［3，5-7］:

觸地點(diǎn)處:撓度 y(0)=0，轉(zhuǎn)角 θ(0)=0，水平力 H0=T0cosθ0，彎矩 M(0)=0。

上部端點(diǎn)處:彎矩M(L)=0，拉力為T0。

3 懸空管道長(zhǎng)度

懸空管道長(zhǎng)度L是求解本問(wèn)題的關(guān)鍵參數(shù)，不能預(yù)先確定。文獻(xiàn)［1］認(rèn)為剛性海床情況下提起段長(zhǎng)度等于2倍的提升力除以管道線重，文獻(xiàn)［2-3］根據(jù)幾何條件確定懸空管道長(zhǎng)度，而文獻(xiàn)［5-8］都用各自的迭代方法確定該參數(shù)。與這些方法不同，將根據(jù)力學(xué)原理提供一個(gè)新的確定懸空管道長(zhǎng)度的公式。

考慮懸空管道的外力對(duì)管道上部端點(diǎn)的彎矩，根據(jù)彎矩平衡原理有:

將式(2)、(3)代入式(6)得:

式(7)中只有一個(gè)獨(dú)立未知量變量L，L一旦確定，x(L)，y(L)和x(S)就隨之確定，因此可以求解該方程來(lái)確定懸空管道長(zhǎng)度。求解方法是迭代法，思路為不斷假設(shè)一個(gè)懸空管道長(zhǎng)度L，求解以式(4)或式(5)為控制微分方程的定解問(wèn)題得到管道形態(tài)，直到下式成立:

此時(shí)即可認(rèn)為該長(zhǎng)度就是管道懸空長(zhǎng)度，其中ε為迭代計(jì)算誤差限。

4 數(shù)學(xué)模型

取微分方程(4)，選擇如下邊界條件:觸地點(diǎn)處的轉(zhuǎn)角為零，上部端點(diǎn)處的彎矩為零，再結(jié)合式(7)，得到分析模型1，即:

取微分方程(4)，選擇如下邊界條件:觸地點(diǎn)和上部端點(diǎn)處的彎矩都為零，再結(jié)合式(7)，得到分析模型2，即:

取微分方程組(5)，選擇如下邊界條件:觸地點(diǎn)處的轉(zhuǎn)角為零，彎矩為零，水平力為T0cosθ0，上部端點(diǎn)處的彎矩為零，再結(jié)合式(7)，得到分析模型3，即:

5 數(shù)值求解算例及結(jié)果分析

模型1和模型2的未知函數(shù)是θ(S)，模型3的未知函數(shù)是θ(S)和T(S)，由于上述微分方程均是高度非線性方程，無(wú)法得到其解析形式的解。為了分析這些模型的差異，并找出準(zhǔn)確描述海洋管道提放物理過(guò)程的數(shù)學(xué)模型，下面對(duì)它們進(jìn)行數(shù)值求解。

采用MATLAB作為求解環(huán)境，利用其中的內(nèi)置函數(shù)bvp4c［9］求解如表1、2所示的一組具體參數(shù)對(duì)應(yīng)的分析模型1、2和3，計(jì)算結(jié)果如圖2所示。

圖2 不同拉力作用時(shí)管道的轉(zhuǎn)角、形態(tài)和彎矩Fig.2 Angle，configuration and bending moment of pipe loading under different tension

從計(jì)算結(jié)果圖2中可以發(fā)現(xiàn)以下主要特點(diǎn):

1)模型3確定的懸空管道長(zhǎng)度在三個(gè)模型中是最大的，模型2確定的懸空管道長(zhǎng)度最小;

2)模型1的彎矩極值在三個(gè)模型中是最大的，發(fā)生在觸地點(diǎn)處，模型2的彎矩極值總是最小，模型2和3的彎矩極值都不發(fā)生在觸地點(diǎn)處;

3)彎矩極值隨著懸空管道長(zhǎng)度的變化而變化，其先隨懸空管道長(zhǎng)度增加而增加，之后隨懸空管道長(zhǎng)度增加而減小，最大值發(fā)生在懸空管道約100 m處;

4)隨著載荷增大三個(gè)模型的解越來(lái)越接近，在懸空管道長(zhǎng)度小于1 000 m時(shí)，三個(gè)模型的轉(zhuǎn)角、形態(tài)和彎矩的結(jié)果差別較大，且長(zhǎng)度越小差別越大。在懸空管道長(zhǎng)度大于1 000 m時(shí)，模型2和模型3的轉(zhuǎn)角、形態(tài)和彎矩的計(jì)算結(jié)果幾乎相同，但是模型1的結(jié)果與其它兩個(gè)模型在彎矩的極值及其發(fā)生位置存在較大差異;

5)從圖2(a)～(d)中可以清楚看出，模型1的彎矩極值發(fā)生在觸地點(diǎn)，這與實(shí)際情況矛盾，而模型2在觸地點(diǎn)的轉(zhuǎn)角值也與該問(wèn)題的邊界條件矛盾。

表1 材料幾何參數(shù)Tab.1 Geometrical and material parameters

表2 載荷參數(shù)Tab.2 Loading parameters

根據(jù)計(jì)算結(jié)果的特點(diǎn)，可以初步判定描述海洋管道提放全程的數(shù)學(xué)模型應(yīng)該采用分析模型3。為了進(jìn)一步證實(shí)該模型的正確性，作者以分析模型3為理論模型，編寫開發(fā)了海洋管道提放靜力學(xué)分析軟件DRICAS，采用上文中的管道材料幾何參數(shù)，重新計(jì)算了如表3所示載荷情況下的管道形態(tài)。

表3 載荷參數(shù)Tab.3 Loading parameters

將其形態(tài)計(jì)算結(jié)果與專業(yè)海洋工程有限元分析軟件OrcaFlex的計(jì)算結(jié)果進(jìn)行對(duì)比，可以看到兩者的計(jì)算結(jié)果在各個(gè)載荷情況下都是一致的，如圖3所示。這進(jìn)一步證實(shí)了模型3的正確性。

圖3 DRICAS和OcraFlex形態(tài)結(jié)果對(duì)比Fig.3 Configuration results comparison between DRICA and OcraFlex

6 結(jié)語(yǔ)

通過(guò)三個(gè)模型的數(shù)值求解和分析表明:

1)對(duì)海洋管道緩慢單點(diǎn)提放進(jìn)行全程分析宜采用模型3，它能準(zhǔn)確描述管道提放的完整過(guò)程;其中該問(wèn)題的邊界條件是觸地點(diǎn)處水平轉(zhuǎn)角為零，彎矩為零，受到提升載荷水平分量大小的拉力，上部端點(diǎn)處彎矩為零;另外懸空管道長(zhǎng)度可以根據(jù)彎矩平衡原理確定。

2)因?yàn)閺澗貥O值位置不合理，模型1不能準(zhǔn)確描述海洋管道單點(diǎn)提放過(guò)程并用于管道的力學(xué)分析。

3)懸空管道長(zhǎng)度小于900 m時(shí)，模型2不能準(zhǔn)確描述海洋管道單點(diǎn)提放過(guò)程;懸空管道長(zhǎng)度大于1 000 m時(shí)，模型2與模型3計(jì)算結(jié)果的相對(duì)誤差很小，且懸空管道越長(zhǎng)相對(duì)誤差越小，因此在懸空管道長(zhǎng)度大于1 000 m時(shí)可用模型2代替模型3，從而減少計(jì)算機(jī)數(shù)值計(jì)算時(shí)間。

4)彎矩極值隨著懸空管道長(zhǎng)度變化而變化，其先隨著懸空管道長(zhǎng)度增加而增加，之后隨著懸空管道長(zhǎng)度增加而減小，最大值往往發(fā)生在管道提升過(guò)程的中前期或者下放過(guò)程的中后期。

［1］ 邢靜忠.海底管道離地點(diǎn)的接觸條件探討［J］.力學(xué)與實(shí)踐，2005，27(5):57-60.

［2］ 宋甲宗，戴英杰.海洋管道鋪設(shè)時(shí)的二維靜力分析［J］.大連理工大學(xué)學(xué)報(bào)，1999，39(1):91-94.

［3］ 戴英杰，宋甲宗，馮 剛.海底管道收棄管作業(yè)分析［J］.海洋工程，2000，18(3):75-78.

［4］ 邢靜忠，柳春圖，曾曉輝.海洋石油管道單點(diǎn)提升分析［J］.海洋工程，2002，20(3):29-33.

［5］ 曾曉輝，柳春圖，邢靜忠.海底管道鋪設(shè)的力學(xué)分析［J］.力學(xué)與實(shí)踐，2000，24(2):19-21.

［6］ Charles P.Sparks.Fundamentals of Marine Riser Mechanics［M］.Tulsa，Oklahoma:PennWell，2007:254-256.

［7］ 詹 侃，陸仁華.海底管道在鋪設(shè)過(guò)程中的二維靜態(tài)分析［J］.海洋工程，1991，9(4):14-20.

［8］ 帥 健，呂英民，張進(jìn)國(guó).用動(dòng)坐標(biāo)迭代法分析海洋管道敷設(shè)時(shí)的應(yīng)力［J］.油氣儲(chǔ)運(yùn)，1996，15(10):21-24.

［9］ Shampine L F，Gladwell I，Thompson S.Solving ODEs with MATLAB［M］.New York:Cambridge University Press，2003:133-211.