顧繼俊 李云龍

(中國(guó)石油大學(xué)(北京))

0 引 言

在深水油氣田勘探開(kāi)發(fā)中，深水立管是必不可少的水下設(shè)備[1]。例如常見(jiàn)的鋼制懸鏈線(xiàn)立管(SCR)可應(yīng)用在水深3 000 m左右的海況中，但因其為鋼制，導(dǎo)致立管觸地區(qū)域(TDP)彎曲剛度非常大[2]，容易產(chǎn)生疲勞損傷。為了解決鋼制懸鏈線(xiàn)立管在深水中存在的問(wèn)題，一種深水緩波形立管(LWR)以其更好的順應(yīng)性得到廣泛應(yīng)用。其原理是在立管上安裝質(zhì)量較輕的浮筒，利用浮力使立管的一段隆起，形成類(lèi)似于波浪的形狀，從而使立管頂部與觸地點(diǎn)分離[3]，以減小外部載荷對(duì)觸地點(diǎn)的影響。

深水油氣田開(kāi)發(fā)及混輸工藝要求立管可以輸送較大流量的氣液混合物[4-5]。在輸送過(guò)程中氣液混合物與立管的流固耦合作用機(jī)理很復(fù)雜，可能會(huì)造成立管的參數(shù)共振或組合共振，這使得流體引起的立管動(dòng)力響應(yīng)問(wèn)題很難預(yù)測(cè)。如果長(zhǎng)時(shí)間振動(dòng)容易造成立管疲勞損傷甚至發(fā)生破裂，導(dǎo)致油氣泄漏，從而給國(guó)家經(jīng)濟(jì)和海洋資源帶來(lái)嚴(yán)重影響[6]。因此，研究氣液兩相流作用下緩波形立管的動(dòng)力響應(yīng)問(wèn)題對(duì)保護(hù)海洋立管系統(tǒng)具有重要意義。J.L.RIVERIN等[7]對(duì)復(fù)雜結(jié)構(gòu)管道進(jìn)行了試驗(yàn)，認(rèn)為雖然內(nèi)流對(duì)管道沖擊力很小，但會(huì)導(dǎo)致其劇烈振動(dòng)。XING L.C.等[8]、YE J.等[9]、LUO X.M.等[10]和LI N.L.等[11]分別對(duì)不同結(jié)構(gòu)管道的段塞流現(xiàn)象展開(kāi)試驗(yàn)，并對(duì)其持液率、周期和壓力等參數(shù)進(jìn)行了分析。B.I.STUART[12]基于能量守恒定理，建立了含有油氣水三相流垂直上升管的壓力梯度表達(dá)式。劉昶等[13]為了識(shí)別海洋立管系統(tǒng)中可能出現(xiàn)的流型及考察立管系統(tǒng)在各流型作用下的動(dòng)態(tài)響應(yīng),利用室內(nèi)試驗(yàn)裝置模擬海洋立管系統(tǒng),對(duì)立管中兩相流動(dòng)壓力波動(dòng)及其引起的管道振動(dòng)位移進(jìn)行了測(cè)試。謝翠麗等[14]為研究氣液兩相混輸管道內(nèi)流導(dǎo)致管道破壞的機(jī)理，利用 ANSYS Workbench 軟件，建立了L形彎管有限元模型和氣液兩相CFD模型，進(jìn)行了L形彎管內(nèi)流致振的流固耦合數(shù)值模擬。楊思齊等[15]對(duì)高壓管匯彎管部位壓裂液固液兩相流動(dòng)及流固耦合效應(yīng)進(jìn)行了研究，認(rèn)為隨支撐劑質(zhì)量濃度和流速增大，流固耦合效果越來(lái)越明顯，彎管剪應(yīng)力和變形程度增加，其中流速能顯著影響彎管流固耦合程度。馬天麒等[16]對(duì)內(nèi)輸多相流立管在內(nèi)流與外部繞流耦合作用下的振動(dòng)響應(yīng)進(jìn)行分析，建立了內(nèi)部多相流-立管-外部繞流的耦合方程,并采用廣義積分變換法將偏微分方程轉(zhuǎn)化為常微分方程進(jìn)行求解。

以上研究均未考慮緩波形立管中氣液兩相流引起的立管流致振動(dòng)問(wèn)題。為此，筆者建立了緩波形立管的兩相流模型，在不同工況下進(jìn)行了數(shù)值模擬分析，模擬結(jié)果可為緩波形立管的振動(dòng)和疲勞特性評(píng)估提供指導(dǎo)。

1 模型建立

1.1 幾何模型建立

緩波形立管主要包括懸垂段、浮子段和下降段。為了保證流體流動(dòng)的完整性，又添加了連接井口的下傾管，它與前面3段共同組成了局部緩波形立管系統(tǒng)模型，具體尺寸如表1所示。模型尺寸均根據(jù)相似理論計(jì)算所得，保證了模型的合理性。立管模型材料為PVC軟管，其參數(shù)為：彎曲剛度36 N·m2，密度1 570 kg/m3，彈性模量3.2 GPa。下傾管的左邊界設(shè)置為兩相流質(zhì)量流量入口，懸垂段的上邊界設(shè)置為壓力出口。為了研究?jī)上嗔髁⒐艿膭?dòng)力響應(yīng)問(wèn)題，分別選取P1、P2、P3和P4等4個(gè)關(guān)鍵點(diǎn)進(jìn)行監(jiān)測(cè)。緩波形立管模型如圖1所示。

表1 緩波形立管模型尺寸 mmTable 1 Dimension of lazy-wave riser mm

圖1 緩波形立管模型Fig.1 Model of lazy-wave riser

1.2 邊界條件設(shè)置

將下傾管左邊界端面設(shè)為固支，懸垂段上邊界端面設(shè)為簡(jiǎn)支，立管內(nèi)壁面設(shè)為流固耦合面，用來(lái)傳遞計(jì)算數(shù)據(jù)。模擬總時(shí)長(zhǎng)根據(jù)實(shí)際情況設(shè)定，時(shí)間步長(zhǎng)設(shè)為0.05 s。當(dāng)氣液相流速較慢時(shí)，為保證流動(dòng)的完整性，將會(huì)適當(dāng)延長(zhǎng)模擬總時(shí)長(zhǎng)，時(shí)間步長(zhǎng)保持不變。

1.3 湍流模型

有限元模擬最常用的模型是湍流模型。湍流模型可細(xì)分為雷諾應(yīng)力模型和渦黏模型，工程應(yīng)用最普遍的是渦黏模型。通過(guò)引入湍流黏度，利用Boussinesq假設(shè)，將雷諾應(yīng)力與湍流黏度聯(lián)系起來(lái)，其關(guān)系式為：

(1)

式中：σij為Kronecker函數(shù)，湍流黏度μt通常使用k-ε模型、k-ω模型和Spalart-Allmaras模型，u為湍流流體的速度，ρ為流體的密度。

(2)

式中：k為湍流動(dòng)能；ε為湍流耗散率；Cμ為經(jīng)驗(yàn)系數(shù)，通常取0.09。其中湍流動(dòng)能以及湍流耗散率可以用以下方程進(jìn)行求解。

湍流動(dòng)能方程為：

(3)

湍流耗散率方程為：

(4)

式中：Gk為層流速度梯度而產(chǎn)生的湍流動(dòng)能；Gb為浮力產(chǎn)生的湍流動(dòng)能；YM為過(guò)度擴(kuò)散產(chǎn)生的波動(dòng)；C1ε、C2ε、C3ε為經(jīng)驗(yàn)系數(shù)，分別取1.44、1.92和0；σk為k方程的湍流普朗特?cái)?shù)，取1.0；σε為ε方程的湍流普朗特?cái)?shù)，取1.0；Sk、Sε為自定義。

上述模型經(jīng)驗(yàn)系數(shù)均來(lái)自于對(duì)前人試驗(yàn)的總結(jié)，可以在大部分工況下使用。

1.4 工況參數(shù)設(shè)置

從文獻(xiàn)[17]中選取出部分試驗(yàn)工況，并結(jié)合實(shí)際生產(chǎn)情況設(shè)置了4組模擬工況，如表2所示。

表2 工況參數(shù)Table 2 Parameters of working conditions

1.5 網(wǎng)格收斂性分析

網(wǎng)格劃分時(shí)，將網(wǎng)格尺寸控制在小于水泡直徑3.5 mm。立管壁面邊界層設(shè)為5層，總網(wǎng)格數(shù)為188 541，如圖2所示。劃分的網(wǎng)格質(zhì)量均大于0.4，其中99.74%以上的網(wǎng)格質(zhì)量在0.93以上，達(dá)到模擬要求。

圖2 有限元模型Fig.2 Finite element model

繼續(xù)將網(wǎng)格加密，將立管壁面邊界層設(shè)置為10層，并對(duì)立管彎曲處進(jìn)行網(wǎng)格加密處理，得出新的數(shù)值模型的網(wǎng)格總數(shù)量為537 482，網(wǎng)格質(zhì)量均大于0.6。利用兩套數(shù)值模型分別針對(duì)工況1進(jìn)行數(shù)值模擬，通過(guò)監(jiān)測(cè)同一位置的振動(dòng)位移，分析網(wǎng)格質(zhì)量對(duì)計(jì)算精度的影響，從而選擇最佳的數(shù)值模型進(jìn)行剩余工況的計(jì)算。

目前，我國(guó)職業(yè)經(jīng)理人市場(chǎng)存在的問(wèn)題是：職業(yè)化程度低、數(shù)量少、法律法規(guī)體系不健全、企業(yè)主和經(jīng)理人的信任體系沒(méi)有建立起來(lái)、認(rèn)證和評(píng)價(jià)體系還沒(méi)有各方的認(rèn)可。這就導(dǎo)致企業(yè)在選聘職業(yè)經(jīng)理人時(shí)獲取信息存在問(wèn)題。

圖3為工況1下未加密的數(shù)值模型在計(jì)算過(guò)程中的迭代收斂情況。圖3縱坐標(biāo)表示以10為底的對(duì)數(shù)收斂范數(shù)。由圖3可以看出迭代至200步左右收斂。圖4為兩套數(shù)值模型在P2監(jiān)測(cè)點(diǎn)處的Y向位移對(duì)比。從圖4可以看出，兩套數(shù)值模型計(jì)算結(jié)果差異很小，說(shuō)明未加密的數(shù)值模型計(jì)算精度也比較高。因此，為了提高計(jì)算效率，使用未加密的數(shù)值模型計(jì)算剩余工況。

圖3 計(jì)算收斂圖Fig.3 Calculation convergence diagram

圖4 兩套模型在P2點(diǎn)Y向位移對(duì)比Fig.4 Y-direction displacement comparison of P2 between two models

2 氣液兩相流緩波形立管數(shù)值模擬

2.1 下傾管內(nèi)氣液相分布

分別模擬入口氣相體積分?jǐn)?shù)為0.75、0.50和0.25的工況，具體參數(shù)如表2中的工況1、2、3所示。模擬得到距離下傾管入口3、6和9 m處截面的氣液相分布，如圖5～圖7所示，分析了不同入口氣相體積分?jǐn)?shù)對(duì)下傾管內(nèi)部?jī)上嗔髁餍偷挠绊憽?/p>

圖5 工況1不同位置截面氣液相分布Fig.5 Gas and liquid distribution at different cross sections under working condition 1

從圖5～圖7可以看出，隨著入口氣相體積分?jǐn)?shù)的變化，管內(nèi)流型也發(fā)生變化。其中距離流體入口越遠(yuǎn)的截面其含氣體積分?jǐn)?shù)越高，但總體來(lái)看，同一工況中下傾管各截面的氣液體積分?jǐn)?shù)波動(dòng)不明顯。工況1下傾管中的液相體積分?jǐn)?shù)都在0.748 9～0.751 7之間，說(shuō)明整個(gè)管內(nèi)氣液相都處于混合狀態(tài)。這是由于工況1下氣液相的入口流速懸殊較大，加之液相體積分?jǐn)?shù)比較低，隨著管內(nèi)氣體高速流過(guò)，帶動(dòng)液相從管道底部脫離，以小液滴的形式與氣體混合向前流動(dòng)，導(dǎo)致下傾管中兩相流呈現(xiàn)混合流動(dòng)狀態(tài)。工況2下管內(nèi)兩相流開(kāi)始出現(xiàn)分層流動(dòng)狀態(tài)，工況3下管內(nèi)的兩相流分層流動(dòng)現(xiàn)象更加明顯。因此，隨著入口氣相體積分?jǐn)?shù)的減小，下傾管中更容易出現(xiàn)分層流動(dòng)現(xiàn)象。

圖6 工況2不同位置截面氣液相分布Fig.6 Gas and liquid distribution at different cross sections under working condition 2

圖7 工況3不同位置截面氣液相分布Fig.7 Gas and liquid distribution at different cross sections under working condition 3

2.2 緩波形立管中的段塞流

通過(guò)對(duì)工況4數(shù)值模擬，得到緩波形立管中的流體流動(dòng)狀態(tài)，如圖8所示。從圖8可以看出，緩波形立管中已經(jīng)發(fā)生段塞流現(xiàn)象，下傾管中的液相流至緩波形立管下降段底部時(shí)，由于慣性力作用液體將快速形成液塞堵住下降段底部，如圖8中的細(xì)節(jié)圖(2)所示。因?yàn)樗容^容易變形，在重力作用下導(dǎo)致液塞回落到下降段底部，此時(shí)氣體將再次進(jìn)入下降段。隨著液體的連續(xù)涌入，下降段底部的液面逐漸升高，促使液塞再次形成。接下來(lái)液塞又會(huì)在重力和氣體沖擊作用下被破壞，如此反復(fù)，直到下降段底部的液體逐漸升高形成穩(wěn)定的液柱，如圖8中細(xì)節(jié)圖(3)所示。此時(shí)，緩波形立管的懸垂段也會(huì)因?yàn)樯鲜銮闆r形成穩(wěn)定的液柱，在浮子段形成氣塞，如圖8中的細(xì)節(jié)圖(4)所示。由于緩波形立管中液柱堵住了氣體的通路，所以造成下傾管中幾乎被氣體充滿(mǎn)，如圖8中的細(xì)節(jié)圖(1)所示，整個(gè)緩波形立管系統(tǒng)就形成了兩段比較嚴(yán)重的段塞流。

圖8 工況4緩波形立管內(nèi)氣液分布Fig.8 Gas and liquid distribution in the lazy-wave riser under working condition 4

由模擬結(jié)果可知，立管系統(tǒng)中最大振動(dòng)位移出現(xiàn)在緩波形立管下降段的底部，因此選擇P2為監(jiān)測(cè)點(diǎn)。工況4下，該點(diǎn)在X、Y、Z方向的振動(dòng)位移和應(yīng)力如圖9～圖12所示。

圖9 P2監(jiān)測(cè)點(diǎn)在X方向的位移Fig.9 Displacement of monitoring point P2 in X direction

圖10 P2監(jiān)測(cè)點(diǎn)在Y方向的位移Fig.10 Displacement of monitoring point P2 in Y direction

圖11 P2監(jiān)測(cè)點(diǎn)在Z方向的位移Fig.11 Displacement of monitoring point P2 in Z direction

圖12 P2監(jiān)測(cè)點(diǎn)所受應(yīng)力Fig.12 Stress on monitoring point P2

綜上所述，工況4條件下立管系統(tǒng)產(chǎn)生比較嚴(yán)重的段塞流流型。由圖9～圖11可以看出，P2監(jiān)測(cè)點(diǎn)在Y方向的振動(dòng)位移遠(yuǎn)大于在X、Z方向的振動(dòng)位移。當(dāng)管內(nèi)兩相流流動(dòng)至142 s左右時(shí)，P2點(diǎn)在各方向上的位移均發(fā)生了不同幅度的波動(dòng)。其中，X方向最大位移從2.60×10-6m增長(zhǎng)到3.42×10-6m，增幅達(dá)到16%。Y方向最大位移從4.2×10-4m增長(zhǎng)到4.7×10-4m，增幅達(dá)到12%。Z方向最大位移從2.6×10-7m增長(zhǎng)到4.55×10-7m，增幅達(dá)到75%。由圖12可以看出，P2監(jiān)測(cè)點(diǎn)應(yīng)力值從400 kPa增長(zhǎng)到460 kPa，增幅達(dá)到15%。這是因?yàn)榱黧w流動(dòng)至142 s左右，段塞流開(kāi)始噴發(fā)，致使立管發(fā)生明顯振動(dòng)。由此可見(jiàn)，段塞流的噴發(fā)會(huì)對(duì)立管振動(dòng)產(chǎn)生比較大的影響。

2.3 入口氣相體積分?jǐn)?shù)對(duì)管道振動(dòng)響應(yīng)的影響

通過(guò)對(duì)比工況1、2、3的模擬結(jié)果發(fā)現(xiàn)，緩波形立管最大振動(dòng)位移同樣出現(xiàn)在P2監(jiān)測(cè)點(diǎn)，其在X、Y、Z方向上的位移時(shí)程曲線(xiàn)如圖13～圖15所示。從圖13～圖15可以看出：在流動(dòng)的初期3種工況下，管內(nèi)流體在X、Y、Z方向均產(chǎn)生較大的振動(dòng)位移，其中Y方向的位移遠(yuǎn)大于X和Z方向，這是因?yàn)楣軆?nèi)兩相流流動(dòng)初期不穩(wěn)定；隨后各工況的振動(dòng)位移逐漸發(fā)生衰減直至穩(wěn)定振動(dòng)，此時(shí)Y方向的振幅均大于X和Z方向，且振動(dòng)頻率基本保持不變，其中工況1下的振幅均遠(yuǎn)大于工況2和工況3。

圖13 P2監(jiān)測(cè)點(diǎn)在X方向的位移Fig.13 Displacement of monitoring point P2 in X direction

圖14 P2監(jiān)測(cè)點(diǎn)在Y方向的位移Fig.14 Displacement of monitoring point P2 in Y direction

圖15 P2監(jiān)測(cè)點(diǎn)在Z方向的位移Fig.15 Displacement of monitoring point P2 in Z direction

工況1、2、3下立管系統(tǒng)所受應(yīng)力分布如圖16～圖18所示。由圖16～圖18可以發(fā)現(xiàn)：3種工況下立管系統(tǒng)的最大應(yīng)力均出現(xiàn)在流體入口的固支端；工況1的最大應(yīng)力為11.250 MPa，工況2的最大應(yīng)力為6.679 MPa，工況3的最大應(yīng)力值為4.629 MPa。對(duì)比3種工況下立管所受應(yīng)力極值可以發(fā)現(xiàn)，隨著氣相體積分?jǐn)?shù)的增大，立管所受應(yīng)力極值持續(xù)增大，而當(dāng)入口氣相體積分?jǐn)?shù)從0.50增加到0.75時(shí)，P2監(jiān)測(cè)點(diǎn)所受的應(yīng)力極值有大幅增大。

圖17 工況2立管應(yīng)力云圖Fig.17 Stress cloud of riser in working condition 2

圖18 工況3立管應(yīng)力云圖Fig.18 Stress cloud of riser in working condition 3

綜上所述，工況1下立管的振幅和應(yīng)力極值較其他兩種工況均有大幅增長(zhǎng)。其主要原因是當(dāng)入口氣相體積分?jǐn)?shù)為0.75時(shí)，立管系統(tǒng)的響應(yīng)頻率與其一階固有頻率吻合，立管發(fā)生了共振現(xiàn)象，同時(shí)也表明改變兩相流的流動(dòng)參數(shù)可以引發(fā)立管共振。

3 結(jié) 論

本文建立了氣液兩相流緩波形立管的雙向流固耦合模型，應(yīng)用k-ε湍流模型，通過(guò)設(shè)置不同入口氣相體積分?jǐn)?shù)和入口質(zhì)量流量模擬了4種工況，分析了不同的入口邊界條件對(duì)緩波形立管系統(tǒng)中流型變化和管道振動(dòng)響應(yīng)的影響，得出以下結(jié)論。

(1)工況1中下傾管內(nèi)兩相流流型為氣液混合流動(dòng)，工況2和工況3為氣液分層流動(dòng)，隨著入口氣相體積分?jǐn)?shù)的減小，下傾管中的分層流動(dòng)現(xiàn)象更明顯。

(2)緩波形立管中出現(xiàn)了兩段比較嚴(yán)重的段塞流。段塞流的噴發(fā)加劇了立管的振動(dòng)，入口氣液相流量較小時(shí)更容易形成嚴(yán)重段塞流。

(3)工況1、2、3下Y向振幅遠(yuǎn)大于X和Z方向振幅，其中工況1導(dǎo)致立管的振動(dòng)強(qiáng)度最大。3種工況下立管所受的最大應(yīng)力均處在流體入口的固支端。

(4)當(dāng)入口氣相體積分?jǐn)?shù)從0.50增加到0.75時(shí)，立管的振動(dòng)位移和最大應(yīng)力均發(fā)生大幅度增長(zhǎng)，這是因?yàn)榱⒐芟到y(tǒng)響應(yīng)頻率與其一階固有頻率吻合，導(dǎo)致立管發(fā)生了共振現(xiàn)象。