天然氣凝析液長(zhǎng)距離管道穩(wěn)態(tài)水力熱力計(jì)算

王 智，宮 敬，尹鐵男，吳長(zhǎng)春

（中國(guó)石油大學(xué)（北京）城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 102249）

天然氣凝析液長(zhǎng)距離管道穩(wěn)態(tài)水力熱力計(jì)算

王 智，宮 敬，尹鐵男，吳長(zhǎng)春

（中國(guó)石油大學(xué)（北京）城市油氣輸配技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室，北京 102249）

針對(duì)長(zhǎng)距離天然氣凝析液管道輸送相間質(zhì)量、動(dòng)量和能量傳遞等問(wèn)題，基于多相流及化工熱力學(xué)理論，建立穩(wěn)態(tài)條件下水力熱力計(jì)算方法．結(jié)合兩相管流的機(jī)理模型與相平衡計(jì)算，采用簡(jiǎn)化的能量方程，耦合水力熱力，通過(guò)逐段循環(huán)迭代確定壓力、溫度、持液率等．通過(guò)閃蒸計(jì)算，得到指定溫度、壓力的氣、液組成．利用簡(jiǎn)化的Barnea流型判別方法判斷流型，由各流型下的流動(dòng)機(jī)理模型計(jì)算相應(yīng)的摩阻，持液率及氣、液相速度．管道模擬結(jié)果表明：文中求解方法穩(wěn)定可靠；精度與現(xiàn)場(chǎng)數(shù)據(jù)接近，能夠滿足大型管道工程需要．

兩相流；水力熱力計(jì)算；穩(wěn)態(tài)計(jì)算；天然氣凝析液；管道流動(dòng)

0 引言

常規(guī)油氣輸送在深海、沙漠腹地等地理環(huán)境復(fù)雜地區(qū)成本較高，因此經(jīng)濟(jì)性好的氣、液混輸技術(shù)得到日益廣泛的應(yīng)用．為了保證輸送安全運(yùn)行，管道設(shè)計(jì)合理，需要以準(zhǔn)確的管道流動(dòng)模擬作為基礎(chǔ)．人們主要以空氣、水體系為研究對(duì)象，而在混輸管道、石油化工、油路系統(tǒng)中大多接觸到的是油氣流體，其特點(diǎn)是組分多、相態(tài)變化復(fù)雜．對(duì)由相變導(dǎo)致的相間質(zhì)量、動(dòng)量和能量傳遞等問(wèn)題，天然氣凝析液管道的流動(dòng)參數(shù)計(jì)算尚不成熟．

早期的穩(wěn)態(tài)計(jì)算多采用相關(guān)式模型，常用的方法有Eaton相關(guān)式方法［1］、Oliemans相關(guān)式方法［2］及Beggs—Brill相關(guān)式方法［3］．當(dāng)管流條件適用時(shí)相關(guān)式模型比較準(zhǔn)確，但油氣兩相管流影響因素較多，純經(jīng)驗(yàn)相關(guān)式的計(jì)算精度和使用范圍受到限制，由此產(chǎn)生穩(wěn)態(tài)機(jī)理模型，如Wallis G B提出分散氣泡流的均質(zhì)無(wú)滑脫模型［4］；Taitel Y等提出分離流的機(jī)理模型［5］；Xiao J J等提出綜合機(jī)理模型，包括流型判定模型與各流型下的機(jī)理模型［6］．流型不同，氣、液兩相流動(dòng)的相互作用機(jī)理不同．Barnea D A綜合流型方面研究成果，提出可適用于各種傾角的流型判別模型［7］，該模型為目前計(jì)算結(jié)果較好的流型判別模型．Taitel Y等［8］認(rèn)為在石油工業(yè)中，瞬態(tài)相對(duì)較慢，應(yīng)該采用準(zhǔn)靜態(tài)動(dòng)量方程和瞬態(tài)連續(xù)方程得到簡(jiǎn)化的模型，并且設(shè)計(jì)相應(yīng)的求解方法．IFP、TOTAL和ELF AQUITAINE合作開(kāi)發(fā)基于漂移流模型的氣、液兩相管流瞬態(tài)模擬軟件TACITE．Tulsa大學(xué)多相流課題組也在無(wú)壓波模型的基礎(chǔ)上推出多相流軟件TUFFP．挪威IFE和SINTEF等機(jī)構(gòu)聯(lián)合開(kāi)發(fā)基于雙流體模型軟件OLGA，并得到行業(yè)內(nèi)認(rèn)可．

在穩(wěn)態(tài)流動(dòng)的前提條件下，筆者針對(duì)天然氣凝析液以兩相管流的機(jī)理模型為主體，結(jié)合相平衡計(jì)算、流型判斷模型及兩相管流相關(guān)式方法，耦合相態(tài)與水力熱力計(jì)算，建立天然氣凝析液管道穩(wěn)態(tài)條件下的水力熱力計(jì)算方法，為長(zhǎng)距離天然氣凝析液輸送管路的設(shè)計(jì)、計(jì)算提供簡(jiǎn)便易行的水力熱力計(jì)算方法．

1 控制方程及其簡(jiǎn)化

1．1 基本假設(shè)

（1）天然氣凝析液為氣、液兩相流動(dòng)狀態(tài)，忽略?shī)A帶；（2）一維流動(dòng)不考慮各參數(shù)在管道截面上分布的不均勻性，即截面上壓力、持液率、氣相速度、液相速度為截面平均值；（3）流動(dòng)是慢瞬變過(guò)程，穩(wěn)態(tài)條件下局部動(dòng)量平衡，忽略動(dòng)量方程中的時(shí)間項(xiàng)、加速度項(xiàng)及相間的壓差［8］．

1．2 連續(xù)性方程及動(dòng)量方程

基于質(zhì)量守恒和動(dòng)量守恒，在一維假設(shè)條件下，沿管長(zhǎng)得到氣、液相連續(xù)性方程組為

式中：Qg、Ql分別為氣、液相體積流量；ρg和ρl分別為氣、液相密度；A為管道截面積；L為管道長(zhǎng)度；˙mlg為單位管長(zhǎng)中液相的蒸發(fā)速率，反映液相向氣相傳遞的質(zhì)量［9］，

其中，G為總質(zhì)量流量下標(biāo)；X為氣相質(zhì)量相含率；p為壓力；T為管道中氣、液混合物的溫度；下標(biāo)T和p分別為溫度或壓力不變時(shí)質(zhì)量含氣率隨壓力和溫度的變化率．

氣、液相的聯(lián)合動(dòng)量方程［6］：

式中：Γgw和Γlw分別為氣、液相與管壁的剪切力；Ag和Al分別為氣、液相所占管道橫截面積；εg和εl分別為氣、液截面相分率；θ為管道傾角；Sg、Sl、Si分別為氣相、液相和相間濕周；fg和gl分別為氣、液相的范寧摩阻系數(shù)；ug和ul分別為氣、液相速度．

1．3 能量方程簡(jiǎn)化

在穩(wěn)態(tài)條件下，忽略瞬時(shí)項(xiàng)和加速度項(xiàng)的影響，能量方程［10］表示為式中，G為總質(zhì)量流量；下標(biāo)hm為單位質(zhì)量混合物的焓；cpm為氣、液混合物的定壓比熱容，cpm為氣、液混合物的湯焦效應(yīng)系數(shù)，為重力加速度；U為氣、液兩相沿著管壁向環(huán)境傳熱的總傳熱系數(shù)；Te為管道周?chē)h(huán)境溫度；D0為管道外徑．

2 封閉關(guān)系及流型判定

基于氣、液兩相流體模型建立控制方程組，摩阻封閉關(guān)系表示氣、液相與管壁間的剪切力，氣、液相間的剪切力，以及氣、液相相平衡關(guān)系．

2．1 組分熱力學(xué)模型

計(jì)算時(shí)天然氣凝析液的組成為已知條件．對(duì)于式（1）—（3）中所含有的流體物性參數(shù)，利用相對(duì)成熟的組分熱力學(xué)模型得到．采用PR狀態(tài)方程及范德華單流體混合規(guī)則，利用逸度模型，在壓力溫度一定的條件下進(jìn)行閃蒸的相平衡計(jì)算，得到氣、液兩相密度、焓、黏度及表面張力、湯焦效應(yīng)系數(shù)等熱力學(xué)參數(shù)［11］．

2．2 流型判定

流型判定是氣、液兩相管流穩(wěn)態(tài)計(jì)算的重要基礎(chǔ)，不僅直接影響相關(guān)式模型中壓降與持液率的計(jì)算，也是剪切力計(jì)算的重要依據(jù)．

在實(shí)際模擬程序中復(fù)雜的流型劃分并未提供給計(jì)算精度，各種商業(yè)軟件在模擬方法設(shè)計(jì)中都根據(jù)氣、液相間作用的本質(zhì)對(duì)流型進(jìn)行大幅簡(jiǎn)化［12］．Barnea D A建立適用于所有管徑、所有傾角的流型判別方法．根據(jù)Barnea流型判別準(zhǔn)則［13］簡(jiǎn)化流型判定［14］，針對(duì)天然氣凝析液管道中氣相占主導(dǎo)的特點(diǎn)，將流型劃分為泡狀流、分層流、環(huán)狀流和間歇流．

2．3 濕周與摩阻封閉關(guān)系

2．3．1 流型與濕周關(guān)系

根據(jù)文獻(xiàn)［15］，濕周與流型的關(guān)系見(jiàn)表1．其中：hl為液位高度；Lu為段塞單元長(zhǎng)度；Lf、Ls分別為液膜、液塞長(zhǎng)度；εg、εl分別為氣、液相的相含率；φ為液膜所對(duì)應(yīng)的半圓心角，

表1 流型濕周關(guān)系Table 1 Flow—pattern and wetted perimeter relationship

分別按泡狀流和分層流計(jì)算段塞流的液塞區(qū)與液膜區(qū)，然后基于各段長(zhǎng)度加權(quán)平均，求解速度界面處段塞單元的平均剪切力．液塞長(zhǎng)度采用斯科特（Scott）關(guān)系式［15］計(jì)算，液膜區(qū)泰勒氣泡的運(yùn)動(dòng)速度采用本迪克森（Bendiksen）關(guān)系式［15］計(jì)算．

2．3．2 范寧摩阻系數(shù)

（1）氣、液相與壁面：層流采用Hagen—Poiseille方程［16］計(jì)算；湍流采用Colebrook—White隱式方程［17］計(jì)算．

（2）氣、液相間：分層流采用Andritsos—Hanratty相關(guān)式［17］計(jì)算；環(huán)狀流采用Wallis［4］相關(guān)式計(jì)算．

泡狀流可視為均相流動(dòng)，兩相間沒(méi)有滑脫，因此不計(jì)相間阻力．對(duì)于氣、液相與壁面的范寧摩阻系數(shù)，在層流時(shí)采用Hagen—Poiseille方程計(jì)算；在湍流時(shí)采用Colebrook—White隱式方程計(jì)算，方程中相關(guān)參數(shù)使用氣、液混合物的平均值．

段塞流的范寧摩阻系數(shù)為液膜段和液塞段范寧摩阻在液塞單元長(zhǎng)度上的加權(quán)平均，液膜段按照分層流計(jì)算，液塞段按照泡狀流計(jì)算．

3 管段劃分及控制方程離散

首先將管道分為若干長(zhǎng)度管段，每個(gè)管段內(nèi)管道具有同一傾角和外部環(huán)境；然后在每一管段內(nèi)再將它等分為長(zhǎng)度近似的微元管段［18］．文中劃分約為1．0km的管段微元，最長(zhǎng)不超過(guò)1．5km．

任選第i個(gè)管段，對(duì)應(yīng)的上游節(jié)點(diǎn)為i—1和i（見(jiàn)圖1），得到2個(gè)節(jié)點(diǎn)的離散控制方程［11］．

氣相連續(xù)方程為

圖1 管道節(jié)點(diǎn)劃分示意Fig．1 Schematic diagram of discretization

液相連續(xù)方程為

聯(lián)合動(dòng)量方程為能量方程為

4 計(jì)算方法

針對(duì)穩(wěn)態(tài)問(wèn)題，管道各截面處的質(zhì)量流量相等，結(jié)合閃蒸計(jì)算從入口邊界推算下游相鄰節(jié)點(diǎn)的流速等變量．因此，采用向下游逐點(diǎn)求解的方法，從入口開(kāi)始向下游推算，一直達(dá)到管道出口［19］（見(jiàn)圖2）．

圖2 穩(wěn)態(tài)水力熱力計(jì)算流程Fig．2 Algorithm flow chart of Steady—state hydraulic thermodynamic

（1）將整條管道劃分為若干管段．

（2）在管段入口溫度壓力條件下進(jìn)行閃蒸，得到兩相的物性、相含率、折算速度等參數(shù)，根據(jù)物性參數(shù)，利用B—B［16］（Beggers—Brill）方法計(jì)算下一節(jié)點(diǎn)的各相流速、壓力初值．

（3）利用離散形式的簡(jiǎn)化能量方程求解下一節(jié)點(diǎn)溫度值．

（4）求解該管段的溫度、壓力平均值，并通過(guò)閃蒸計(jì)算得到該段的平均兩相的物性參數(shù)．

（5）根據(jù)速度初值，兩相物性等參數(shù)判斷流型．

（6）根據(jù)各流型下的機(jī)理模型計(jì)算該管段的剪切力、持液率及下一節(jié)點(diǎn)的壓力、折算速度．

（7）將上一步得到的壓力與該節(jié)點(diǎn)壓力初值比較，如滿足精度要求，返回步驟（2）計(jì)算下一個(gè)管段；反之，返回步驟（3），將得到的壓力帶入簡(jiǎn)化的能量方程，繼續(xù)循環(huán)，直至滿足精度要求．

對(duì)于步驟（6），在流型為分層流、段塞流、環(huán)狀流時(shí)，通過(guò)2個(gè)子步驟完成，首先求解各個(gè)流型下的平衡方程［16］；其次求解聯(lián)合動(dòng)量方程，得到下一節(jié)點(diǎn)的壓力．對(duì)于平衡方程，根據(jù)持液率、濕周、流型、流速及摩阻之間的關(guān)系，可以轉(zhuǎn)換為持液率的隱式方程，通過(guò)二分法或者牛頓法求解［20］．在泡狀流時(shí)，直接利用混合物的參數(shù)，按照單相流動(dòng)的方法求解壓降．

5 模擬計(jì)算

5．1 組分及基礎(chǔ)數(shù)據(jù)

為驗(yàn)證算法的適用性和精確程度，選擇5條天然氣凝析液管道進(jìn)行計(jì)算．管道為油氣田地面集輸管道，用P1—P5表示，均為起伏管道，傾角為0°～15°，全長(zhǎng)分別為74．0、50．0、43．4、6．8和4．3km，組分和基礎(chǔ)參數(shù)見(jiàn)表2和表3．

表2 P1—P5管道組分摩爾分?jǐn)?shù)Table 2 Mixture component in pipeline P1—P5 %

表3 P1—P5管道基礎(chǔ)參數(shù)Table 3 Basic parameter in pipeline P1—P5

5．2 實(shí)例驗(yàn)證

5條管道有實(shí)際運(yùn)行數(shù)據(jù)，為現(xiàn)場(chǎng)SCADA系統(tǒng)實(shí)時(shí)采集所得到的數(shù)據(jù)平均值．對(duì)比現(xiàn)場(chǎng)與計(jì)算數(shù)據(jù)可以看出，計(jì)算值和實(shí)測(cè)值絕對(duì)誤差上，壓降相差0．3～0．8MPa，溫度相差2～7℃不等；相對(duì)誤差上，壓降為3%～11%，溫降為7%～14%（見(jiàn)圖4）．產(chǎn)生差異的原因：對(duì)于大型管道，基礎(chǔ)參數(shù)相對(duì)宏觀和粗略，如P2管道實(shí)測(cè)入口溫度為21．7℃，出口溫度為29．6℃，而提供的環(huán)境溫度為28℃；另外，文中采用各類(lèi)文獻(xiàn)中較認(rèn)可的小口徑多相實(shí)驗(yàn)環(huán)道的經(jīng)驗(yàn)關(guān)系，與實(shí)際大口徑管道相比存在一定偏差［21］．總體上，對(duì)于大型天然氣凝析液管道工程，精度可以滿足實(shí)踐需要．

表4 現(xiàn)場(chǎng)與計(jì)算數(shù)據(jù)對(duì)比結(jié)果Table 4 Comparison result between calculation and field

根據(jù)壓力、溫度及持液率的沿程分布，以P1管道為例，計(jì)算結(jié)果表明，流體的溫度隨著輸送距離延伸不斷下降，在管道的上游下降較快，隨著管道延伸溫降趨于緩和（見(jiàn)圖3）．這是流體不斷與周?chē)橘|(zhì)換熱造成的，流體最初進(jìn)入管道時(shí)流體溫度與周?chē)h(huán)境溫度相差較大，換熱相對(duì)劇烈，因此出現(xiàn)較大的溫降；隨著管道流動(dòng)距離的延伸，流體與環(huán)境的溫差不斷減小，換熱趨于緩和，因此管道下游的溫降較上游的要小．

根據(jù)管道溫度壓力分布，由相平衡理論，整條管道的壓力溫度分布應(yīng)處于相圖中的反凝析區(qū)，因此無(wú)論壓力還是溫度降低都使液相含量增加．根據(jù)持液率分布，盡管局部受到地形影響（表現(xiàn)為上坡段液相積累，下坡段液相加速）（見(jiàn)圖3），總體上持液率增加．另外，由于溫降總體上比較平穩(wěn)，所以液體的析出速率總體變化不大；持液率的變化更多地受到地形影響．

圖3 P1管道流動(dòng)參數(shù)沿程分布Fig．3 Flow parameter distribution in pipeline P1

6 結(jié)束語(yǔ)

在多相流動(dòng)及熱力學(xué)理論基礎(chǔ)上，提出一種穩(wěn)態(tài)條件下的天然氣凝析液長(zhǎng)距離管道輸送的水力熱力計(jì)算方法，引入簡(jiǎn)化的Barnea流型判別方法，使各種流型下的機(jī)理模型統(tǒng)一于流型判斷中；在兩相管流機(jī)理模型基礎(chǔ)上，結(jié)合熱力計(jì)算及相平衡計(jì)算，采用相關(guān)式方法的部分成果為計(jì)算提供初值和部分經(jīng)驗(yàn)關(guān)系．

在對(duì)天然氣凝析液混輸管道進(jìn)行水力熱力計(jì)算時(shí)，計(jì)算結(jié)果在各參數(shù)的沿線趨勢(shì)上符合基本物理規(guī)律，數(shù)值與現(xiàn)場(chǎng)實(shí)時(shí)數(shù)據(jù)的平均值接近．總體上，該算法簡(jiǎn)便易行，能夠滿足大型管道工程精度要求．

［1］ Eaton B A，Andrews D E，Knowles C R，et al．The prediction of flow patterns，liquid holdup and pressure losses occurring during continuous two—phase in horizontal pipelines［J］．J．P．T．，1967，（6）：815—923．

［2］ Oliemans R V A，Pots B F M，Trompe N．Modeling of annular dispersed two—phase flow in vertical pipes［J］．Intl．J．Multiphase Flow，1986，12（5）：711—732．

［3］ Beggs H D，Brill J P．A study of two phase flow in inclined pipes［J］．J．P．T．，1973（5）：607—617．

［4］ Wallis G B．One dimension two—phase flow［M］．New York：McGraw—hill，1969．

［5］ Taitel Y，Lee N，Dukler A E．Transient gas—liquid flow in horizontal pipes：modeling the flow pattern transitions［J］．AIChE Jour—nal，1978，24（5）：920—934．

［6］ Xiao J J，Shoham O，Brill J P．A comprehensive mechanistic model for two—phase flow in pipelines［C］．New Orleans：The 65th An—nual Technical Conference and Exhibition of the Society of Engineers，SPE 20631，1990．

［7］ Barnea D A．Unified model for predicting flow pattern transitions for the whole range of pipe inclinations［J］．Intl．J．Multiphase Flow，1987，13（1）：1—12．

［8］ Taitel Y，Barnea D．Simplified transient simulation of two phase flow using quasi—equilibrium momentum balance［J］．Int．J．Multi—phase Flow 1997，23（3）：493—501．

［9］ 李玉星．濕天然氣管道輸送工藝計(jì)算研究［D］．北京：中國(guó)石油大學(xué)（北京），1997：9—10．LI Yuxing．Study on wet gas transportation process calculation［D］．Beijing：China University of Petroleum（Beijing），1997：9—10．

［10］ 鄧道明．天然氣——凝析液混輸管流特性模擬研究［D］．北京：中國(guó)石油大學(xué)（北京），2005：50—82．Deng Daoming．Modeling gas—condensate two—phase flow in pipeline［D］．Beijing：China University of Petroleum（Beijing）：2005：50—82．

［11］ 黃小紅，左龍．基于計(jì)算方法的地層油注氣后相態(tài)特征［J］．大慶石油學(xué)院學(xué)報(bào)，2006，30（3）：54—56．Huang Xiaohong，Zuo Long．Phase state characteristic after gas injection into formation oil based on computational method［J］．Journal of Daqing Petroleum Institute，2006，30（3）：54—56．

［12］ Dhulesia H，Lopez D，Aquitaine E．Critical evaluation of mechanistic two—phase flow pipeline and well simulation models［C］．SPE 36611，1996．

［13］ Barnea D．A unified model for predicting flow pattern transitions for the whole range of pipe inclinations［J］．Intl．J．Multiphase Flow，1987，13（1）：1—12．

［14］ 韓洪生，王忠信，楊樹(shù)人，等．圓管中氣、液相流動(dòng)空隙率數(shù)學(xué)模型［J］．大慶石油學(xué)院學(xué)報(bào)，2002，26（4）：19—21．Han Hongsheng，Wang Zhongxin，Yang shuren，et al．Mathematical model of void fraction for gas—liquid two—phase flow in pipelines ［J］．Journal of Daqing Petroleum Institute，2002，26（4）：19—21．

［15］ Ovadia Shoham．Mechanistic modeling of gas—liquid two—phase flow in pipes［M］．USA：Society of Petroleum Engineers，2006．

［16］ 陳濤平，張權(quán)．高凝油抽油機(jī)井井筒壓力——溫度分布預(yù)測(cè)［J］．大慶石油學(xué)院學(xué)報(bào)，2011，35（2）：62—67．Chen Taoping，Zhang Quan．Prediction of wellbore pressure and temperature distribution in the pumping wells producing high pour point oil［J］．Journal of Daqing Petroleum Institute，2011，35（2）：62—67．

［17］ Andritsos N，Hanratty T J．Influence of interfacial waves in stratified gas—liquid flows［J］．AIChE J，1987，33（3）：444—454．

［18］ 吳國(guó)忠，陳超．埋地管道傳熱數(shù)值模擬網(wǎng)格劃分方法［J］．大慶石油學(xué)院學(xué)報(bào)，2005，29（2）：82—84．Wu Gongzhong，Chen Chao．Grid division method for buried pipelines during numerical simulation［J］．Journal of Daqing Petroleum Institute，2005，29（2）：82—84．

［19］ 賈永英，劉揚(yáng)，劉曉燕，等．含蠟原油管道再啟動(dòng)壓力計(jì)算［J］．大慶石油學(xué)院學(xué)報(bào)，2008，32（5）：45—54．Jia Yongying，Liu Yang，Liu Xiaoyan，et al．Pressure calculation of the restart of waxy crude oil pipelines［J］．Journal of Daqing Pe—troleum Institute，2008，32（5）：45—47．

［20］ 賈永英，劉揚(yáng)，王玉潔，等．非牛頓原油層流埋地?zé)嵊凸艿姥鼐€溫度的確定［J］．大慶石油學(xué)院學(xué)報(bào)，2006，30（6）：51—52．Jia Yongying，Liu Yang，Wang Yujie，et al．Temperature determination of laminar of non—Newtonian crude oil in buried pipelines ［J］．Journal of Daqing Petroleum Institute，2006，30（6）：51—52．

［21］ Bendiksen K H，Malnes D，Moe R，et al．The dynamic two—fluid model OLGA：theory and application［C］．SPE 19451，1991．

TE832

A

2095—4107（2013）04—0060—07

DOI 10．3969／j．issn．2095—4107．2013．04．009

2013—01—17；編輯：張兆虹

國(guó)家科技重大專(zhuān)項(xiàng)（2011ZX05026—004—03）

王 智（1982—），男，博士研究生，主要從事多相管流及油氣田集輸技術(shù)方面的研究．