劉德鵬， 艾尚茂， 孫麗萍， 賈魯生

(1.哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001； 2.中海油研究總院有限責(zé)任公司 工程研究設(shè)計院， 北京 100028)

海洋立管的作用是將油氣資源從海底輸送到海上平臺，是海洋油氣行業(yè)中的關(guān)鍵水下設(shè)備。立管內(nèi)部通常為混合了多種氣體與液體的多相流體，該種流動情況被稱為段塞流。由于立管結(jié)構(gòu)的幾何非線性[1-2]以及段塞流不均勻的質(zhì)量分布，使得輸送段塞流的立管動態(tài)響應(yīng)預(yù)報十分復(fù)雜。

考慮到流固耦合方法計算復(fù)雜且耗時較長，一般在立管設(shè)計中需要建立合適的段塞流數(shù)學(xué)模型來開展海洋立管總體響應(yīng)預(yù)報。段塞流具有速度和密度隨時間變化的特征[3]，當內(nèi)部流動速度一致且每部分段塞長度相等時，即為穩(wěn)態(tài)段塞流。早期研究往往將立管內(nèi)部流動簡化為均勻流考慮[4]。Patel 等[5]提出了一種密度隨位置和時間的變化而呈正弦變化的穩(wěn)態(tài)段塞流模型。當前國內(nèi)對于段塞流的研究主要集中在嚴重段塞流的成型機理以及影響[6]，穩(wěn)態(tài)段塞流模型相關(guān)研究較少。

Shabana[7]提出了一種有限元建模方法，即絕對節(jié)點坐標法(absolute nodal coordinate formula, ANCF)。該方法中節(jié)點坐標在全局坐標系中定義，有效地避免了復(fù)雜的坐標變換。此外，ANCF采用斜率向量代替角坐標的方式來描述結(jié)構(gòu)的局部旋轉(zhuǎn)[8]，可以準確地描述剛體位移與柔性變形之間的耦合關(guān)系。通過將ANCF與任意拉格朗日-歐拉(arbitrary Lagrangian-Eulerian, ALE)描述結(jié)合，可以精確高效地描述具有大變形特征的介質(zhì)[9]。

為預(yù)測內(nèi)部段塞流引起的立管結(jié)構(gòu)振動，本研究在ALE-ANCF框架下建立了海洋立管與內(nèi)部段塞流的有限元分析模型。以鋼懸鏈立管為例，分析了內(nèi)部穩(wěn)態(tài)段塞流對立管動態(tài)響應(yīng)的影響。

1 ALE-ANCF有限元理論

1.1 控制體ALE-ANCF描述

ALE采用物質(zhì)坐標來描述控制體的邊界和沿介質(zhì)的物質(zhì)流動，介質(zhì)上任一點的物質(zhì)坐標為該點與起始點之間的介質(zhì)長度。由于控制邊界s1和s2是隨時間變化的，控制區(qū)間的位置及長度隨之變化。為便于后續(xù)控制區(qū)間的描述，定義局部坐標ξ為[10]：

(1)

在任意時刻t時控制體內(nèi)部任意物質(zhì)點的空間坐標為：

r=r(ξ,t)

(2)

r=Nw,r′=N′w

(3)

(4)

1.2 穩(wěn)態(tài)段塞流單元

采用Patel等[5]提出的穩(wěn)態(tài)段塞流模型來模擬內(nèi)部流體密度的時空變化，如圖1所示，其中t、n分別為管道切向與法向，Di為立管內(nèi)徑，vi為段塞速度，ls為段塞長度。該數(shù)學(xué)模型中段塞流密度分為恒定分量與時變分量：

ρf=ρ0(1+ηei(ks-ωt+ψ))

(5)

式中：ρf和ρ0分別為為內(nèi)部流體密度和恒定密度分量；k=2π/ls為段塞流波數(shù)；η為波動幅值系數(shù)；ω為段塞圓頻率；ψ為初始相位，用于指定初始時刻的段塞流質(zhì)量分布。

圖1 段塞流模型Fig.1 Configuration of slug flow

段塞流單元可以在管道結(jié)構(gòu)的切向上任意移動，形函數(shù)僅約束其在管道結(jié)構(gòu)法向的位移。假定內(nèi)部流動為穩(wěn)態(tài)段塞流，段塞速度vi=ω/k，此時對于任意段塞流單元，物質(zhì)的流入速度與流出速度相等。段塞流單元邊界的物質(zhì)坐標為段塞速度的函數(shù)：

s1=-vit,s2=Le-vit

(6)

式中Le為單元長度。將式(6)代入式(4)，即可得到段塞流單元內(nèi)部任意點的速度與加速度。

段塞流單元的平衡方程為：

(7)

式中：Qi為廣義慣性力；Mf為單元質(zhì)量陣。根據(jù)達朗貝爾原理，Qi與Mf可由虛功方程得到：

(8)

(9)

1.3 立管單元

(10)

式中：Mr、Qf與Qe分別為管道單元的質(zhì)量矩陣、廣義外力以及廣義粘彈性力。根據(jù)達朗貝爾原理，Mr、Qf與Qe由虛功方程得：

(11)

(12)

(13)

式中：ρr為立管材料密度；Ar為立管截面積；f為單元外力矢量；ε和κ分別為單元的軸向應(yīng)變和曲率，定義為：

(14)

單元的剛度矩陣可由Qe對w直接求導(dǎo)得到[11-12]:

(15)

式(15)右端兩部分分別表示剛度矩陣的拉伸分量和彎曲分量。

施加在立管單元上的外力包括重力、浮力、拖曳力和慣性力，后兩項水動力通過莫里森公式計算：

(16)

(17)

1.4 數(shù)值求解方法

靜態(tài)平衡方程的求解采用Newton-Raphson方法，迭代形式為[13]：

w(n+1)=w(n)+Δw(n)

Δw(n)=[K(w(n))](Qe-Qf+Qr)

(18)

式中：n為迭代次數(shù)，當Δw滿足：

‖Δw(n)‖≤w0

(19)

時迭代終止。

在t+Δt時刻，輸送段塞流的立管的平衡方程為：

(20)

式(20)可采用Newmark法進行求解，在該方法中，t+Δt時刻的位移、速度與加速度為：

(21)

將式(21)代入t+Δt時刻的平衡方程，可以得到遞推公式：

(22)

其中：

2 穩(wěn)態(tài)段塞流對立管的影響分析

2.1 立管與段塞流參數(shù)

以鋼懸鏈立管(SCR)為例，分析穩(wěn)態(tài)段塞流對立管的影響。圖2給出了鋼懸鏈立管模型的示意圖，該立管總長1 260 m，工作水深為900 m，水平跨距700 m，頂部懸掛點位于水面下10 m。立管外徑D=0.4 m與內(nèi)徑Di=0.36 m，空氣中單位長度的質(zhì)量為189.1 kg/m，法向拖曳力系數(shù)與慣性力系數(shù)分別為1與2。

圖2 鋼懸鏈立管Fig.2 A schematic presentation of the SCR model

鋼懸鏈立管內(nèi)部輸送流體為API 48原油和甲烷氣體，內(nèi)部流體平均密度選為632.135 kg/m3，液相密度與氣相密度分別為790 kg/m3與0.675 kg/m3，波動幅值系數(shù)為0.25，段塞速度為2.5 m/s，段塞長度Ls=80 m。

為了探討穩(wěn)態(tài)段塞流對鋼懸鏈立管的影響，本研究還對輸送內(nèi)部均勻流的立管進行了模擬，該內(nèi)部均勻流模型具有與穩(wěn)態(tài)段塞流相同的平均密度與流動速度，2種內(nèi)部流動情況下的內(nèi)部流體總質(zhì)量相同。

2.2 靜態(tài)分析

在分析穩(wěn)態(tài)段塞流對立管的動態(tài)響應(yīng)影響之前，先簡要分析內(nèi)部均勻流與段塞流情況下立管的靜態(tài)特性。依據(jù)單個段塞的長度，將鋼懸鏈立管劃分為315個等長單元。采用土壤剛度來描述立管與海床之間的相互作用[14]。圖3給出了2種內(nèi)部流動情況下的立管張力與彎矩特性曲線。

通過觀察圖3可以發(fā)現(xiàn)2種內(nèi)部流動情況下的立管張立曲線幾乎完全重合，這說明段塞流的不均勻密度分布對立管的張力影響很小。2種內(nèi)部流動情況下的彎矩曲線具有相似的趨勢，但是內(nèi)部穩(wěn)態(tài)段塞流情況下的彎矩曲線呈現(xiàn)出明顯的局部波動。其原因是鋼懸鏈立管的截面抗彎剛度極大，導(dǎo)致很小的曲率變形就會引起很大的彎矩變化。

圖3 立管靜態(tài)張力與彎矩Fig.3 Static tension and bending moment along the SCR

2.3 動態(tài)分析

在動態(tài)仿真中，通過對立管的頂部節(jié)點施加X方向的強迫位移來模擬頂部平臺的波頻運動。激勵位移rx可以表示為：

(23)

式中：頂部激勵位移幅值λ和周期Τm分別取為0.5 m和9 s。動態(tài)仿真時長設(shè)置為2 000 s, 時間步長為0.05 s。

圖4給出了鋼懸鏈立管頂部張力的時歷曲線。可以看出在內(nèi)部穩(wěn)態(tài)段塞流的情況下，立管頂部張力的變化呈現(xiàn)出雙頻疊加的特性。與內(nèi)部均勻流情況相比，穩(wěn)態(tài)段塞流情況下頂部張力的波動幅值有所增大，2種情況下頂部張力幅值之間的差值為3.78 kN，僅為內(nèi)部均勻流情況下頂部靜態(tài)張力的0.27%，這說明穩(wěn)態(tài)段塞流對鋼懸鏈立管的動態(tài)張力影響很小。

圖4 頂部張力時歷曲線Fig.4 Time histories of the top tension

圖5給出了2種內(nèi)部流動情況下沿管長彎矩分布時歷圖。比較圖5(a)、(b)可以發(fā)現(xiàn)彎矩在穩(wěn)態(tài)段塞流情況下的變化幅值明顯大于內(nèi)部均勻流情況。立管觸地段的彎矩在2種流動情況下的變化趨勢相近，而在立管懸掛段，彎矩的變化頻率與段塞流頻率相近。

圖5 沿管長彎矩分布時歷Fig.5 Distribution of bending moments along the arc length of the riser

為更明顯地觀察彎矩變化規(guī)律，分別選擇觸地點附近的節(jié)點35(距錨點140 m處)與懸掛段彎矩最大處附近的節(jié)點75(距錨點300 m處)作為觀測點，圖6給出了上述2節(jié)點處彎矩變化(距錨點140 m處)的時歷曲線。可以看出，穩(wěn)態(tài)段塞流情況下上述兩節(jié)點的彎矩變化均呈現(xiàn)出雙頻疊加的特性，節(jié)點35處的彎矩變化以頂部激勵頻率為主，而節(jié)點75處的彎矩變化以穩(wěn)態(tài)段塞流的頻率為主，這說明穩(wěn)態(tài)段塞流對懸掛段彎矩的影響更為顯著。

圖6 節(jié)點35與75處彎矩時歷曲線Fig.6 Time histories of bending moment at node 35 and 75

位移計算結(jié)果如圖7所示。結(jié)果顯示，段塞流對節(jié)點35的X向位移幾乎沒有影響，而其Z向位移會隨著段塞流的不均勻質(zhì)量移動發(fā)生波動，2種內(nèi)部流動情況下的最大Z向位移值相差0.01D。穩(wěn)態(tài)段塞流情況下節(jié)點75的X與Z向位移均出現(xiàn)了雙頻疊加的特性，與均勻流動情況下的最大位移值相差分別為0.07D和0.08D。上述分析表明，穩(wěn)態(tài)段塞流對立管觸地段影響極小，但是會改變其Z向的運動特性。而對于立管懸掛段，穩(wěn)態(tài)段塞流的影響顯著，這也是該段彎矩變化特性發(fā)生改變的原因。

圖7 節(jié)點35與75處位移時歷曲線Fig.7 Time history of displacements at node 35 and 75

3 段塞長度參數(shù)分析

為分析段塞長度對立管動態(tài)響應(yīng)的影響，對輸送不同長度段塞流的鋼懸鏈立管進行了時域仿真。段塞長度與管徑之間的比值選取為Ls/D=100，160，220，280。內(nèi)部流體平均密度與流動速度保持不變圖8和圖9分別給出了不同段塞長度下的節(jié)點35與節(jié)點75處的位移時歷曲線。

圖8 不同段塞長度下節(jié)點35處位移時歷曲線Fig.8 Time histories of displacements at node 35 under different slug length

由圖8(a)可知，段塞長度的改變對觸地點X方向的位移沒有明顯影響。對于圖8(b)，可以發(fā)現(xiàn)隨著段塞長度的逐漸增加，觸地點的Z向位移幅值顯著增大。此外，在液塞段經(jīng)過時，觸地點的Z向位移受到抑制，而氣塞段經(jīng)過時則與之相反，這一現(xiàn)象是由液塞與氣塞交替流動導(dǎo)致的立管內(nèi)部流體質(zhì)量不斷變化而引起的。由圖9(a)、(b)可以看出，段塞長度的改變會導(dǎo)致立管懸掛段位移幅值的變化，隨著段塞長度的逐漸增加，懸掛段X與Z方向位移的雙頻疊加特性愈為明顯。

圖9 不同段塞長度下節(jié)點75處位移時歷曲線Fig.9 Time histories of displacements at Node 75 under different slug length

不同段塞長度下的彎矩時歷曲線如圖10所示。由圖10(a)可以看出，隨著段塞長度的增加觸地點彎矩明顯增大，其變化的雙頻疊加特性更加明顯，當Ls/D=100,280時，彎矩平均值分別為101.06 kN與127.53 kN,兩者之間差值為觸地點靜態(tài)彎矩的36%。對于圖10(b)，隨著段塞長度的增加，懸掛段彎矩的波動幅值逐漸增大，其變化的主導(dǎo)頻率與段塞流頻率相接近。

圖10 不同段塞長度下的彎矩時歷曲線Fig.10 Time histories of bending moment under different slug lengths

4 結(jié)論

1) 穩(wěn)態(tài)段塞流的不均勻密度分布對立管的張力影響很小，但是對立管懸掛段的彎矩影響顯著，彎矩變化主導(dǎo)頻率與段塞流頻率相近。

2) 穩(wěn)態(tài)段塞流會導(dǎo)致鋼懸鏈立管位移幅值的增大，液塞與氣塞交替流動導(dǎo)致的內(nèi)流質(zhì)量變化會改變觸地點垂向的運動特性。

3) 隨著段塞長度的增加，立管懸掛段位移的雙頻疊加特性更為明顯，觸地點處的彎矩顯著增大。