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

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

0 引 言

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

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

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

因此，本文擬建立上凸型內(nèi)孤立波作用下頂張力立管運動方程，運用有限元方法求解方程，并編制程序?qū)α⒐苓\動過程進行模擬。采用我國南海西北部海域的觀測數(shù)據(jù)，研究上凸型內(nèi)孤立波引起立管的位移、彎矩、應(yīng)力和頂?shù)锥宿D(zhuǎn)角等動力響應(yīng)規(guī)律，并討論不同因素對立管響應(yīng)的影響，為海洋立管的工程設(shè)計提供參考。

1 數(shù)值模型

1.1 頂張力立管在內(nèi)孤立波作用下控制方程

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

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

1.2 內(nèi)孤立波數(shù)值模型

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

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

由水動力學(xué)、流體連續(xù)方程結(jié)合波面方程可以得到內(nèi)孤立波引起的水平方向速度和加速度。

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

內(nèi)孤立波引起的水平速度：

內(nèi)孤立波引起的水平加速度：

1.3 運動方程的邊界條件

立管模型頂?shù)锥司鶠殂q接，因而邊界條件分別為：

1.4 有限元離散及方程求解

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

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

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

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

2 計算過程與程序驗證

基于上述計算方法，用MATLAB編制相應(yīng)計算程序，輸入上凸型內(nèi)孤立波要素和立管的相關(guān)參數(shù)，即可得到立管的響應(yīng)。計算流程如圖2所示。

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

因上凸型內(nèi)孤立波發(fā)生時間不可預(yù)見并且觀測難度較大，目前尚無上凸型內(nèi)孤立波經(jīng)過時立管響應(yīng)的實測數(shù)據(jù)，所以用文獻結(jié)果與本文進行對比驗證。

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

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

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

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

3 算例分析

3.1 算例數(shù)據(jù)

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

表1 立管參數(shù)Tab.1 Dimensions and coefficients of the TTR

表2 上凸型內(nèi)孤立波參數(shù)Tab.2 Properties of elevation ISW

3.2 上凸型內(nèi)孤立波引起的立管響應(yīng)

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

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

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

為防止節(jié)點破壞，立管設(shè)計規(guī)范API中規(guī)定鉆井模式下的立管頂?shù)撞咳嵝越宇^平均轉(zhuǎn)角不大于2°，最大轉(zhuǎn)角不大于4°[18]。當上凸型內(nèi)孤立波經(jīng)過時立管頂?shù)锥宿D(zhuǎn)角隨時間變化的情況如圖7所示。由圖7可知，頂端轉(zhuǎn)角先增大后減小為0°，方向有一次反轉(zhuǎn)，這與上層流體中的立管位移方向有一次變化相對應(yīng)；底端轉(zhuǎn)角經(jīng)歷由0°逐漸增大到極值再緩慢減小至0°的過程，頂端極值為0.56°，底端極值為1.54°。如果考慮內(nèi)孤立波經(jīng)過時上部浮體運動所引起的立管運動，立管轉(zhuǎn)角可能會超過設(shè)計允許值，應(yīng)當引起注意。

圖6 頂張力立管應(yīng)力包絡(luò)圖Fig.6 Envelopes of stress for the TTR

圖7 立管頂?shù)锥宿D(zhuǎn)角時程曲線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均處于下層流體，其位移隨時間的變化基本同步，均在波峰經(jīng)過時達到極值，并且都沿-x方向發(fā)生位移；65 m處為上下兩層流體的分界線，此處位移較小，在波峰經(jīng)過前其位移方向基本同下層流體中的立管（-x方向），但在波峰將要離開時發(fā)生反方向位移；91 m處為上層流體的中間，其位移（-x方向）首先增大后減小，在波峰到達之前位移減為0，接著發(fā)生反方向位移；很明顯，位于兩層流體分界面及其以上部位的立管位移方向發(fā)生了一次逆轉(zhuǎn)。

從圖8可看出，立管的位移主要發(fā)生在400-800 s之間，為更全面地研究立管位形變化情況，圖9呈現(xiàn)了這段時間內(nèi)的立管全長位移。由圖9分析可知，在上凸型內(nèi)孤立波逐漸靠近立管直至離開的過程中，位于下層流體中的立管位移逐漸增加到最大，再逐漸減小至0，位移方向始終與下層流速方向相同（-x方向）；由于下層流速略大于上層流速，位于上層流體中的立管首先發(fā)生與下層流速同向的位移，隨著上層流速的逐漸增大，流體作用力導(dǎo)致立管上半部分-x方向位移逐漸減小并發(fā)生反向位移（x方向），整個立管呈現(xiàn)“S”型；隨著內(nèi)孤立波逐漸遠離，最終全長位移歸為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 不同因素對頂張力立管極值響應(yīng)的影響

不同振幅內(nèi)孤立波引起的流場不同，進而導(dǎo)致立管響應(yīng)有差異。分別選取內(nèi)孤立波振幅等于15 m、25 m、35 m、45 m和55 m，得到立管全長的位移狀態(tài)，如圖10所示（Amplitude表示內(nèi)孤立波的振幅）。當振幅等于15 m、25 m、35 m、45 m和55 m時，對應(yīng)的立管全長最大位移分別為0.30D、0.90D、1.77D、2.93D和4.37D?？梢娢灰祈憫?yīng)隨著內(nèi)孤立波振幅的增大而增大，位于下層流體中的立管位移增幅大于上層流體中的。

通常海水密度并不恒定，上下層流體之間的密度差也會影響立管極值響應(yīng)。圖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）中也可看出，密度差增大會導(dǎo)致上下層流速同時增大，從而流體對立管作用力也增大。

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

圖10 不同內(nèi)孤立波振幅引起的全長位移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為考慮立管壁厚對立管位移響應(yīng)的影響。壁厚分別選取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。立管的壁厚是一個重要設(shè)計參數(shù)，設(shè)計可能有內(nèi)孤立波經(jīng)過海域的立管時，在允許范圍內(nèi)，立管壁厚應(yīng)盡可能厚一些，以增加其剛度，減小位移。

4 結(jié) 論

本文基于KdV方程，參照南海實測數(shù)據(jù)，對上凸型內(nèi)孤立波作用下頂張力立管的動力響應(yīng)做了數(shù)值模擬，采用Newmark-β法實現(xiàn)有限元模擬的時域分析，得到了上凸型內(nèi)孤立波作用下頂張力立管的無量綱位移、彎矩、應(yīng)力包絡(luò)及頂?shù)锥宿D(zhuǎn)角；并對上凸型內(nèi)孤立波振幅、上下層流體密度之差、頂張力和壁厚對立管響應(yīng)的影響進行了分析。

研究結(jié)果表明：上凸型內(nèi)孤立波使位于上下兩層流體中的頂張力立管位移、彎矩、應(yīng)力及頂?shù)锥宿D(zhuǎn)角增大，但在兩層海水分密度界面附近影響相對較小，另外，內(nèi)孤立波對立管中下部的受力及響應(yīng)有顯著影響。位于下層流體中的立管僅發(fā)生-x方向位移，位于上層流體中的立管先發(fā)生短時小幅度-x方向位移再反方向運動。由于這一密度分層海域上下層水深相當，上凸型內(nèi)孤立波引起的剪切作用劇烈，導(dǎo)致立管位形呈現(xiàn)"S"型，而增大立管頂部的張力或者增大立管壁厚能減輕上凸型內(nèi)孤立波的作用。