顧繼俊 李云龍

(中國石油大學(北京))

0 引 言

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

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

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

1 模型建立

1.1 幾何模型建立

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

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

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

1.2 邊界條件設置

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

1.3 湍流模型

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

(1)

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

(2)

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

湍流動能方程為：

(3)

湍流耗散率方程為：

(4)

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

上述模型經(jīng)驗系數(shù)均來自于對前人試驗的總結，可以在大部分工況下使用。

1.4 工況參數(shù)設置

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

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

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

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

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

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

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

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

圖3 計算收斂圖Fig.3 Calculation convergence diagram

圖4 兩套模型在P2點Y向位移對比Fig.4 Y-direction displacement comparison of P2 between two models

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

2.1 下傾管內氣液相分布

分別模擬入口氣相體積分數(shù)為0.75、0.50和0.25的工況，具體參數(shù)如表2中的工況1、2、3所示。模擬得到距離下傾管入口3、6和9 m處截面的氣液相分布，如圖5～圖7所示，分析了不同入口氣相體積分數(shù)對下傾管內部兩相流流型的影響。

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

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

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

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

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

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

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

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

圖12 P2監(jiān)測點所受應力Fig.12 Stress on monitoring point P2

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

2.3 入口氣相體積分數(shù)對管道振動響應的影響

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

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

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

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

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

圖17 工況2立管應力云圖Fig.17 Stress cloud of riser in working condition 2

圖18 工況3立管應力云圖Fig.18 Stress cloud of riser in working condition 3

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

3 結 論

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

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

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

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

(4)當入口氣相體積分數(shù)從0.50增加到0.75時，立管的振動位移和最大應力均發(fā)生大幅度增長，這是因為立管系統(tǒng)響應頻率與其一階固有頻率吻合，導致立管發(fā)生了共振現(xiàn)象。