濕氣管道積液臨界氣速預測的新模型

李國豪，鄧道明，宮敬

（中國石油大學（北京）油氣管道輸送安全國家工程實驗室，石油工程教育部重點實驗室，城市油氣輸配技術北京市重點實驗室，北京102249）

引 言

國內外天然氣集輸經常采用濕氣集輸的方式。Meng[1]指出標況（101.3 kPa, 15℃）下液相體積流量與氣相體積流量之比小于1100 m3/106m3時為低液相負荷流動。濕氣管道一般為低液相負荷流動，例如，我國某頁巖氣田濕氣管道氣液分離之前的平均液相負荷一般為200 m3/106m3左右，最大的平均液相負荷為1000 m3/106m3。但是即便液相負荷如此低和井場集氣站進行了氣液分離，該頁巖氣田的一些氣井和早期集輸管線仍存在較嚴重積液問題。

早期的一些學者認為，低液相負荷下天然氣兩相管流為彌散霧狀流，即認為低液相負荷流動不會產生嚴重積液現象。例如Adewumi 等[2]在美國國家自然科學基金、工程基金和氣體研究所(GRI)的資助下，對低液相負荷天然氣-凝析液管流進行了一系列研究，其模擬計算理論主要基于漂移的彌散霧狀流模型。但是其模型計算結果沒有得到實驗或生產數據的廣泛支持。顯然濕氣管道的彌散霧狀流模型難以解釋經過分離器分離的濕氣在下游集輸管道里仍然積液嚴重的原因。

正常情況下，濕氣管道內低液相負荷的流動為氣液兩相分層流。當氣田開發初期開井數較少或氣田開發后期氣井產量較低時，集輸管內氣體流速較低，如果氣相對液相的拖曳力不足以拖動液相向前流動，氣液分層流便開始變得不穩定，最終液體會在上傾管道低洼處積集，積液現象因此產生。另外過大的集輸管道直徑也會引起積液。本文將上傾管道液膜不穩定的最大氣相表觀速度（或者說液膜穩定的最小氣相表觀速度）作為臨界氣速。工程實踐表明，起伏的濕氣集輸管道積液是普遍現象。

濕氣集輸管道積液的存在會導致諸多風險：（1）誘使管內產生內腐蝕，進而影響生產運行[3]；（2）集輸管道積液意味著存在更大的相間滑脫損失，這會導致井口回壓增加，從而降低氣井產量，甚至可能造成氣井積液；（3）集輸管道積液還伴隨著管道集輸效率降低、超壓安全風險、下游處理量不穩定等問題。因此，研究濕氣管道的積液對集輸管道的設計和運行都具有重要的指導意義。

1 積液預測模型評述

一些學者采用商業多相流軟件OLGA 對濕氣管道進行積液研究。例如劉建武等[4]、張愛娟等[5]建立起伏的集輸管道的OLGA 計算模型，經過OLGA 計算獲得了管道全線的積液量等相對宏觀的參數。但是，OLGA 軟件作為一個黑箱，使用者并不知曉其具體的計算模型和方法，難以利用OLGA 建立積液機理計算模型。

目前主流的積液機理模型有最小壓降模型、持液率突變模型、零液壁剪切應力模型和最小界面剪切應力模型4種。

1.1 最小壓降模型

最小壓降模型是指：逐漸減小氣相流量，當氣液兩相管流的總壓降最小時，管內由摩阻壓降主導轉換為重力壓降主導，管內發生積液。Belt[6]、Yuan等[7]均支持該模型。Fan 等[8]的實驗也表明，當水平傾角大于5°時，最小壓降與臨界氣速相對應。

但是，最小壓降并不能完全解釋積液的機理。盡管臨界氣速測量方法相同，Alsaadi 等[9]和Wang等[10]的實驗結果卻截然不同：Alsaadi等[9]的實驗管徑是76.2 mm，在實驗中觀察到最小壓降對應的氣速要小于臨界氣速；Wang 等[10]的實驗管徑是30 mm 和40 mm，實驗結果卻顯示最小壓降對應的氣速要略微大于臨界氣速。這表明最小壓降可能與臨界氣速并不完全對應。Skopich 等[11]也觀察到最小壓降有時存在于穩定流中，有時又存在于不穩定流中。

1.2 持液率突變模型

對于近水平氣液兩相流流動，僅改變氣相流量，由動量方程計算出的持液率可能會如圖1所示。圖1 中陰影區域對應的持液率有多個，學術上將該區域稱為持液率多解區。許多學者從持液率解穩定性的角度對多解區進行研究：Landman[12]利用KH 波理論對界面波穩定性進行分析，得出持液率的高解是不存在的，中間解可能存在，低解是始終存在的；Barnea 等[13]在Landman[12]的基礎上考慮了分層流的結構穩定性，指出中間解不滿足結構穩定性，只有最低解是穩定的。這表明在圖1多解區的左邊界處，持液率由A點突變到B點。

Espedal[14]、Langsholt 等[15]、Kjolaas 等[16]、Fan[17]均在實驗中觀察到了持液率突變現象。Biberg 等[18]將分層流持液率突變點作為臨界積液點。沈偉偉[19]綜合采用“最小滑脫”準則和分層流多解方法判斷持液率突變，從而計算相應的臨界氣速。該模型的缺點是適用范圍可能較窄，因為當管道傾角和液相負荷較大時，持液率可能不突變。

圖1 持液率多解區Fig.1 Regional schematic diagram for multiple liquid holdup solutions

1.3 零液壁剪切應力模型

零液壁剪切應力模型最早由Turner 等[20]提出，其認為環狀流液膜轉向的臨界狀態是液壁剪切應力為零，將液壁剪切應力為零的狀態作為氣井積液的臨界點；然而他們認為零液壁剪切應力模型的預測值與現場數據不符，采用液滴模型預測臨界氣速。

Fan[17]將該模型引入近水平管的分層流模型，將液膜最底部的壁面剪切應力為零作為積液的臨界判據。Alsaadi[21]在Fan[17]的基礎上進行簡化，將液膜的平均壁面剪切應力為零作為積液的臨界判據，進而把混合動量方程中的氣速和持液率解耦，節約了計算時間。若將分層流的液膜單獨作為分析對象，氣液界面剪切應力是主動力，液壁剪切應力是被動力，液壁剪切應力為零是液膜轉向的現象，其轉向的本質原因是氣體攜帶能力不足，從界面剪切應力的角度出發能夠進一步揭示積液的原因。

1.4 最小界面剪切應力模型

Barnea[22]提出了垂直管環狀流的最小界面剪切應力模型，其假設環狀流液膜沿管道周向均勻分布，得出對于一定的液相表觀速度，界面剪切應力僅為液膜厚度/持液率的函數，如圖2 所示。圖2 中界面剪切應力曲線存在最小值，Barnea[22]認為在最小界面剪切應力左側的環狀流是穩定的，其右側的環狀流是不穩定的。Luo 等[23]和Shekhar 等[24]均認為液膜沿管道周向均勻分布的假設并不合理，并做了相應修正：Luo等[23]用管道底部的最大液膜厚度進行計算；Shekhar 等[24]根據Paz 等[25]的實驗數據，考慮了最大液膜厚度與平均液膜厚度的關系，最終通過平均液膜厚度求出臨界氣速。沈偉偉等[26]在Shekhar等[24]的研究基礎上考慮了液滴夾帶的影響。

圖2 最小界面剪切應力示意圖Fig.2 Schematic diagram of minimum interfacial shear stress

Brito[27]將最小界面剪切應力模型從環狀流擴展到了分層流，并認同最小界面剪切應力左側對應穩定流，右側對應不穩定流的觀點。本文對穩定區與不穩定區的劃分也是如此。Brito[27]通過窮舉的算法尋找最小界面剪切應力對應的臨界氣速。 Brito 模型的不足是通過Blasius 經驗關系式來計算液壁剪切應力，因為Kowalski[28]測量氣液兩相分層流壁面剪切應力的實驗顯示：氣壁剪切應力與Blasius 關系式吻合得較好，但液壁剪切應力與Blasius 關系式明顯不符。

2 模型建立

Langsholt 等[15]氣液兩相流的實驗顯示：氣液界面波在高氣速下是3D小尺度波（滾波），隨著氣速減小，界面波的波幅逐漸變大，依次出現2D 短波和2D長波，最終流型由分層流變為假段塞流。Fan 等[8]將實驗中出現的分層流流型分為孤波（光滑分層流）和滾波兩種，孤波在高氣速下出現，滾波在接近流型轉變時出現。Alsaadi[21]的實驗中所有的分層流均是滾波分層流和長波分層流，長波分層流在接近流型轉變時出現。以上實驗現象均表明分層流在失穩過程中，界面波的波幅會變大。

可以推測，積液和界面波的變化有關，界面波的變化又直接與界面摩擦因子和界面剪切應力相聯系。可以認為積液是由于氣相給予液相的拖曳力不夠，即認為氣液界面剪切應力為最小值時為積液臨界狀態。

石油工業濕氣管道的管徑一般在2"以上，例如我國某頁巖氣開發采氣管線公稱直徑為DN65，濕氣集輸管道的最大管徑為DN550。Fan[29]在3"管道的低液相負荷氣液兩相流動的實驗中觀察到氣液界面是接近水平的。鄧道明等[30]曾建立高壓大管徑天然氣兩相流動計算模型，通過將模型計算結果與生產數據比較，認為平界面分層流模型與高壓大直徑天然氣管道工藝計算更為貼合。這里假設分層流的界面為水平界面，即管道橫截面液膜分布如圖3 所示。利用Biberg[31]基于雙極坐標系導出的液相速度分布解析式，將分層流界面剪切應力最小作為積液判據，建立新的濕氣管道積液預測模型。

圖3 濕氣管道內分層流液膜分布Fig.3 Film distribution of stratified flow in wet gas pipeline

2.1 動量方程

由于分層流流型是典型的分離流，故采用經典的雙流體模型；該模型的主要特點是考慮了相間的相互作用力。氣、液相的動量方程為：

式中，下角標g 和l分別代表氣相和液相，i表示界面；幾何參數A 代表各相所占管道的橫截面積，S代表各相的濕周；剪切應力τw表示各相壁面剪切應力，τi表示界面剪切應力；dp/dz表示各相沿流動方向（管道軸向）壓力梯度；ρgsinθ 表示各相重力沿軸向的分力。

2.2 分層流液膜區液相流量公式化

由于濕氣管道氣液兩相流動是低液相負荷的流動，其液相Reynolds 數一般要比轉捩Reynolds 數（Re=2000）低，因此其液膜流動可認為是層流。Biberg[31]曾基于雙極坐標系，假設氣液兩相分層流的氣相為湍流、液相為層流、氣液界面為平界面，并且假設氣液界面上的剪切應力均勻分布，利用N-S 方程推導出了雙極坐標系下液膜區的速度分布。對液膜區的液相速度求面積分，可以得到液相流量：

φB(φ)和φi(φ)是液面角φ 的函數。如圖4 所示，當φ∈[0,π/2]時，φB(φ)和φi(φ)均單調遞增，且當φ→0時，φB(φ)和φi(φ)→0。

圖4 權函數φB(φ)和φi(φ)Fig.4 Weight functions φB(φ)and φi(φ)

2.3 界面剪切應力

忽略液相斷面上由位置高度的變化而引起的壓力梯度，則氣、液兩相的壓力梯度dp/dz應相等，將式(3)代入氣相動量方程式(1)，可以得到界面剪切應力τi與液面角φ和氣壁剪切應力τwg之間的關系：

式中，氣壁剪切應力τwg的計算式為式(8)，其中氣壁摩擦因子fg的計算采用Blasius關系式[32]。

將式(8)代入式(7)，并利用Biberg[31]推出的液面角φ與持液率Hl的顯式表達式，得式(9)

圖5 τi與Hl和usg的關系[式(9)]Fig.5 Interfacial shear stress τivs holdup Hl and superficial gas velocity usg[Eq.(9)]

表1 3″上傾空氣-水管道積液實驗基本參數Table 1 Air-water experimental parameters in 3″upwardly inclined pipeline

由式(9)可知，當D、θ、ρg、ρl、usl已知時，τi僅是Hl和usg的函數；圖5 為在表1 參數下τi與Hl和usg關系的曲面圖，工程上濕氣集輸管道因其較低的液相負荷，氣液界面剪切應力與持液率和氣相折算速度關系也基本如圖5 所示。由圖5 可知，當usg不變、Hl∈[0,0.5]（即φ∈[0,π/2]）時，由于式(9)等號右邊第一項隨Hl的減小而遞增，第二項隨Hl的減小而遞減，使得圖5 呈現下凹的特性，在極值點?τi(Hl)/?Hl=0 處取得最小界面剪切應力τi_min，對應的持液率稱為臨界持液率Hl_c。另外，從圖5 中可知，Hl_c對usg的變化不敏感（圖中與usg軸基本平行的深色帶區域），這為后續快速求解臨界氣速創造了條件。

τi還可通過封閉關系式(10)得出，式中的界面剪切因子fi計算值對臨界氣速的預測至關重要；但是目前的實驗手段并不能直接測出τi的大小，fi的作用機理仍然是一個未完全解決的問題，目前似乎還沒有一種普適性的fi經驗關系式。本文模型中的fi采用Alsaadi[21]提出的關系式(11)。

圖6 與τi的交點Fig.6 Intersection of τi from Eq.(9)and from Eq.(10))。

2.4 臨界氣速預測

3 模型評估

圖7 臨界氣速計算流程Fig.7 Flow chart for calculating critical gas velocity

為了對新模型進行評估，本文收集了現有文獻中可以獲得的88 組實驗數據用于評估模型。數據來源于表2 中的文獻，其中Alsaadi[21]和Rodrigues[33]提供了分層流向假段塞流轉變的數據，其余文獻給出了臨界氣速的實驗數據。但是各研究者對于臨界氣速的識別方法有所不同。Langsholt 等[15]將持液率突變對應的氣速作為臨界氣速；Brito[27]、Fan[17]、Nair[34]、Alsaadi等[9]將液膜回流對應的氣速作為臨界氣速。上述幾位研究者對液膜回流的檢測方法又有所不同，Brito[27]和Nair[34]往液膜中注入鹽水，然后通過探針來感知液膜有無回流；Fan[17]則是往液膜注入顏料，通過注入點上游顏色變化來判斷液膜有無回流；Alsaadi 等[9]僅通過高速攝像機對液膜回流進行感知，由于攝像機只能捕捉到相對宏觀的現象，因此其探測液膜回流對應的氣速精度可能比前兩種方法低。

如圖8所示，對于接近水平管道來說，不論是實驗數據還是模型預測值都存在如下趨勢：隨著傾角和液相負荷的增加，臨界氣速也隨之增加。

就臨界氣速而言，新模型預測值普遍高于Alsaadi 等[9]的實驗值，但略微低于Fan[17]的實驗值；整體而言新模型預測值更接近Fan[17]的實驗結果。對此有以下解釋。

（1）Alsaadi 等[9]與Fan[17]實驗的管徑均為3″，壓力為常壓，流體介質均為空氣/水，當液相表觀流速為0.01 m/s 時，兩者的實驗條件一致，但是兩者的臨界氣速實驗結果有差異。這可能是由于檢測液膜回流的方法不同所致，前者僅采用高速攝像機，后者注入了顏料，前者感知力要弱于后者，故前者的臨界氣速實驗值要小于后者的實驗值。

表2 積液實驗條件參數匯總Table 2 Parameters of liquid accumulation experiments

圖8 新模型與文獻[9,17]的實驗數據對比Fig.8 Comparison of predictions from new model with experimental data Ref.[9,17]studies

（2）界面摩擦因子對于臨界氣速的預測至關重要，但目前沒有一種界面摩擦因子的經驗關系式能夠適用所有工況。所以模型的預測值與實驗值的偏差也可能是由于界面摩擦因子的不準確所致。

圖9展示了各模型預測的臨界氣速與來源于表2 中88 個相應實驗數據點之間的偏差。其中Alsaadi[21]和Rodrigues[33]的實驗臨界氣速取的是分層流最小氣速和假段塞流最大氣速的平均值。所選取的模型為各判據中最具有代表性的模型，其中新模型和Brito[27]模型均是以分層流的最小界面剪切應力為積液判據，Alsaadi[21]模型以分層流的零液壁剪切應力為積液判據，而沈偉偉[19]模型基于持液率突變準則。圖9的主圖表示各模型預測值與全部實驗數據的偏差。整體上看，圖9(b)Brito[27]模型的預測效果最差，Brito[27]模型的預測值較多位于±25%相對偏差線之外；另外3 組模型的預測值則大部分在±25%偏差線之內。圖9的縮略圖為各模型預測值與來源于Rodrigues[33]、Langsholt 等[15]、Espedal[14]的比較接近高壓力、大管徑實驗值的對比；其中Langsholt 等、Espedal 用較重的SF6 代替空氣，用于模擬高壓力下氣相。從縮略圖中可明顯看出只有新模型的預測值基本在±10%相對偏差線之內。

表3所示，與全部實驗數據點相比，新模型預測值的平均相對偏差為5%，標準差為16%；而與接近高壓力、大管徑的數據對比，新模型的平均相對偏差為-1%，標準差為8%。綜合來說，本文模型最優，文獻[19]模型次之，文獻[27]模型的預測結果與實驗數據偏差最大。

表3 各模型預測偏差Table 3 Prediction deviations for new and other models

4 結 論

積液現象發生的本質是氣相對液相的攜帶能力不足。基于最小界面剪切應力的積液判據，引入平氣液界面分層流液膜區的速度場描述，建立了近水平濕氣管道積液臨界氣速預測新模型，提出了求解臨界氣速的簡捷算法。本研究結論如下。

（1）對于近水平管道，模型預測的臨界氣速隨著傾角和液相負荷的增加而增加，這與實驗的變化趨勢一致。

圖9 各模型預測的臨界氣速與表2中實驗值的偏差Fig.9 Comparison of predicted critical gas velocities from new and other typical models with experimental data

（2）利用已有的實驗測量臨界氣速數據，對新模型和Brito、Alsaadi、沈偉偉三種近水平管積液模型進行評估，結果顯示新模型預測的臨界氣速與已有的實驗數據最為吻合，沈偉偉的模型次之，Brito模型的預測結果與實驗數據偏差最大。新模型與全部實驗數據的整體平均相對偏差為5%，標準差為16%。

（3）與其中的接近高壓力、大管徑數據相比，新模型預測精度更加優于其他模型，此時新模型平均相對偏差為-1%，標準差為8%。

符 號 說 明

A——面積，m2

d——管內徑，m

Fr——Froude數

f——摩擦因子

g——重力加速度，m/s2

Hl——持液率

hl——液膜厚度，m

p——壓力，Pa

dp/dz——壓降梯度，Pa/m

Ql——液相體積流量，m3/s

R——管道半徑，m

Re——Reynolds數

S——濕周，m

s——標準差

ug——氣相實際流速，m/s

ul——液相實際流速，m/s

usg——氣相表觀流速，m/s

usg_c——臨界氣相表觀流速，m/s

usl——液相表觀流速，m/s

ε——絕對偏差

θ——管道軸向與水平方向的夾角，rad

μ——動力黏度，Pa·s

ρ——密度，kg/m3

τi——界面剪切應力，Pa

τw——壁面剪切應力，Pa

φ——液面角，即氣液界面周長對應的半圓心角，rad

下角標

g——氣相

i——界面

l——液相

w——壁面