上凸型內孤立波場中頂張力立管極值響應分析

張 莉 ， 郭海燕 ， 李 朋 ， 王增波

（1.山東科技大學 土木工程與建筑學院，山東 青島266590；2.山東省土木工程防災減災重點實驗室，山東 青島266590；3.中國海洋大學 工程學院，山東 青島 266100；4.海洋石油工程（青島）有限公司，山東 青島 266520）

0 引 言

內孤立波是存在于海洋密度躍層中的非線性大振幅重力波，在其傳播過程中能引起具有較大剪切作用的水平流，從而對柔性海工結構物構成嚴重威脅。海洋立管是海上油氣開發中不可或缺的關鍵構件[1-2]，通常承擔鉆井、完井、修井、輸送油氣等重要作業任務，常見的一種結構形式是頂張力立管（Top Tensioned Riser，簡稱TTR）。頂張力立管的頂部一般會被施加預張力，其上部連接海面平臺，下部連接海底井口的錐形節點或應力節點。貫通海底和海面的頂張力立管，極易受到內孤立波所引起的強烈剪切作用的影響。

目前，已有許多學者就下凹型內孤立波對海工結構作用問題進行了研究，但對上凸型內孤立波的影響研究很少。蔡樹群等[3-4]首次采用Morison公式，計算了內孤立波對小直徑圓柱形樁柱的作用力和力矩，并指出其作用力和力矩遠大于表面波的。謝皆爍等[5]基于二維歐拉方程建立流體非線性連續分層模型，計算了內孤立波作用下小直徑樁柱上的力和力矩，認為不論內孤立波與海流的方向是否一致，前者作用力都大于后者作用力。黃文昊等[6]針對內孤立波與圓柱型結構的相互作用特性開展了系列實驗并建立了內孤立波荷載模型，為常用海洋平臺在內孤立波荷載作用下計算提供了一定依據。蔣武杰等[7]用振型疊加法研究了頂張力立管在下凹型內孤立波與非均勻海流共同作用下的多模態振動，認為下凹型內孤立波與海流耦合效應不可忽視。郭海燕等[8]分析了頂張力立管在下凹型內孤立波、表面波和頂部浮式裝置運動聯合作用下的動力響應情況，指出內孤立波的影響遠比表面波和浮體運動的影響大。樊洪海等[9]基于Euler-Bernoulli理論，計算了內孤立波作用和頂部浮體運動導致自由懸掛立管的動力響應，同樣只考慮了下凹型內孤立波。

然而，實際層化海水中存在上凸型和下凹型兩種內孤立波。在深水區，由于存在淺躍層，內孤立波大多是下凹型；當傳播到淺水區時，在非線性作用下，下凹型的內孤立波會轉化為上凸型[10]。早前在各海域觀測到的下凹型內孤立波較多，近期部分文獻報道我國南海北部發現了大量上凸型內孤立波[11-13]。觀測資料表明，強非線性上凸型內孤立波最大振幅可達40 m[14]，對海工結構物有顯著影響，所以只考慮下凹型內孤立波的影響并不全面。

因此，本文擬建立上凸型內孤立波作用下頂張力立管運動方程，運用有限元方法求解方程，并編制程序對立管運動過程進行模擬。采用我國南海西北部海域的觀測數據，研究上凸型內孤立波引起立管的位移、彎矩、應力和頂底端轉角等動力響應規律，并討論不同因素對立管響應的影響，為海洋立管的工程設計提供參考。

1 數值模型

1.1 頂張力立管在內孤立波作用下控制方程

建立如圖1所示坐標系。以立管未變形的位形為z軸，向上為正，取立管底部為坐標原點，x軸水平向右為正，內孤立波沿ox軸傳播。僅考慮立管受內孤立波來流向的流體力，因而立管的變形和內孤立波引起的流速均在oxz平面發生。頂張力立管裝有張緊設備的上端與上部浮體鉸接，下端與井口的萬向節鉸接，可以將其視為受張力的簡支梁。由于立管頂部張力作用，建模時可以假定立管線彈性、小變形。根據達朗貝爾原理可以得到頂張力立管順流向控制方程[15-16]

其中：mr為單位長度立管質量；mi為單位長度管內流體質量，ma為外部流體附加質量，為附加質量系數，ρe為外部海水的密度，D 為立管外徑；C 為結構阻尼；C′為水動力阻尼，C′=ρeCdDu，Cd為曳力系數，u為內孤立波引起的水平流速；為立管的加速度，為立管的運動速度；E為立管的彈性模量，I為截面慣性矩；Te為有效張力；f′為流體作用力為內孤立波引起的海水加速度，為海水相對于立管結構的速度。

1.2 內孤立波數值模型

假定海水密度的垂向分布是兩層結構，這種近似雖然簡單，但卻是一般中、低緯度淺海區域春夏秋三季層化狀況的粗略近似[10]，具有實際意義。設上層海水的密度和深度分別為ρ1和h1，下層海水的密度和深度分別為ρ2和h2，水深h=h1+h2，內孤立波振幅為η0。坐標如圖1所示，根據兩層密度分布假設，有KdV理論解內孤立波界面位移表達式[11]（不計淺水項和耗散項）：

其中：η0為內孤立波振幅；φ是內孤立波相位角，φ=為內孤立波相速度，，Δρ、ρ是上下兩層流體的密度差與平均密度，分別為：；l為特征長度，l≈

由水動力學、流體連續方程結合波面方程可以得到內孤立波引起的水平方向速度和加速度。

圖1 頂張力立管分析模型及受力情況Fig.1 Sketch of a TTR under loads

內孤立波引起的水平速度：

內孤立波引起的水平加速度：

1.3 運動方程的邊界條件

立管模型頂底端均為鉸接，因而邊界條件分別為：

1.4 有限元離散及方程求解

采用Galerkin方法對立管控制方程進行離散，最終得到單元剛度方程：

將各單元剛度矩陣組合至整體剛度矩陣，即得到頂張力立管在內波場中的有限元方程組：

其中：[M]、[C]、[K]分別為立管系統的整體質量矩陣、整體阻尼矩陣和整體剛度矩陣；為立管每一個節點上的順流向加速度、速度、位移列陣；｛F｝為整體荷載列陣。

在時域內求解方程，采用Newmark-β法進行逐步積分求解，即只要給出ti時刻的位移xi、速度和加速度，可計算出經過Δt時間后 ti+1時刻的解 xi+1、和，進而可推算出內孤立波作用下頂張力立管所有時刻的位移、彎矩、應力、轉角等動力響應。

2 計算過程與程序驗證

基于上述計算方法，用MATLAB編制相應計算程序，輸入上凸型內孤立波要素和立管的相關參數，即可得到立管的響應。計算流程如圖2所示。

圖2 上凸型內孤立波作用下頂張力立管分析計算流程圖Fig.2 Flow chart of the riser dynamic analysis under elevation ISW

因上凸型內孤立波發生時間不可預見并且觀測難度較大，目前尚無上凸型內孤立波經過時立管響應的實測數據，所以用文獻結果與本文進行對比驗證。

（1）采用文獻[14]中的內孤立波參數，將本文程序計算得到的水平流速與文獻進行對比，如圖3所示。其中（a）為文獻結果，（b）為本文計算結果，由圖可得，上下兩層流體中水平流速隨時間變化的情況同步，程序計算所得最大流速（上層0.57 m/s、下層0.71 m/s）與文獻結果基本相等，流速方向也完全一致。

（2）為驗證流體力作用下頂張力立管動力響應計算的正確性，采用文獻[17]中的數據對線性波浪作用下頂張力立管位移進行分析，并將計算結果與文獻進行對比，如圖4所示。其中，（a）圖中帶“×”標記的線條為文獻中的時域分析結果，本文時域計算結果（b）圖與其振動形態近似，無量綱最大位移包絡圖也能較好地吻合，從而說明本文方法及計算程序可靠。

圖3 內孤立波引起的水平流速與文獻的對比Fig.3 Comparison of the simulated horizontal flow velocity induced by ISW with the result of other paper

圖4 波浪作用下立管響應與文獻的對比Fig.4 Comparison of the riser response under wave with the result of other paper

3 算例分析

3.1 算例數據

假設立管均質且壁厚沿管長不變，沿用上文中各變量的定義，立管參數見表1（表中CM和Cd的選取依據規范API-RP-16Q）。流場參數采用文獻[14]中實際觀測到的上凸型內孤立波數據，見表2。

表1 立管參數Tab.1 Dimensions and coefficients of the TTR

表2 上凸型內孤立波參數Tab.2 Properties of elevation ISW

3.2 上凸型內孤立波引起的立管響應

圖5為計算得到的上凸型內孤立波經過時立管無量綱位移包絡圖和彎矩包絡圖。由圖5可知，立管在上下層流體中段附近各有一個極值，下層流體中的響應大于上層流體中的響應，兩層流體分界面附近位移響應最小。全長最大位移在距海底32.17 m處，其絕對值為1.77D；全長最大彎矩處于距海底27.79 m 處，為 33.45 kN·m。

圖5 頂張力立管無量綱位移包絡圖及彎矩包絡圖Fig.5 Envelopes of dimensionless displacement and moment for the TTR

立管設計時要求最大等效應力不大于屈服應力的2/3[18]，因而本文對上凸型內孤立波作用下頂張力立管的應力進行了計算，得到應力包絡如圖6所示。由圖6分析可知，立管應力分別在上下層有一個極值，最大應力發生在距海底27.79 m處，為36.56 MPa。立管張力從頂端向下逐漸減小，因而位移、彎矩和應力極值均出現在立管偏下部分。

為防止節點破壞，立管設計規范API中規定鉆井模式下的立管頂底部柔性接頭平均轉角不大于2°，最大轉角不大于4°[18]。當上凸型內孤立波經過時立管頂底端轉角隨時間變化的情況如圖7所示。由圖7可知，頂端轉角先增大后減小為0°，方向有一次反轉，這與上層流體中的立管位移方向有一次變化相對應；底端轉角經歷由0°逐漸增大到極值再緩慢減小至0°的過程，頂端極值為0.56°，底端極值為1.54°。如果考慮內孤立波經過時上部浮體運動所引起的立管運動，立管轉角可能會超過設計允許值，應當引起注意。

圖6 頂張力立管應力包絡圖Fig.6 Envelopes of stress for the TTR

圖7 立管頂底端轉角時程曲線Fig.7 Time history of the angles at top and bottom ends

圖8為距海底分別為16 m、32 m、48 m、65 m及91 m處立管位移的時程曲線。其中16 m、32 m和48 m均處于下層流體，其位移隨時間的變化基本同步，均在波峰經過時達到極值，并且都沿-x方向發生位移；65 m處為上下兩層流體的分界線，此處位移較小，在波峰經過前其位移方向基本同下層流體中的立管（-x方向），但在波峰將要離開時發生反方向位移；91 m處為上層流體的中間，其位移（-x方向）首先增大后減小，在波峰到達之前位移減為0，接著發生反方向位移；很明顯，位于兩層流體分界面及其以上部位的立管位移方向發生了一次逆轉。

從圖8可看出，立管的位移主要發生在400-800 s之間，為更全面地研究立管位形變化情況，圖9呈現了這段時間內的立管全長位移。由圖9分析可知，在上凸型內孤立波逐漸靠近立管直至離開的過程中，位于下層流體中的立管位移逐漸增加到最大，再逐漸減小至0，位移方向始終與下層流速方向相同（-x方向）；由于下層流速略大于上層流速，位于上層流體中的立管首先發生與下層流速同向的位移，隨著上層流速的逐漸增大，流體作用力導致立管上半部分-x方向位移逐漸減小并發生反向位移（x方向），整個立管呈現“S”型；隨著內孤立波逐漸遠離，最終全長位移歸為0。

圖8 不同深度處立管位移時程圖Fig.8 Time history of displacement in different depth

圖9 400 s、500 s、600 s、700 s和 800 s時立管全長位移圖Fig.9 Displacements in 400 s,500 s,600 s,700 s and 800 s

3.3 不同因素對頂張力立管極值響應的影響

不同振幅內孤立波引起的流場不同，進而導致立管響應有差異。分別選取內孤立波振幅等于15 m、25 m、35 m、45 m和55 m，得到立管全長的位移狀態，如圖10所示（Amplitude表示內孤立波的振幅）。當振幅等于15 m、25 m、35 m、45 m和55 m時，對應的立管全長最大位移分別為0.30D、0.90D、1.77D、2.93D和4.37D。可見位移響應隨著內孤立波振幅的增大而增大，位于下層流體中的立管位移增幅大于上層流體中的。

通常海水密度并不恒定，上下層流體之間的密度差也會影響立管極值響應。圖11討論了密度差的影響，上層密度1 022 kg/m3保持不變，令下層密度分別為1 024 kg/m3、1 025 kg/m3、1 026 kg/m3、1 027 kg/m3和 1 028 kg/m3（即上下層密度差為 2 kg/m3、3 kg/m3、4 kg/m3、5 kg/m3和 6 kg/m3），得到立管的最大位移分別為0.88D、1.33D、1.77D、2.20D和2.66D。很明顯，隨著上下層流體密度差的增大，立管全長位移也跟著增大，上下層流體中的位移方向依然相反。從公式（3）中也可看出，密度差增大會導致上下層流速同時增大，從而流體對立管作用力也增大。

作為頂張力立管設計過程中的重要參數之一，頂部張力的大小直接影響立管的響應大小。圖12表示頂部施加張力不同時立管全長位移，當頂部張力分別為1.1Te、1.3Te、1.5Te、1.7Te和1.9Te時，立管最大位移分別為1.77D、1.57D、1.41D、1.29D和1.19D，立管下半段對頂張力的變化更敏感。隨著頂部張力的增大，相當于增加了立管的剛度，位移相應減少。因而在內波頻發海域，可通過增加張緊器中的張力來改善立管的性能，但張力的增加會引起頂部應力集中。

圖10 不同內孤立波振幅引起的全長位移Fig.10 Displacements under different amplitudes of ISW

圖11 不同海水密度比下立管的全長位移Fig.11 Displacements under different density ratio

圖12 不同頂部張力下立管的全長位移Fig.12 Displacements with different top tension

圖13 不同壁厚下立管的全長位移Fig.13 Displacements with different wall thickness

圖13為考慮立管壁厚對立管位移響應的影響。壁厚分別選取3.0 cm、2.5 cm、2.0 cm、1.5 cm和1.0 cm，立管最大位移分別為0.86D、1.00D、1.28D、1.77D和3.00D。立管的壁厚是一個重要設計參數，設計可能有內孤立波經過海域的立管時，在允許范圍內，立管壁厚應盡可能厚一些，以增加其剛度，減小位移。

4 結 論

本文基于KdV方程，參照南海實測數據，對上凸型內孤立波作用下頂張力立管的動力響應做了數值模擬，采用Newmark-β法實現有限元模擬的時域分析，得到了上凸型內孤立波作用下頂張力立管的無量綱位移、彎矩、應力包絡及頂底端轉角；并對上凸型內孤立波振幅、上下層流體密度之差、頂張力和壁厚對立管響應的影響進行了分析。

研究結果表明：上凸型內孤立波使位于上下兩層流體中的頂張力立管位移、彎矩、應力及頂底端轉角增大，但在兩層海水分密度界面附近影響相對較小，另外，內孤立波對立管中下部的受力及響應有顯著影響。位于下層流體中的立管僅發生-x方向位移，位于上層流體中的立管先發生短時小幅度-x方向位移再反方向運動。由于這一密度分層海域上下層水深相當，上凸型內孤立波引起的剪切作用劇烈，導致立管位形呈現"S"型，而增大立管頂部的張力或者增大立管壁厚能減輕上凸型內孤立波的作用。