鋼懸鏈線立管初始位形及內(nèi)力分析?

李效民， 馬芳俊， 郭海燕

(中國海洋大學(xué)工程學(xué)院，山東 青島 266100)

鋼懸鏈線立管(SCR)出現(xiàn)于1990年代中期，由于成本相對較低以及對浮體升沉和漂移運(yùn)動較好的順應(yīng)能力等特點(diǎn)而成為深水油氣開發(fā)的首選立管形式[1]。隨著油氣資源勘測開發(fā)逐漸步入深海，鋼懸鏈線立管的設(shè)計、施工和安裝方法已經(jīng)成為油氣勘探開發(fā)行業(yè)重要的技術(shù)問題。

鋼懸鏈線立管(SCR)初始位形和初始內(nèi)力分析是對其進(jìn)行動力響應(yīng)分析的基礎(chǔ)，當(dāng)前國內(nèi)外學(xué)者對SCR的靜力位形的建模和求解等已進(jìn)行了大量研究。Ghadimi[2]采用集中質(zhì)量法建立立管模型，對立管進(jìn)行靜態(tài)分析和動態(tài)分析，并與實(shí)驗結(jié)果對比分析，表明該模型的合理性。O’Brien等[3]用類似懸鏈線索元確定了立管的平衡位形。Jain[4]用有限差分法對lazy-S型立管懸垂段部分進(jìn)行了靜力分析。Daniel等[5]基于懸鏈線位形附近小位移假設(shè)，利用線性理論導(dǎo)出了平衡方程，研究了在靜力荷載作用下鋼懸鏈線立管響應(yīng)的分析方法。Chai等[6]考慮了懸鏈線立管的彎矩、扭矩、剪力，提出了鋼懸鏈線立管的三維集中質(zhì)量公式。唐友剛等[7]采用集中質(zhì)量法研究了深海平臺系纜形狀和張力分析方法，指出需要考慮海底地形對系纜形狀和張力的影響。Santillan等[8]基于幾何-彎矩-曲率關(guān)系建立了立管平衡方程并采用有限差分法對陡峭型S立管和惰性S立管靜力位形進(jìn)行了數(shù)值模擬和實(shí)驗對比。白興蘭等[9]基于大撓度柔性索理論，采用具有彎曲剛度的細(xì)長梁模型模擬鋼懸鏈線立管。李艷等[10]采用三維集中質(zhì)量法，建立了立管構(gòu)型及張力的計算模型。

綜上所述，當(dāng)前鋼懸鏈線立管初始位形確定的方法主要有：根據(jù)懸鏈線方程建立鋼懸鏈線立管模型，忽略管道的彎曲剛度，只考慮自重和浮力，用有限元靜力分析迭代求解獲得立管的初始位形；用集中質(zhì)量法來模擬鋼懸鏈線立管，通過在節(jié)點(diǎn)之間增加伸長彈簧和旋轉(zhuǎn)彈簧模擬立管的軸向及彎曲變形；根據(jù)平衡原理對微元體列控制方程，采用有限元法或有限差分法求解方程等?？傊蠖际峭ㄟ^各種結(jié)構(gòu)動力學(xué)方法建立立管控制方程，然后基于有限元方法或有限差分法等對控制方程進(jìn)行數(shù)值求解，控制方程和數(shù)值計算是彼此獨(dú)立分離的概念。而對于復(fù)雜的立管幾何構(gòu)型來說控制方程的建立是一個及其復(fù)雜的過程，而且一旦邊界條件改變，有可能需要重新建立計算模型。而基于向量式有限元(VFIFE)的分析方法是與上述方法截然不同的一種新型建模和數(shù)值計算方法。

由丁承先[11]教授提出的向量式有限元，又稱有限質(zhì)點(diǎn)法，以向量力學(xué)作為分析的理論基礎(chǔ)，并與數(shù)值計算結(jié)合，是求解結(jié)構(gòu)的大變形、大變位、彈塑性、碰撞、倒塌等非線性或不連續(xù)性力學(xué)行為的新方法[12]。向量式有限元基于結(jié)構(gòu)的物理模式，其計算步驟可概括為：將結(jié)構(gòu)離散為滿足牛頓第二定律的質(zhì)點(diǎn)；時間歷程被分為一系列的途徑單元；在途經(jīng)單元內(nèi)，采用虛擬的逆向運(yùn)動，通過單元的純變形，計算結(jié)構(gòu)的內(nèi)力；采用中央差分的顯示積分法求解運(yùn)動公式。

本文基于向量式有限元提出了一種新的確定鋼懸鏈線立管初始位形求解方案?；谠摪阜桨笇⒘⒐茈x散成一系列質(zhì)點(diǎn)的集合，質(zhì)點(diǎn)間用只承受內(nèi)力而無質(zhì)量的彎曲桿件元連接，采用牛頓第二定律來表達(dá)質(zhì)點(diǎn)的運(yùn)動行為，并用中央差分的顯式積分法求解質(zhì)點(diǎn)各時刻的控制方程，編制相應(yīng)Matlab求解程序。分析了三種不同立管形式，并與他人研究結(jié)果進(jìn)行對比分析，結(jié)果表明基于向量式有限元的位形求解方案是可行和正確的，最后指出了向量式有限元法相較于傳統(tǒng)有限元法存在的優(yōu)勢。

1 鋼懸鏈線立管初始位形分析模型

1.1 初始位形建模方案

要對立管展開靜力和動力分析，首先要確定立管的初始位形，以初始形態(tài)作為計算位移、變形和應(yīng)力等的基礎(chǔ)架構(gòu)。對于懸鏈線立管初始位形的確定可分為兩類問題。一是初值問題：已知某一端的所有邊界條件，包括該端的位置坐標(biāo)、張力的水平分量和豎直分量。二是邊值問題：即兩端的位置坐標(biāo)均已知，而其兩端的張力分量均未知。本文采用如圖1所示方案解決立管的這兩類問題。

圖1 向量式有限元法建模流程圖Fig.1 Flow scheme for VFIFE

1.2 立管模型建立

基于上述方案，假設(shè)立管在初始位置時是以一定角度傾斜的直管，有質(zhì)量且暫不考慮重力和浮力作用，如圖1所示。確定立管整體坐標(biāo)系，坐標(biāo)系定義如下：取立管底部為坐標(biāo)原點(diǎn)，x軸與來流方向平行，y軸為水深方向，向上為正，z軸垂直于x軸和y軸。

圖2 立管的變形及坐標(biāo)示意圖Fig.2 Diagram of the riser deformation and coordinate

采用N+1(從底向上編號依次為1，2，…，N，N+1)個等間距的質(zhì)點(diǎn)模擬立管的直線狀態(tài)，質(zhì)點(diǎn)之間利用平面彎曲桿單元相連接，單元的質(zhì)量平均分配到兩端的質(zhì)點(diǎn)上，立管的初始張力通過單元的初始軸力來實(shí)現(xiàn)。立管底端鉸接，將N+1號質(zhì)點(diǎn)在rtime時間內(nèi)(假設(shè)位置到指定位置的模擬時間)水平運(yùn)動到指定位置L1點(diǎn)，同時在rtime時間內(nèi)將重力和浮力以斜坡函數(shù)加載的方式作用在各質(zhì)點(diǎn)上，立管在重力、浮力與單元作用力下最終達(dá)到平衡狀態(tài)，即可得到如圖2(a)所示的鋼懸鏈線立管的初始位形。

1.3 單元節(jié)點(diǎn)內(nèi)力計算

(1)

(2)

(3)

圖3 單元的逆向運(yùn)動Fig.3 The retrogression of element

(4)

(5)

將兩個節(jié)點(diǎn)內(nèi)力分量分別轉(zhuǎn)換成域坐標(biāo)分量，

(6)

最后進(jìn)行內(nèi)力集成。

1.4 質(zhì)點(diǎn)控制方程

對于空間桿系結(jié)構(gòu)，質(zhì)點(diǎn)有3個位移分量，考慮阻尼的影響，質(zhì)點(diǎn)控制方程式為：

(7)

1.5 控制方程式的求解

向量式有限元求解時采用中央差分的顯示時間積分法，這是為了避免隱式解法帶來的復(fù)雜迭代和收斂問題。但為了計算的準(zhǔn)確性，必須控制時間步長h。

(8)

其中C1=1/(1+ζh/2)；C2=C1(1-ζh/2)。

計算步驟可總結(jié)為：

之后，用n步代替n-1步，用n+1步代替n步，更新位置向量，重復(fù)上述步驟，循環(huán)計算，直到完成分析。

2 算例分析

基于上述理論及計算方法，本文編制相應(yīng)的Matlab求解程序，對鋼懸鏈立管初始位形及內(nèi)力進(jìn)行分析，為驗證程序的可靠性，本文選取三種立管進(jìn)行計算。

選取的立管一為高秦嶺[13]所采用的簡單鋼懸鏈線立管，具體參數(shù)為：水深1 000 m，長1 298 m，外徑0.274 m，壁厚0.01 m，海水密度為1 025 kg·m-3，內(nèi)部流體密度為1 000 kg·m-3，立管管材密度7 850 kg·m-3，彈性模量2.07 e11N·m-2；本程序中rtime=450.0 s。

圖4繪制了立管一在初始時刻、rtime/2時刻、rtime時刻以及最終時刻的位形圖，顯示了立管一從假設(shè)的虛擬斜直狀態(tài)在重力和浮力等載荷作用下達(dá)到最終尋找的初始位形的形成過程。圖5是分別通過向量式有限元、經(jīng)典懸鏈線方程以及Ansys軟件分析計算出的立管一整體位形，圖6是立管一在觸地點(diǎn)附近的局部位形圖。由圖5可看出向量式有限元法得出的整體位形與其他方法得出的結(jié)果相當(dāng)吻合。從圖6可以看出三種方法得到的初始位形在觸地點(diǎn)局部略有差異，原因在于懸鏈線法忽略了彎曲剛度的影響，文獻(xiàn)[13]中通過設(shè)定位形判別條件，一步步迭代，使立管最終達(dá)到平衡狀態(tài)即SCR立管的初始位形，而采用向量式有限元計算時，結(jié)構(gòu)荷載可以在計算分析的初始步全部加在立管上，立管在力的作用下形成初始位形，向量式有限元是對立管的力學(xué)行為建立的模型，位形的形成中不需要迭代求解，所得結(jié)果相對精確。

圖4 不同時段立管一位形圖Fig.4 The first riser`s configuration at different times

圖5 三種模型得到的立管一整體位形比較Fig.5 The configuration`s comparison of three calculation models about the first riser

本文立管二采用Chen[14]、Ghadimi[2]的長度為540 m的簡單懸鏈線立管，具體參數(shù)為：水平投影距離255 m，豎直投影距離375 m，外徑0.276 6 m，海水密度為1 025 kg·m-3，立管干重102 kg·m-2，軸向剛度3.27e8N，抗彎剛度3.4e4N·m2；本程序中rtime=550.0 s。

通過向量式有限元法、細(xì)長柔性桿模型[14]、理論計算和集中質(zhì)量法[2]得到頂部張力分別為：163.535、166.6、165和167 kN，誤差不大于2.07%。得到的拖地段長度分別為：120、123.5、121.5和120 m，誤差不大于2.83%。

圖6 三種模型在立管一觸地點(diǎn)附近位形Fig.6 The configurations near TDP with three calculation models about the first riser

圖7 不同時段立管二位形Fig.7 The second riser′s configuration at different times

如圖7，顯示了立管二從假設(shè)的虛擬斜直立管在重力和浮力等載荷作用下達(dá)到初始位形的過程。圖8是分別通過向量式有限元法、細(xì)長柔性管模型[14]、理論計算和集中質(zhì)量法[2]得到簡單懸鏈線立管的整體位形，由圖可知四種方法得到的整體位形基本吻合。圖9是觸地點(diǎn)附近局部位形，從圖中看出四種方法計算結(jié)果略有差異，究其原因還是由于計算模型的差異所導(dǎo)致的，集中質(zhì)量法采用伸長彈簧和旋轉(zhuǎn)彈簧分別模擬立管的軸向及彎曲變形，而向量式有限元法采用質(zhì)點(diǎn)及質(zhì)點(diǎn)間的平面彎曲桿件模擬，兩者不同導(dǎo)致位形差異。

圖8 四種模型得到的立管二整體位形比較Fig.8 The configuration`s comparison of four calculation models about the second riser

圖9 四種模型在立管二觸地點(diǎn)附近位形Fig.9 The configurations near TDP with four calculation models about the second riser

圖10繪制了540 m管道張力分布圖，與姬鸞[15]的結(jié)果比較可以看出兩者差異很小。集中質(zhì)量法是通過桿單元模擬單元節(jié)點(diǎn)，用旋轉(zhuǎn)彈簧-阻尼系統(tǒng)等效彎曲剛度，并需要建立方程，通過迭代求解；VFIFE使用的是梁單元，用質(zhì)點(diǎn)軌跡來模擬整體運(yùn)動行為，無需迭代求解。兩者之間計算方式的不同造成計算結(jié)果有所差異。VFIFE是基于物理模型，模擬立管的力學(xué)行為，更符合實(shí)際受力情況。

本文選取的立管三為Wang[16]所采用的緩波形鋼懸鏈線立管，具體參數(shù)為：水深1 850 m，立管長度2 600 m，懸垂段1 690 m，浮子段520 m，下降段390 m，外徑0.203 2 m，壁厚0.019 1 m，海水密度為1 024 kg·m-3, 內(nèi)部流體密度為998 kg·m-3，立管管材密度7 860 kg·m-3，彈性模量2.06e11N·m-2；本程序中rtime=2200.0 s。

Wang[16]給出的數(shù)值計算和OrcaFlex計算得到的頂端張力分別是820.75、820.94 kN，向量式有限元法得到的頂端張力為819.216 kN，誤差不大于0.21%。Wang[16]兩種方法得到觸地點(diǎn)處張力分別是75、75.056 kN，向量式有限元法得到的結(jié)果是75.159 kN，相比誤差不大于0.212%。Wang[16]兩種計算方法得到的最大彎矩分別是67.441、67.074 kN·m，本文得到的最大彎矩為68.215 kN·m，誤差不大于1.70%。

圖10 兩種計算模型得到的立管二張力比較Fig.10 The tension comparison of two calculation models about the second riser

圖11 不同時段立管三位形Fig.11 The third riser′s configuration at different times

圖12 四種模型得到的立管三整體位形比較Fig.12 The configuration′s comparison of four models about the third riser

圖13 兩種計算模型得到的立管三張力比較Fig.13 The tension comparison of two calculation models about the third riser

圖14 兩種計算模型得到的立管三彎矩比較Fig.14 The bending moment comparison of two calculation models about the third riser

圖11顯示的是立管三從假設(shè)的虛擬斜直立管在重力和浮力等載荷作用下達(dá)到初始位形的過程。圖12是通過向量式有限元法計算出的位形與Wang[16]給出的的位形比較圖，由圖12可看出向量式有限元法得出的位形圖與其他方法得出的位形基本吻合。圖13和圖14分別是采用本文方法得到的緩波形立管的張力和彎矩，與文獻(xiàn)對比圖。從圖中可以看出結(jié)果相當(dāng)吻合，證明了向量式有限元法分析緩波形鋼懸鏈線立管的可行性及正確性。

3 結(jié)論與展望

本文提出了一種確定鋼懸鏈線立管初始位形新的求解方案，基于向量式有限元對鋼懸鏈線立管的整體位形和初始內(nèi)力進(jìn)行分析，對比傳統(tǒng)有限元分析過程和結(jié)果，可以看出：

(1)傳統(tǒng)有限元方法基于變分原理，需要由單元剛度矩陣集成為整體剛度矩陣，計算過程中不容易增減單元和改變邊界條件。向量式有限元采用的是物理模式，對立管的力學(xué)行為建立計算理論，計算中不需集成結(jié)構(gòu)的剛度矩陣，可方便更改單元和邊界條件，也不需迭代求解非線性方程組，從而簡化了立管分析的難度。

(2)向量式有限元對立管的分析即使是初始靜力位形分析都是動態(tài)的分析。不同于傳統(tǒng)有限元方法的靜態(tài)—動態(tài)分開分析，向量式有限元分析初始時可以將荷載直接加載到立管上，而且分析過程中可以任意改變荷載，可以更好的模擬真實(shí)海洋環(huán)境下立管的運(yùn)動，使分析更加準(zhǔn)確。

(3)向量式有限元法計算程序更加簡潔。向量式有限元法計算程序只存在質(zhì)點(diǎn)位置計算和單元內(nèi)力計算兩個循環(huán)，整體步驟呈系統(tǒng)化；采用中央差分的顯式時間積分求解，避免了迭代和收斂問題，從而使分析更加簡潔。

結(jié)果表明，基于向量式有限元法并采用本文提出的位形求解方案是完全可行的，為懸鏈線立管初始位形和內(nèi)力計算提供了一種條理清晰和高效的數(shù)學(xué)模型及可靠穩(wěn)定的數(shù)值模擬方法?？傊蛄渴接邢拊ㄔ诜治鰬益溇€立管時具有很大的優(yōu)勢，為進(jìn)一步對各種形式的海洋立管復(fù)雜行為尤其是動力分析提供了一種新的可行的方法，其結(jié)果亦可與其他求解方案對比分析以此保證立管服役期間的安全可靠運(yùn)行。