隨機荷載下陡波形立管的非線性動力分析

顧洪祿， 李效民， 郭海燕， 崔 鵬， 劉 震， 李福恒

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

近年來，隨著海洋油氣開采逐漸向深水、超深水過渡，簡單鋼懸鏈線立管(steel catenary riser，SCR)頂部張力水平顯著升高，此外，劇烈的頂部浮體運動導致懸掛點、觸地點處極易產生疲勞破壞[1]。因此，SCR可能已經不再適用于深海油氣開采[2]。為改善SCR的疲勞和強度問題，波形立管(緩波形立管、陡波形立管)應運而生，通過在SCR中間管段安裝一定數量的浮力塊，形成拱彎段，波形立管大幅度降低了頂張力的水平，解耦了頂部浮體與觸地點之間的運動。其中，陡波形立管構型陡、跨距小，近地段與海床接近垂直，能有效避免海床穩定性問題。并且，陡波形立管對管內流體變密度的情況適應性較高，在深海油氣開采中具有顯著的優越性[3]。

然而，陡波形立管分析設計需要面臨整體構型復雜、邊界條件處理難度高、計算模型難收斂、內部流體非線性耦合等眾多突出難點。此外，陡波形立管在環境荷載下極易產生大位移、大轉角，具有顯著的幾何非線性，進一步增加了模型建立以及結構動力分析的難度。因此，對陡波形立管進行非線性動力分析具有重要的學術價值和工程實際意義。

國內外學者對頂張力立管、SCR的靜動力響應、管土相互作用等進行了廣泛的研究[4-12]。但是對于波形立管的研究相對較少，且主要集中于緩波形立管。

Li等[13]忽略抗彎剛度的影響，對緩波形立管動力響應進行了初步分析，該研究將各管段簡化為懸鏈線，不能夠考慮海流以及內流的影響。Santillan等[14-15]將陡波形立管簡化為細長彈性桿，利用有限差分法對立管進行了靜動力分析。基于相同的模型，通過在緩S形和陡S形立管拱彎點處施加集中力模擬浮子段，對兩種立管進行了參數敏感性分析，這與陡波形立管強幾何非線性以及浮子段均勻受力的實際情況存在差異。Sun等[16]利用OrcaFlex對陡波形立管的動力響應進行了參數敏感性分析。Yang等[17]對不同參數下的緩波形立管疲勞壽命進行了一系列參數敏感性分析。Wang等[18-19]考慮內部流體流動和管土相互作用，基于小變形梁、非線性大變形梁理論建立了緩波形立管二維非線性動力模型，利用有限差分法求解立管動力響應，然后對緩波形立管靜力性能進行了一系列參數敏感性分析。此外，對緩波形立管在簡諧頂部運動下的動力響應進行了研究[20]，但是以上研究僅僅局限在二維平面內。李艷等[21]應用三維集中質量模型，對緩波形立管在頂部浮體作簡諧運動時的動力響應進行了研究。宋磊建等[22]計算了規則波作用下深海臍帶纜在懸鏈線布局和緩波形布局下張力、曲率的動態響應。

上述研究均基于傳統桿、梁理論建立波形立管力學模型。眾所周知，在傳統桿、梁理論中，需要將物理參數由局部坐標轉換到整體坐標，這會顯著增加計算時間和對計算機CPU的要求。近年來，一種基于絕對坐標的特殊彈性桿理論被廣泛應用于波形立管的數值模擬中。

基于大變形桿理論和有限單元法，Ruan等[23]對彈性海床上的緩波形臍帶纜的靜力力學性能進行了一系列參數分析。基于相同的理論，Ruan等[24]和Cheng等[25]對緩波形立管在頂部浮體運動、波流荷載、內部流體和海床摩擦共同作用下的靜動力響應進行了一系列參數敏感性分析。Kim等[26]對簡單懸鏈線立管與緩波形立管在相同環境荷載下的整體靜動力響應進行了比較分析，但并未對緩波形立管進行深入研究。Qiao等[27]基于大變形梁理論和力平衡原則建立了陡波形立管運動方程，結合有限差分法和打靶法對陡波形立管在海流作用下的靜力性能進行了分析。Liu等[28]基于柔性桿理論建立了陡波形輸流立管的運動方程，并對陡波形立管的靜動力響應進行了一系列參數敏感性分析。

上述研究中，直接在全局坐標系中得到包含所有幾何非線性的立管控制方程，避免了繁瑣的坐標轉換，極大節省了計算資源。本文正是基于此類模型——柔性桿模型[29]建立陡波形立管運動方程。

此外，從上述關于波形立管研究的討論可以看出，目前大多數研究忽略了內部流體的作用，將頂部浮體運動、波浪簡化為簡諧運動，這與實際海況不符。并且目前對于陡波形立管的動力分析十分缺乏，關于陡波形立管在隨機波浪、頂部浮體隨機運動和內部流體共同激勵下的非線性動力分析尚未見報道。

本文基于柔性桿模型，充分考慮陡波形立管幾何非線性，建立了深水陡波形立管運動方程，采用有限元法離散立管運動方程，并利用Newton-Raphson法和Newmark-β法分別計算立管的靜、動力響應，進一步計算得到了陡波形立管在內部流體、隨機波浪、頂部浮體靜態偏移及動態響應等多種因素共同作用下非線性動力響應，并對立管應力水平進行一系列參數敏感性分析。

1 基本理論

1.1 柔性桿模型

如圖1所示，桿的位形由桿軸線位置表示，在三維笛卡爾坐標系中，空間曲線r(s,t)是弧長與時間的函數，表示柔性桿變形后軸線位置狀態。

圖1 柔性桿坐標系示意圖

假設變形前后桿的弧長不發生改變；桿上任意一點的內力狀態可以完全由作用在軸線上的合力和力偶表示；忽略轉動慣量、剪切變形、均布扭矩和分布外力矩的作用；假定桿可伸長且伸長量為小量。根據動量守恒與動量矩守恒可得[29-30]

(1)

M′+r′×F+m=0

(2)

式中:ρ、mf為桿、內部流體單位質量；q為桿單位長度所受外力；m為桿單位長度所受外力矩。

由Euler-Bernoulli梁理論，內力矩M可以表示為

M=r′×EIr″+Hr′

(3)

式中:EI為截面彎曲剛度；H為扭矩，在實際海況下，柔性立管不受均布扭矩作用，其自身扭矩也可忽略。因此其運動方程可以寫為

(4)

當考慮桿的伸長量為小量，可以認為桿件變形條件如下

(5)

式中:EA為桿件軸向剛度；λ為拉格朗日乘子。

綜上，式(4)與式(5)構成柔性桿基本控制方程。

1.2 隨機波浪與頂部浮體運動模擬

由于JONSWAP譜適用范圍更廣，適用于不同成長階段的波浪，本文選取DNV RP-C205 (2007)[31]中推薦的JONSWAP譜

(6)

(7)

式中:Aγ=1-0.287ln(γ)為正則化因子；HS為特征波高；ωp=2π/TP為峰值角頻率,TP為峰值周期；γ為無量綱峰形參數，當γ=1時，JONSWAP譜簡化為P-M譜，本文按照如下公式進行取值

然后基于Longuet-Higgins模型[32]對隨機波浪進行模擬。

頂部浮體運動與海面運動狀況直接相關，頂部浮體對于特征波高的頻率響應一般用運動幅值響應算子進行估算。頂部浮體響應采用Sexton等[33]推薦的運動模型表示

(8)

式中:等式右端第一項表示頂部浮體平均靜偏移;第二項表示一階低頻響應;第三項表示二階高頻響應；SL、TL分別為頂部浮體慢漂單邊幅值、周期；慢漂運動與波浪之間相位差αL一般取0；Sn為頂部浮體對特定波浪(周期Tn=2π/ωn，幅值An)的幅值響應；kn、φn分別為波數、波浪相位角；αn為頂部浮體與波浪周期之間的相位差。

隨機波浪荷載通過莫里森方程求解；由于頂部浮體的尺寸與水深相比可以忽略不計，因此將頂部浮體運動轉化為立管動邊界問題。

1.3 有限元模型

將柔性桿運動方程和變形條件分別寫成張量形式，利用三次Hermit插值函數進行離散，運用Galerkin方法得到矩陣形式的微分方程如下

(9)

(10)

通過Newton-Raphson法和Newmark-β法求解式(9)、式(10)便可以得到立管的靜動力響應。

立管截面應力由下式進行計算

(11)

式中:A為立管橫截面積；E為立管彈性模量。

2 模型驗證

基于第1章所述理論，利用MATLAB編寫陡波形立管非線性動力響應分析程序DRSWR。為驗證本文程序的有效性，選取文獻[34]陡波形立管模型，具體參數如圖2、表1所示。在本文中下降段、浮子段、懸掛段的劃分單元個數分別為10、30、40，進一步細化單元數并不能顯著提高計算精度，并且會極大地消耗計算時間。數值模擬持續時間為500 s，時間步為0.25 s，以上參數在后續計算中不發生變化。利用本文程序對陡波形立管進行分析計算，將本文結果與相同參數下OrcaFlex計算結果進行對比驗證。

表1 陡波形立管參數

圖2 陡波形立管構型(m)

2.1 隨機波浪動力分析驗證

本文選取文獻[35]中波浪數據如下：特征波高Hs=6.5 m，峰值周期TP=12.82 s。利用DRSWR對陡波形立管在隨機波浪作用下的應力包絡圖進行計算，將本文計算結果與同條件下OrcaFlex計算結果進行對比。如圖3所示，可以看出，兩者計算結果吻合程度較好。

(a) 最大應力包絡圖

2.2 頂部浮體動力分析驗證

利用DRSWR對立管在頂部浮體運動下的頂張力時程曲線進行計算。頂部浮體振幅為10 m，周期為27 s，沿X方向運動。本文計算結果與OrcaFlex計算結果進行對比，如圖4所示，兩者計算結果基本一致。

圖4 頂部張力時程變化圖

經過上述對比分析，一定程度上驗證了DRSWR的有效性。

3 陡波形立管非線性動力分析

本文對陡波形立管在隨機波浪、頂部浮體運動共同激勵下的動力響應進行了分析計算。

選用Li等研究中的環境荷載數據，具體參數如下：隨機波浪特征波高HS=6.5 m，峰值周期TP=12.82 s；頂部浮體平均靜偏移S0=5 m，慢漂單邊幅值SL=5 m，慢漂周期TL=200 s，內部流體密度ρf=0，內部流體速度v=0。以上數據為本文計算的基本環境荷載數據，后續計算分析中，無特殊說明，均采用以上環境參數。圖5繪制了DRSWR模擬的隨機波浪高度、頂部浮體運動時程曲線，可以看出兩者運動均具有明顯的隨機性。

(a)

圖6給出了陡波形立管在隨機波浪、頂部浮體作用下的應力包絡圖，可以看出懸垂段、浮子段應力變化幅度較大，而下降段變化幅度較小。這主要是由于波浪入射速度隨水深呈指數衰減，海面處波浪最為劇烈，此外，頂部浮體運動直接作用于懸掛點，因此懸垂段應力變化較為劇烈；由于添加了浮力塊，導致浮子段受到更高水平的水動力荷載，因此該段應力變化幅值較大；浮子段解耦了頂部浮體運動與立管底端之間的動態響應，導致下降段應力幅值變化并不明顯。

由以上分析可以看出陡波構型降低懸掛點處張力水平，減輕頂部浮體負荷的同時，還可以減小井口附近立管應力變化幅值，提高井口處立管疲勞壽命。除此之外，沿立管長度存在兩個應力極大值點，分別位于拱彎點(760 m)、垂彎點(1 192.5 m)處，這主要是兩點初始曲率較大所致。

應力水平是立管強度校核的重要依據，并對立管疲勞損傷具有重要的參考價值。為進一步增進對陡波形立管在隨機波浪、頂部浮體共同作用下動力響應的理解，本文對隨機波浪(特征波高、峰值周期)、頂部浮體運動(平均靜偏移、慢漂運動)、內部流體(內流密度、內流速度)等因素對立管應力的影響進行了參數敏感性分析。

3.1 隨機波浪參數敏感性分析

實際工程中，隨機波浪參數由長期觀測結果整理成的波浪散布圖確定。本文為便于研究波浪參數對動力響應的影響，對參數的選取進行了簡化。

3.1.1 特征波高參數敏感性分析

保持其他參數不變，選取特征波高HS=6.5 m、16.5 m、26.5 m，分別計算陡波形立管最大應力包絡圖。

如圖7所示，分別給出了沿立管長度的最大應力包絡圖以及浮子段、懸垂段最大應力包絡圖的放大圖，由于下降段應力水平較低并且變化幅值較小，不再針對該管段進行相關分析。由圖7可以看出立管整體應力水平隨特征波高的增加而升高，波高越大該趨勢越為明顯；各管段中，懸垂段對該參數變化敏感度最高，并且越靠近懸掛點敏感度越高；此外，隨特征波高的增加，懸垂段應力變化的非線性特性逐漸增強。

(a)

3.1.2 峰值周期參數敏感性分析

保持其他參數不變，分別計算峰值周期TP=2.82 s、7.82 s、12.82 s時立管的最大應力包絡圖，得到峰值周期對立管應力的影響，如圖8所示。

(a)

可以看出立管應力隨峰值周期增加逐漸降低，在低周期下，該趨勢更為顯著；在低周期波浪作用下，立管應力水平急劇升高，且最大應力包絡圖曲線不再圓滑并出現震蕩，應力的非線性明顯增強；懸垂段對峰值周期變化的敏感度仍高于其他管段，且懸掛點附近敏感度最高。

由以上分析可以看出，陡波形立管對隨機波浪的變化敏感度較高，大波高、低周期的隨機波會導致立管整體應力水平急劇升高且表現出顯著的非線性特性；懸垂段對隨機波浪的變化敏感度最高，并且越靠近懸掛點，敏感度越高，主要是由于越接近海面，波浪作用越顯著。因此在對立管進行強度校核時，應充分考慮立管所處海域波浪狀況，并采取必要的消波措施。

3.2 頂部浮體運動參數敏感性分析

3.2.1 平均靜偏移參數敏感性分析

保持其他參數不變，當頂部浮體靜偏移S0=5 m、55 m、105 m時，分別繪制立管最大應力包絡圖，如圖9所示。

(a)

由圖9可以看出，關鍵點(拱彎點、垂彎點、懸掛點)附近應力對該參數變化較為敏感；隨靜偏移的增加，拱彎點、垂彎點附近應力逐漸降低，主要是由于這兩處具有較高的初始曲率，應力水平主要由曲率控制，靜偏移增大會導致這兩點處構型更為平緩，曲率降低；隨靜偏移的增加，懸掛點附近應力升高，是由于懸掛點應力水平主要由張力控制，靜偏移增大，立管整體張緊，導致懸掛點產生更高水平的張力。

3.2.2 慢漂運動參數敏感性分析

圖10、圖11分別繪制了不同慢漂幅值、不同慢漂周期下立管最大應力包絡圖。保持其他參數不變，分別選取慢漂幅值SL=5 m、20 m、35 m，慢漂周期TL=100 s、150 s、200 s。

(a)

(a)

由圖10、圖11可以看出頂部浮體慢漂運動只對懸掛點附近應力水平影響較大，應力水平隨慢漂幅值SL的增加而升高，隨周期TL的增加逐漸降低，主要原因是：頂部浮體直接作用于懸掛點，劇烈的頂部浮體慢漂運動會導致懸掛點附近應力水平升高。

由以上分析可以看出劇烈的隨機波浪、頂部浮體運動都會造成懸掛點處高水平的應力，因此實際工程中應在懸掛點處設置彎曲加強構件或喇叭口等輔助裝置，增強該處的抗彎剛度；此外，頂部浮體慢漂運動對立管整體應力的影響程度遠低于平均靜偏移，這是因為陡波構型提高了對頂部浮體瞬時運動的適應性。

3.3 內部流體參數敏感性分析

3.3.1 內流密度參數敏感性分析

立管在工作期間內部通常會充滿多種流體介質，如石油、天然氣、海水、泥沙等，不同密度的內部流體會對立管產生不同程度影響。本文選取3組不同內流密度：1.29 kg/m3、800.00 kg/m3、998.00 kg/m3，得到了不同內流密度下立管應力響應，如圖12所示。

(a)

由圖12可以看出，內流密度的變化對立管整體應力水平影響顯著；隨內流密度的增加，立管各段均出現明顯的極值點；拱彎點附近區域應力水平隨內流密度的增加逐漸減小；觸地點、垂彎點附近區域變化趨勢與拱彎點變化趨勢正好相反。主要是因為內流密度增加使得立管濕重增加，導致浮子段拱彎部分構型更為平緩，拱彎點曲率降低；懸垂段整體構型彎曲程度變高，垂彎點曲率增大；由于下降段整體濕重的增加，出現一定幅度屈曲，曲率升高。

3.3.2 內流速度參數敏感性分析

立管在運行期間，內部流體以一定速度流動時與立管會產生耦合振動現象，對立管動力響應產生一定程度的影響。本文選取內流密度為800 kg/m3，內流速度為0、5 m/s、10 m/s時，分別計算立管最大應力包絡圖，如圖13所示。

(a)

可以看出，立管整體應力水平隨內流速度的增加而升高，尤其是當內流在高速流動時，趨勢更為明顯。懸垂段對于內流速度變化的敏感度明顯高于其他管段，并且內流速度對立管彎曲段應力存在較為明顯的影響，主要是由于內部流體流動會對彎曲管段內壁產生一定的流動壓力，導致管段產生小幅度的彎曲振動耦合。

由以上分析可以看出，內部流體對各關鍵點處的應力影響顯著，因此在對陡波形立管進行強度校核時，應對內流作用給予足夠的重視。

3.4 多因素變化下陡波形立管動力響應分析

為進一步探究隨機波浪和頂部浮體運動對陡波形立管動力響應的影響，本文選取3組不同的隨機波浪工況(如表2所示)和頂部浮體慢漂幅值(SL=5 m、20 m、35 m)，對多因素共同變化下的關鍵點的最大應力進行了計算分析，計算結果如表3～表5所示。

表2 隨機波浪工況

表3 SL=5 m時陡波形管關鍵點最大應力

表4 SL=20 m時陡波形管關鍵點最大應力

表5 SL=35 m時陡波形管關鍵點最大應力

在不同隨機波浪條件下，關鍵點最大應力都隨慢漂幅值增加逐漸增加；同樣在不同慢漂幅值下，關鍵點最大應力隨特征波高增加逐漸增加。這與3.1.1節單因素變化的結論是一致的。因此，可以看出多因素變化下，特征波高、慢漂幅值對陡波形立管動力響應的影響并沒有發生改變。但是，關鍵點最大應力隨峰值周期變化趨勢與3.1.2節不同，是因為各因素對動力響應的影響作用相反時，各因素之間存在競爭機制，該算例中峰值周期對動力響應的影響明顯小于特征波高和慢漂幅值。

此外，通過對比可以看出，特征波高、慢漂幅值對各個關鍵點影響的程度均為：懸掛點>拱彎點>垂彎點，這與3.1節～3.3節單因素變化下的結論也是一致的。

通過以上分析，我們可以看出，在多因素變化下，如果各因素變化同時增強動力響應，因素之間是相互獨立的，變化趨勢與單因素分析結果是一致的。反之，如果各因素變化趨勢對動力響應的作用相反，單因素變化分析結果在多因素變化下則需要進一步分析確定。

4 結 論

本文基于柔性桿模型，充分考慮陡波形立管的幾何非線性、內部流體非線性耦合，利用MATLAB編寫計算程序DRSWR，對陡波形立管在隨機波浪、頂部浮體共同作用下的非線性動力響應進行了計算分析，并得到了隨機波浪(特征波高、峰值周期)、頂部浮體運動(平均靜偏移、慢漂運動周期、慢漂運動幅值)、內部流體(內流密度、內流速度)等對陡波形立管應力的影響規律，主要結論如下：

(1) 在隨機波浪、頂部浮體共同作用下，陡波形立管浮子段整體應力水平、應力變化幅值最高；應力極值點位于拱彎點、垂彎點。

(2) 陡波形立管對于隨機波浪變化的敏感度較高，懸掛點處應力對隨機波浪變化最為敏感，大波高、低周期的波浪會導致立管整體應力水平急劇升高。

(3) 陡波形立管對于頂部浮體運動敏感度明顯低于隨機波浪，浮體平均靜偏移對立管的影響大于頂部浮體慢漂運動；隨平均靜偏移增大，立管應力極值減小；大幅值、低周期的慢漂運動會引起懸掛點應力顯著增大。

(4) 陡波形立管對內流密度、內流速度表現出較高的敏感度，尤其是各關鍵點處敏感度最高；高密度、高流速的內部流體會使陡波形立管處于高應力狀態。