李效民， 張　林， 牛建杰， 韓永慶， 郭海燕

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

基于向量式有限元的深水頂張力立管動力響應分析

李效民， 張林， 牛建杰， 韓永慶， 郭海燕

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

向量式有限元法亦稱有限質(zhì)點法是以向量式力學與數(shù)值計算為基礎(chǔ)的一種結(jié)構(gòu)行為分析的創(chuàng)新性方法。該方法將結(jié)構(gòu)離散為一系列質(zhì)點的集合，質(zhì)點間采用單元連接，單元沒有質(zhì)量，只承受內(nèi)力，質(zhì)點的運動則直接采用牛頓第二定律來表達，采用中央差分的顯式積分法求解質(zhì)點群各時刻運動方程。首次將該方法應用于頂張力立管的動力行為分析中，采用平面彎曲桿件元模擬質(zhì)點間的相互作用力，推導了立管質(zhì)點運動公式，編制相應Matlab求解程序，計算了立管在海流、波浪作用下的動力響應特性并與基于傳統(tǒng)有限元法的結(jié)果進行了對比分析，表明該方法在立管行為分析中的正確和有效性，并指出了該方法在處理其它立管形式行為分析中的優(yōu)勢。

向量式有限元(VFIFE)；有限元法(FEM)；向量式力學；頂張力立管(TTR)

頂張力立管(TTR)是各種形式的海洋立管中既可用于鉆井、完井及修井又可以用于注水、生產(chǎn)等功能的立管，是油氣勘探開發(fā)中必不可少的立管形式。近幾十年來受到了國內(nèi)外科技界和工業(yè)界的廣泛關(guān)注，對其展開了深入研究?？傮w來說，研究分析主要集中在渦激振動、參數(shù)激勵(穩(wěn)定性)、順流向振動響應研究等幾方面，并采用理論分析、試驗研究與數(shù)值模擬相結(jié)合的方式展開相關(guān)研究。

Morkookaza等[1]在時域和頻域內(nèi)對比分析了頂張力立管在不同波浪力和海洋環(huán)境條件下的動力響應。Chatjigeorgiou等[2]研究了垂直海洋立管在其頂端承受參數(shù)激勵下的非線性響應特性，詳細分析了水動力阻尼對系統(tǒng)動力不穩(wěn)定區(qū)域的影響。Kuiper等[3]研究了深水立管參數(shù)激勵的穩(wěn)定性，表明浮式結(jié)構(gòu)的升沉運動能夠?qū)е铝⒐艹霈F(xiàn)不穩(wěn)定性。Chaplin等[4]為增加立管渦激振動物理特性的理解并得到詳實的試驗數(shù)據(jù)來驗證渦激振動預測模型，在Delta 水槽中設(shè)計了一精巧的頂張力立管渦激振動實驗研究方案。高云等[5]比較系統(tǒng)的研究了不同鰭高、不同螺距的Strake下立管狀態(tài)下的響應特性。王騰等[6]模擬了一鉆井隔水管在波、流和平臺運動聯(lián)合作用下的非線性動力響應，研究了平臺的靜偏移和動態(tài)運動對隔水管最大彎矩分布的影響。暢元江等[7]考慮波流聯(lián)合作用下隔水管與導向架之間的間隙接觸條件，用ANSYS研究了非線性邊界條件下隔水管的動態(tài)響應。楊和振等[8]采用Floquet理論分析了深水立管參數(shù)穩(wěn)定性，指出了參數(shù)共振區(qū)域，分析了阻尼對參數(shù)激勵不穩(wěn)定性的影響。唐友剛等[9]研究了頂張力立管在剪切流中的渦激振動響應和參激-渦激耦合振動響應特性，分析參數(shù)激勵對立管橫向渦激振動的影響。陳偉民等[10]采用有限元法對浮體升沉與立管渦激振動的耦合動力響應進行了分析。唐國強等[11]建立了一種立管渦激振動的經(jīng)驗模型，研究了立管頂部張力、水流速度、立管外徑以及立管管壁厚度的變化對于立管動力響應的影響。

綜上所述，當前人們對頂張力立管的研究大都是采用基于小變形小位移理論的梁單元的傳統(tǒng)有限元方法進行模擬。在立管的運動響應計算過程中剛度矩陣不隨時間變化，也即始終以其初始垂直位形作為立管的計算剛度矩陣。對于深水立管來說，當其相對位移較小時該方法是足夠精確的，然而當立管在海流和頂端較大浮體偏移下的運動響應計算誤差則相對偏大，因為嚴格意義來說，此時立管運動行為屬于小變形大變位的幾何非線性問題。對于立管的大變位問題可以采用在應變計算中包含高階項、推導平衡方程時考慮變形的影響、Lagrange(T.L.)法、修正的Lagrange(U.L.)法和動坐標法/共旋坐標法等諸多方法來求解，不過對其進行建模求解需要特殊技巧和方案，同時計算中難以區(qū)分立管單元的剛體運動和純變形，此時對立管動態(tài)響應求解將變得及其復雜。

向量式有限元法(VFIFE)亦被稱作有限質(zhì)點法(FPM)是由丁承先等[12]基于向量力學與數(shù)值計算相結(jié)合提出的一種新型結(jié)構(gòu)行為分析方法，是求解結(jié)構(gòu)大變形、大變位、彈塑性、碰撞等非線性或不連續(xù)性力學行為的新方法。與傳統(tǒng)有限元需要建立完整的數(shù)學理論描述不同，向量式有限元法以結(jié)構(gòu)的物理模式為基礎(chǔ)，將結(jié)構(gòu)離散為質(zhì)點群，質(zhì)點間采用單元連接，單元沒有質(zhì)量，只承受內(nèi)力，利用牛頓第二定律描述這些質(zhì)點的運動，采用途徑單元來描述質(zhì)點的運動，采用虛擬的逆向運動計算結(jié)構(gòu)的純變形，可以同時計算結(jié)構(gòu)的運動和變形，而不采用傳統(tǒng)有限元方法使用的連體控制方程描述結(jié)構(gòu)。倪秋斌等[13]基于向量式有限元法考慮拉索幾何和材料非線性對斜拉橋振動控制進行了數(shù)值模擬。喻瑩等[14-15]基于有限質(zhì)點法對結(jié)構(gòu)的倒塌破壞過程進行了模擬和分析。羅堯治等[16]闡述了有限質(zhì)點法這種新的數(shù)值分析方法在空間結(jié)構(gòu)復雜行為研究領(lǐng)域的優(yōu)勢及應用，并對該方法的發(fā)展趨勢作出了展望。班自愿等[17]對懸臂梁分別承受集中載荷和彎矩下的大變形進行算例分析。袁峰[18]采用彈性海床和塑性海床模擬管土相互作用，將向量式有限元法初步應用于海洋管土動力循環(huán)中，取得了一些有益的分析結(jié)果。

本文首次將向量式有限元法應用于深水頂張力立管的動力響應研究中。將立管離散為一群質(zhì)點的集合，采用平面彎曲桿件元模擬質(zhì)點間的相互作用力，建立了立管質(zhì)點運動計算公式，編制相應Matlab求解程序，計算了深水立管在海流、波浪作用下的動力響應特性并與基于傳統(tǒng)有限元法(FEM)的結(jié)果進行了對比分析。通過對比分析，驗證向量式有限元法在深水頂張力立管動力行為分析中的有效性及其優(yōu)勢，并為進一步將這種新型數(shù)值計算方法在深水各種形式的海洋立管平面內(nèi)外復雜環(huán)境載荷條件下的行為分析中的應用前景作出展望。

1立管質(zhì)點運動方程的建立及載荷

1.1立管質(zhì)點的點值描述

立管整體坐標系及其振動變形及坐標示意圖如圖1所示。假設(shè)頂張力立管在初始位置時完全垂直。設(shè)定一組域坐標(x，y，z)，基底為(i，j，k)，如圖1(a)所示。坐標系定義如下：以立管未發(fā)生變形時垂直軸為y軸，向上為正，并取立管底部為原點坐標；x軸與來流方向平行，z軸垂直于x和y軸。 將立管分為N段共N+1(從底向上編號依次為1, 2，… ，I,J, …,N,N+1)個空間點來進行描述。等效節(jié)點質(zhì)量用tn-1時立管的密度、橫截面積、跨度以及內(nèi)部流體的密度、斷面面積、計算跨度確定。初始時刻各點沿x,y向平移速度、轉(zhuǎn)角和旋轉(zhuǎn)速度均為0。

圖1　立管振動變形及坐標示意圖Fig.1 Diagram of the riser vibration amd coordinate

(1)

1.2節(jié)點變形

(2)

(3)

(4)

θI=βI-ΔθθJ=βJ-Δθ

(5)

1.3內(nèi)力計算

然后將I、J節(jié)點內(nèi)力向量分別轉(zhuǎn)換成域坐標分量，節(jié)點彎矩依然不需要轉(zhuǎn)換。

(6)

最后對每一個運動質(zhì)點作內(nèi)力的集成。

1.4質(zhì)點控制方程

考慮結(jié)構(gòu)阻尼，各質(zhì)點運動遵循牛頓第二定律，控制方程為：

(7)

對于立管平面振動問題，其各質(zhì)點有三個位移分量包括x，y軸向的平移分量和繞z軸的轉(zhuǎn)動。上式具體展開表達式為：

(8)

1.5水動力載荷

實際作用于立管上的波浪、海流載荷是一種分布載荷，可以采用傳統(tǒng)有限元的處理方法將其等效到立管各質(zhì)點上。對于單位長度立管所受到的作用力可以采用Morison方程進行計算。

(9)

2數(shù)值求解

采用中央差分的顯式時間積分法進行求解，可以避免隱式解法帶來的迭代和收斂問題。但是時間步長必須進行控制，否則也會導致計算錯誤。對于每一個質(zhì)點運動公式的計算采用第n-1步的運動公式，代入差分公式計算第n步的點位置值xn。然后通過計算出的第n步點位置值xn以及第n步的內(nèi)力公式來計算第n步的節(jié)點內(nèi)力和等效力，如此循環(huán)往復直到模擬計算時刻最后一步。

假設(shè)時間步長為h，那么質(zhì)點的速度和加速度分別可以表示為

(10)

將其代入質(zhì)點運動公式，得到其對應的差分計算公式為：

(11)

式中C1=1/(1+ζh/2)，C2=C1(1-ζh/2)。

3程序編制及算例分析

基于上述理論和計算方法，編制相應的Matlab求解程序，對深水立管動態(tài)響應進行數(shù)值模擬，并將其結(jié)果與以小變形理論為基礎(chǔ)的傳統(tǒng)有限元模擬結(jié)果進行對比分析。傳統(tǒng)有限元模擬中選取梁單元模擬立管，對立管運動方程采用Galerkin法進行有限元離散和Newmark-β法進行數(shù)值求解。本文選取兩種具有不同長度和頂張力的立管進行模擬。立管一：水深300 m，長度300 m，頂端張力500 kN，立管二：水深1 000 m，長度1 000 m，頂端張力1 500 kN，其他主要參數(shù)相同：外徑0.26 m，內(nèi)徑0.2 m，海水密度為1 025 kg/m3, 內(nèi)部流體密度為800 kg/m3，立管管材密度7 850 kg/m3，彈性模量2.07×11N/m4，曳力系數(shù)0.7，慣性系數(shù)2。分析中兩種立管上下兩端均假設(shè)為鉸接，劃分單元數(shù)N都為60。

圖2和圖3繪制了兩種立管在海流流速分別為0.1 m/s，0.2 m/s，0.3 m/s和0.4 m/s四種工況下的最大位移圖。從兩圖中可以看出當流速較小時基于向量式有限元方法(VFIFE)計算的結(jié)果與基于傳統(tǒng)有限元法(FEM)計算的結(jié)果相當吻合，而當流速增大時，兩者之間的差異逐漸變大，并且前者始終小于后者。其主要原因在于向量式有限元方法在變形計算中考慮了立管張力的變化，隨著海流流速的增加，立管位移增大，管內(nèi)軸力也隨之增加，四種工況下立管一的頂端張力分別為500.3 kN、504 kN、518 kN和544 kN，立管二頂端張力分別為1 501 kN、1 510 kN、1 541 kN和1 598 kN，軸力的增大使得立管位移有所減小。

圖2　不同流速下立管一最大位移圖Fig.2 The maximum displacement envelope of the first model riser under different currentvelocity

圖3　不同流速下立管二最大位移圖Fig.3 The maximum displacement envelope of the second model riser under different current velocity

圖4和圖5分別繪制了立管一和立管二在波高3 m周期10 s線性波浪激勵下的最大位移包絡(luò)圖。從圖4中可以看出立管一基于向量式有限元方法計算的結(jié)果與基于傳統(tǒng)有限元法計算的結(jié)果偏小，而從圖5中可以看出立管二基于向量式有限元方法計算的結(jié)果與基于傳統(tǒng)有限元法計算的結(jié)果差異可以忽略。

圖4　波浪作用下立管一最大位移包絡(luò)圖Fig.4 The maximum displacement envelope of the first model riser under wave excited

圖5　波浪作用下立管二最大位移包絡(luò)圖Fig.5 The maximum displacement envelope of the second model riser under wave excited

圖6　立管一100 m和200 m處位移時程及其頻譜Fig.6 Displacement history and frequency spectrum of the first model riser at 100 m and 300 m

圖7　立管一底端和頂端張力時程及其頻譜圖 Fig.7 Bottom and top end tension history and frequency spectrum of the first model riser

圖6～9分別繪制了立管一和立管二兩個節(jié)點處位移時程曲線及其頻譜圖和立管底端和頂端張力時程曲線及其頻譜圖。從圖6和圖8中可以看出立管兩個節(jié)點處的水平振動頻率和橫向激勵頻率保持一致，從圖7和圖9中可以看出軸力是隨著時間周期變化的，頻率為0.2 Hz是其橫向激勵頻率的兩倍，這與有限元分析中軸力保持不變是不同的，顯然軸力的變化與立管實際受力狀態(tài)更加符合。立管一初始狀態(tài)下立管底端軸力是85.4 kN 頂端是500 kN，振動過程中立管底端和頂端軸力最大分別都增加12 kN左右，分別相對增大14%和2.4%，正是軸力的較大變化導致基于向量式有限元計算得到的振動包絡(luò)圖比基于傳統(tǒng)有限元的圖4結(jié)果偏小。立管二初始狀態(tài)下底端軸力是117.5 kN頂端是1 500 kN，振動過程中立管底端和頂端軸力最大分別都增加5.5 kN左右，分別相對增大4.7%和0.37%，此時立管軸力變化較小從而基于向量式有限元計算得到的振動包絡(luò)圖與基于傳統(tǒng)有限元的圖5結(jié)果基本一致。

圖8　立管二300 m和700 m處位移時程及其頻譜圖Fig.8 Displacement history and frequency spectrum of the second model riser at 300 m and 700 m

圖9　立管二底端和頂端張力時程及其頻譜圖Fig.9 Bottom and top end tension history and frequency spectrum of the second model riser

4結(jié)論

對于深水頂張力立管的動力特性及動力響應分析方面，與以小變形理論為基礎(chǔ)的傳統(tǒng)有限元分析相比，本文基于向量式有限元方法分析有以下結(jié)論：

(1) 向量式有限元方法分析中考慮了立管在變形過程中軸向運動及其引起的軸向張力的變化，其振動頻率為橫向激勵頻率的兩倍，這與立管實際受力狀況更加符合。而有限元分析中一般不考慮立管變形過程中軸向張力的變化，也就是始終等于立管靜止狀態(tài)下初始有效張力。

(2) 當海流流速較小時，立管整體變形較小，軸力改變很小，向量式有限元法與有限元結(jié)果基本一致。隨著流速增大，立管整體變形增大，軸力也隨之增大，導致向量式有限元法結(jié)果比有限元結(jié)果偏小。

(3) 對于波浪激勵下的立管，當立管相對較短時，軸力波動幅值相對有效張力較大，向量式有限元方法比有限元結(jié)果偏小。而當立管相對較長時，軸力波動幅值相對有效張力較小，向量式有限元法與有限元結(jié)果基本一致。

通過本文對深水頂張力立管的向量式有限元理論分析和數(shù)值模擬，可以看出向量式有限元法能夠真實準確地模擬立管的振動行為，在立管動態(tài)行為分析中是正確和有效性，并且具有其獨特的優(yōu)勢，點運動公式的描述方式使得控制方程大大簡化，無需求解復雜的非線性聯(lián)立方程組，采用中央差分的顯式時間積分求解避免了迭代和收斂問題，能夠同時計算立管的內(nèi)力(應力應變)和變形。

總之向量式有限元法在處理其它各種形式的懸鏈線立管和柔性立管中將同樣具有顯著優(yōu)勢，可以采用統(tǒng)一的架構(gòu)處理不同的立管結(jié)構(gòu)形式及其復雜的力學特征(小變形大變位、大變形大變位、接觸碰撞等)，為進一步對各種形式的海洋立管復雜行為分析提供了一種新的可行的求解方案。

[ 1 ] Morkookaza C K, Coelho F M, Shiguemoto D A. Dynamic behavior of a top tensioned riser in frequency and time domain [C]//Proceedings of the Sixteenth International Offshore and Polar Engineering Conference, ISOPE, San Francisco, 2006, 2:31-36.

[ 2 ] Chatjigeorgiou I K, Mavrakos S A. Bounded and unbounded coupled transverse response of parametrically excited vertical marine risers and tensioned cable legs for marine applications [J].Applied Ocean Research, 2002, 24(6):341-354.

[ 3 ] Kuiper G L, Brugmans J, Metrikine A V. Destabilization of deep-water risers by a heaving platform[J]. Journal of Sound and Vibration, 2008, 310(3):541-557.

[ 4 ] Chaplin J R, Bearman P W, Huera Huarte F J,et al. Laboratory measurements of vortex-induced vibrations of a vertical tension riser in a stepped current[J]. Journal of Fluids and Structures, 2005, 21(1):3-24.

[ 5 ] 高云,付世曉,宋磊建.柔性立管渦激振動抑制裝置試驗研究[J].振動與沖,2014,33(14):77-83.GAO Yun,FU Shi-xiao,SONG Lei-jian.Experimental investigation on the suppression device of VIV of a flexible riser[J].Journal of Vibration and Shock,2014,33(14):77-83.[ 6 ] 王騰，張修占，朱為全. 平臺運動下深水鉆井隔水管非線性動力響應研究[J]. 海洋工程，2008, 26(3):21-26.WANG Teng, ZHANG Xiu-zhan, ZHU Wei-quan. Study on nonlinear dynamic response of deepwater drilling riser with vessel motion[J].The Ocean Engineering,2008, 26(3):21-26.

[ 7 ] 暢元江，陳國明，許亮斌. 導向架隔水管在波流聯(lián)合作用下的非線性動力響應[J].中國石油大學學報,2006, 30(5):74-77.

CHANG Yuan-jiang, CHEN Guo-ming, XU Liang-bin. Nonlinear dynamic response induced by wave-current for marine risers with guide-frames[J].Journal of China University of Petroleum, 2006, 30(5):74-77.

[ 8 ] Yang He-zhen, Li Hua-jun. Instability Assessment of deep-sea risers under parametric excitation[J].China Ocean Engineering, 2009, 23(4):603-612.

[ 9 ] 唐友剛，邵衛(wèi)東，張杰，等.深海頂張力立管參激-渦激耦合振動響應分析[J].工程力學，2013,30(5):282-286.

TANG You-gang, SHAO Wei-dong, ZHANG Jie,et al. Dynamic response analysis for coupled parametric vibration and vortex-induced vibration of top-tensioned riser in deep-sea[J]. Engineering Mechanics,2013,30(5):282-286.

[10] Chen Wei-min, Li Min, Guo Shuang-xi,et al. Dynamic analysis of coupling between floating top-end heave and riser’s vortex induced vibration by using finite element simulations[J]. Applied Ocean Research, 2014,48:1-9.

[11] 唐國強，滕斌，呂林，等.深海柔性立管渦激振動經(jīng)驗模型建立及應用[J].中國海洋平臺, 2010，25(3):12-16.

TANG Guo-qiang, TENG Bin, Lü Lin,et al. Development and application of an empirical model for vortex-induced vibration of flxible deep water riser[J].China Offshore Platform, 2010, 25(3):12-16.

[12] 丁承先，段元鋒，吳東岳.向量式結(jié)構(gòu)力學[M].北京:科學出版社，2012.

[13] 倪秋斌，段元鋒，高博青.采用向量式有限元的斜拉索振動控制仿真[J].振動工程學報，2014, 27(2):238-245.

NI Qiu-bin，DUAN Yuan-feng，GAO Bo-qing.Vector form intrinsic finite element (VFIFE) based simulation on vibration control of stay cables[J]. Journal of Vibration Engineering,2014, 27(2):238-245.

[14] 喻瑩，羅堯治.基于有限質(zhì)點法的結(jié)構(gòu)倒塌破壞研究Ⅰ：基本方法[J].建筑結(jié)構(gòu)學報，2011,32(11):17-26.

YU Ying, LUO Yao-zhi.Collapse analysis based on finite particle method I: basic approach[J].Journal of Building Structures, 2011,32(11):17-26.

[15] 喻瑩，羅堯治.基于有限質(zhì)點法的結(jié)構(gòu)倒塌破壞研究Ⅱ：關(guān)鍵問題與數(shù)值算例[J].建筑結(jié)構(gòu)學報，2011,32(11):27-35.

YU Ying, LUO Yao-zhi.Collapse analysis based on finite particle methodmethodⅡ:key problems and numerical examples[J]. Journal of Building Structures,2011,32(11):27-35.

[16] 羅堯治，鄭延豐，楊超，等.結(jié)構(gòu)復雜行為分析的有限質(zhì)點法研究綜述[J].工程力學，2014，31(8):1-7.

LUO Yao-zhi, ZHENG Yan-feng, YANG Chao, et al. Review of the finite partical method for complex behavior of structures[J]. Engineering Mechanics,2014,31(8):1-7.

[17] 班自愿，張若京. 懸臂梁大變形的向量式有限元分析[J].計算機輔助工程，2010,19(4): 25-28.

BAN Zi-yuan,ZHANG Ruo-jing.Vector form intrinsic finite element analysis on large deformation of cantilever beam[J]. Coputer Aided Engineering, 2010,19(4): 25-28.

[18] 袁峰.深海管道鋪設(shè)及在位穩(wěn)定性分析[D].杭州：浙江大學，2013.

Dynamic response of a deep-sea top tensioned riser based on vector form intrinsic finite element

LI Xiao-min, ZHANG Lin, NIU Jian-jie, HAN Yong-qing, GUO Hai-yan

(College of Engineering, Ocean University of China, Qingdao 266100, China)

The vector form intrinsic finite element (VFIFE) method also called the finite particle method (FPM) is an innovative method in the structural behavior analysis based on the vector mechanics and numerical calculations. With this method, a structure was discretized into a series of mass particles connected by massless elements only bearing internal forces. The motion of mass paticles obeyed Newton’s second law. Their motion equations were solved using the explicit time domain integration and the central difference method. VFIFE was adoped for the first time to analyze the dynamic behavior of a top tensioned riser (TTR) and the plane beam element was adoped to simulate interaction forces between mass particles. The particle motion equation of the riser was derived and the corresponding solving program was compiled with Matlab. The dynamic behavior of the riser under sea current and wave excitations was calculated and the results were compared with those with the finite element method (FEM). The results showed that the VFIFE method is correct and effective to analyze the dynamic behavior analysis of a riser and it has more advantages than traditional methods do in riser complex behavior analyses.

vector form intrinsic finite element (VFIFE) method; finite element method (FEM); vector mechanics; top tensioned riser (TTR)

10.13465/j.cnki.jvs.2016.11.035

國家自然科學基金(51279187);山東省優(yōu)秀中青年科學家科研獎勵基金(BS2013HZ014)

2015-04-23修改稿收到日期：2015-06-23

李效民 男，博士，講師，1982年3月生

TU311.3

A