王武杰，崔國民，魏耀奇，潘杰

（1.上海理工大學(xué)能源與動力工程學(xué)院，上海市動力工程多相流動與傳熱重點實驗室，上海 200093；2.國家石油天然氣管材工程技術(shù)研究中心，陜西寶雞 721008；3.西安石油大學(xué)石油工程學(xué)院，西安 710065）

0 引言

準確預(yù)測臨界攜液流速對診斷氣井積液狀況、解決含水氣井穩(wěn)產(chǎn)問題、挖掘煤炭地下氣化開發(fā)潛力至關(guān)重要[1-3]，同時，其對促進致密油氣快速規(guī)模發(fā)展，早日實現(xiàn)中國“能源獨立”戰(zhàn)略具有積極意義[4-5]。越來越多的實驗數(shù)據(jù)表明，傾斜氣井的攜液能力最差，預(yù)測也最為復(fù)雜和困難[6-7]。因此，有必要對傾斜氣井臨界攜液流速的預(yù)測展開深入研究。分層流作為傾斜管道氣液兩相流動的主要流型被普遍認可。分層流流型下，井筒內(nèi)氣液兩相的周向分布對流動特性有至關(guān)重要的影響，尤其對液相流量較小而氣相流量較大的情形更甚[8-9]。目前對臨界攜液狀態(tài)下的氣液相分布缺乏量化研究，且尚無傾斜管道氣液兩相臨界流動條件下的相界面摩擦因子計算模型?；诖耍疚姆謩e考察單位管長氣液兩相系統(tǒng)重力勢能、表面自由能和氣相動能的改變量，并根據(jù)能量最小原理確定臨界攜液狀態(tài)下氣液相分布的結(jié)構(gòu)參數(shù)；同時采用實驗數(shù)據(jù)修正相界面摩擦因子計算模型，使預(yù)測模型封閉；并將模型應(yīng)用于傾斜氣井臨界攜液流速的預(yù)測，以評估模型的應(yīng)用效果。

1 實驗系統(tǒng)

本文通過可視化實驗，模擬氣井?dāng)y液的動態(tài)過程，實驗裝置如圖1所示。

圖1 實驗裝置圖

分別采用水和甘油、水的混合溶液（二者體積分數(shù)均為 50%）模擬井下積液，利用空氣模擬天然氣。實驗條件下空氣密度ρa=1.184 kg/m3，空氣動力黏度μa=1.849×10-5Pa·s；水密度ρw=997.05 kg/m3，水動力黏度μw=8.900 8×10-4Pa·s；水、甘油混合溶液密度ρgw=1 139.56 kg/m3，水、甘油混合溶液動力黏度μgw=6.882×10-3Pa·s；空氣、水表面張力σaw=0.071 97 N/m，空氣與水、甘油混合溶液表面張力σagw=0.067 11 N/m。

分別選用外徑50 mm、壁厚5 mm和外徑70 mm、壁厚5 mm的有機玻璃管進行實驗，步驟為：①在管內(nèi)低洼處注入一定體積的液體，通過風(fēng)機增壓向玻璃管內(nèi)注入壓縮空氣，計量壓力和流量；②采用變頻器實現(xiàn)氣體流速的連續(xù)可調(diào)，直至低洼處的積液完全鋪展于上傾段并與管壁相對靜止（此時達臨界氣體流速）。

2 模型建立與求解方法

2.1 受力分析與控制方程

臨界氣體流速條件下，傾斜管內(nèi)流型為波狀分層流，液膜在重力、相界面和壁面切應(yīng)力的作用下處于平衡狀態(tài)。此時，近氣相側(cè)的液膜在相界面切應(yīng)力的作用下向下游運動，而近壁面?zhèn)鹊囊耗ぴ谥亓ψ饔孟禄亓?，即液膜?nèi)部存在循環(huán)流動。傾斜井筒微元段內(nèi)液膜受力如圖2所示。上述流動過程可簡化為穩(wěn)態(tài)一維流動，并忽略加速壓降，可對氣、液相分別建立動量方程：

圖2 傾斜井筒內(nèi)液膜流動與受力分析

流動過程中近似認為氣相與液相的平均壓力梯度相等，聯(lián)立（1）式和（2）式可得控制方程：

上式中，幾何參數(shù)與氣液相分布有關(guān)，應(yīng)力參數(shù)與持液率、各相速度和摩擦因子有關(guān)。通過閉合幾何參數(shù)和應(yīng)力參數(shù)可使控制方程閉合。

2.2 氣液相分布

根據(jù)Brauner等[10]的研究，液膜的周向分布如圖3所示。其中，O和O1分別表示井筒流道的截面軸心和虛擬圓圓心?？刂品匠讨械膸缀螀?shù)采用下式計算：

圖3 井筒周向的液膜分布

不同相界面假設(shè)下的幾何參數(shù)的計算式見表1。

表1 不同相界面下的幾何參數(shù)計算關(guān)系式

2.3 能量最小原理

在相界面形狀由水平變?yōu)閺澢倪^程中，重力、切應(yīng)力和表面張力將引起系統(tǒng)勢能、動能和表面自由能的變化。臨界條件下，單位管長氣液兩相系統(tǒng)總能的改變量可表示為：

基于能量最小原理，單位管長系統(tǒng)總能量最小時，相界面結(jié)構(gòu)最穩(wěn)定，即：

①單位管長重力勢能的改變量[10]

②單位管長表面自由能的改變量

③單位管長氣相動能改變量

根據(jù)七分之一次方定律，圓管內(nèi)湍流氣相的速度分布表示如下：

單位管長氣相動能改變量采用下式計算：

2.4 剪切應(yīng)力計算

Taitel等[11]給出了相關(guān)應(yīng)力計算表達式：

Banafi等[12]實驗研究指出 Taitel模型可準確預(yù)測氣液兩相低持液率流動條件下的氣壁面摩擦因子，表達式如下：

層流時，Cg=16，m=1；湍流時，Cg=0.046，m=0.2。

Biberg[13]指出，傾斜管道光滑分層流條件下，平均液壁面切應(yīng)力與平均相界面切應(yīng)力滿足下式：

上式中f（φ0P）是 難以獲得解 析 解 的積分項。Biberg[13]給出以下近似計算關(guān)系式：

式中，an和bn的取值見表 2[13]。該近似值的最大絕對誤差不超過2×10-6。

表2 近似關(guān)系式（（26）式）中參數(shù)an和bn的取值

準確計算相界面摩擦因子是預(yù)測相界面切應(yīng)力的關(guān)鍵。傾斜管道中，界面波和持液率對相界面摩擦因子的影響較大。界面波動特性主要取決于機械能從氣相向液相的傳遞速率以及液體中黏性耗散速率，而這一過程受傾角的影響[14]。隨著管道傾角的增加，重力對界面波生長的約束減弱，使得傾斜管中的液層表面更易于觀察到界面波。本文認為在傾斜管道中，重力對傾斜管內(nèi)臨界狀態(tài)下的界面波動作用更加顯著。因此，一方面將界面波等效為管壁的相對粗糙度，另一方面，引入傾角對相界面摩擦因子進行修正[14]，具體過程如下：

其中，無因次相界面濕周定義為：

再考慮傾角的影響對其進行修正：

式中，A和B與流體物性和管道條件有關(guān)[15]。

2.5 模型求解方法

控制方程（（3）式）可以視為持液率的函數(shù)，即：f（ε）=0。零始終是方程f（ε）=0的一個解。的函數(shù)圖像如圖4所示。當(dāng)usg=usg,cr時，f（ε）=0只有一個非零解，且該非零解在（0，0.1）區(qū)間求得，求解步驟如圖5所示。

圖 4 函數(shù) f（ε） 1/3在不同流動條件下的變化曲線

圖5 模型求解過程示意圖

3 結(jié)果分析與實例驗證

3.1 氣液相分布預(yù)測結(jié)果與分析

對控制方程無因次化處理后可得：

其中

由上式可知，臨界狀態(tài)下的氣液相分布與邦德數(shù)、氣相表觀韋伯?dāng)?shù)、持液率和井斜角有關(guān)?；谀芰孔钚≡恚炕瘹庖簝上嘞到y(tǒng)總能量隨相界面結(jié)構(gòu)的變化，并以空氣和水為例，分別考察上述 4個因素對管內(nèi)氣液相分布的影響，計算結(jié)果見圖6（圖中紅色數(shù)據(jù)點為各工況下的能量最小值點）?？梢钥吹剑孩匐S著管徑的增加，最穩(wěn)定的氣液相分布（對應(yīng)能量最小值）所對應(yīng)的相界面彎曲角從125°（微凸）逐漸趨于180°（水平）；反之，管徑越小，在潤濕性和表面張力的影響下相界面彎曲角越小，呈現(xiàn)凸型結(jié)構(gòu)的液相分布（見圖6a）；②隨著井斜角的增加，重力抵消潤濕性和表面張力的作用效果更加顯著，使最穩(wěn)定的氣液相分布所對應(yīng)的相界面彎曲角從120°（微凸）趨近于180°（水平）（見圖 6b）；③結(jié)合圖 6c、圖 6d可知，氣體流速和氣相密度增加，使氣相動能也增加，在潤濕性和表面張力的影響下相界面彎曲角減?。ㄗ兺梗┯欣诮档拖到y(tǒng)表面自由能；④隨著持液率的增加，液相重力勢能的影響增大，有利于穩(wěn)定原來水平相界面的形狀，同時界面微凸也可減小持液率增加導(dǎo)致的系統(tǒng)表面自由能的增加。但是在液相密度遠大于氣相密度的低壓系統(tǒng)下，液相重力勢能起主導(dǎo)作用（見圖6e），隨著持液率的增加，最穩(wěn)定的氣液相分布所對應(yīng)的相界面彎曲角從105°逐漸趨于165°（近水平）。綜上所述，低持液率流動條件下，受表面張力和潤濕性的影響，相界面呈現(xiàn)微凸，符合如圖7所示的實驗觀察結(jié)果。

圖6 液相分布的模型預(yù)測結(jié)果

圖7 實驗觀察的臨界狀態(tài)下的氣液相分布示意圖

3.2 實驗結(jié)果與分析

圖 8為臨界攜液氣體流速和臨界壓力梯度的實驗結(jié)果。由圖可知，臨界攜液流速和臨界壓力梯度隨著傾角的增加呈非線性增加。在相同的管道傾角下，液相密度或者黏度越大，臨界攜液流速越大。這與Birvalski等[16]的實驗結(jié)果一致。同時，管徑增加使臨界攜液氣體流速增加。這是因為在相同持液率下，界面形狀的變化不及Bo和Wesg數(shù)對臨界攜液流速的影響顯著，要使控制方程依然成立，必須提高氣體流速，以抵消管徑增加引起的Bo/Wesg的增量。

圖8 臨界攜液流速和臨界壓力梯度的實驗結(jié)果

3.3 傾斜管內(nèi)相界面摩擦因子計算模型

采用Birvalski等[16]和本文的實驗數(shù)據(jù)，基于簡化后的水平相界面形狀封閉幾何參數(shù)，對不同傾角下的相界面摩擦因子修正參數(shù)進行反算，并通過數(shù)據(jù)擬合得到相界面摩擦因子修正參數(shù)（f（α）），計算結(jié)果如表3所示。采用表 3中相界面摩擦因子修正參數(shù)對相界面摩擦因子計算公式（（31）式）進行修正后，模型的預(yù)測值達到了非常高的精度。皮爾遜相關(guān)系數(shù)相對較高。

表3 不同實驗條件下的相界面摩擦因子修正參數(shù)

由表3可知，f（α）的斜率為正值，即表明相界面摩擦因子隨著傾角的增大整體呈現(xiàn)遞增的趨勢。主要是因為傾角的增加使得約束界面波動的重力分量減小，導(dǎo)致界面波動更加劇烈和明顯。同時，f（α）小于1，表明在微傾管道內(nèi)低持液率流動條件下，相界面摩擦因子要小于相同氣體流速下水平管道相界面摩擦因子預(yù)測值，主要原因是當(dāng)管道接近水平時，傾角對界面波的影響不如持液率的影響顯著。本文基于實驗數(shù)據(jù)，結(jié)合無因次化分析，給出相界面摩擦因子修正參數(shù)的計算式：

3.4 模型對比與工程實例驗證

將模型預(yù)測結(jié)果與實驗數(shù)據(jù)進行對比（見圖9）。本文模型預(yù)測的臨界攜液流速隨管徑的增加而增加，隨液相密度和黏度的增加而增加（同管徑下，甘油溶液的臨界攜液流速大于水的臨界攜液流速），這一規(guī)律與符合 Birvalski等[16]和 Rastogi等[17]的實驗結(jié)果；臨界壓力梯度的預(yù)測結(jié)果隨管徑的增加而減小，符合實驗結(jié)果，管徑的增加使液相水力直徑增加是壓力梯度減小的重要原因。對比本文實驗數(shù)據(jù)，新模型預(yù)測臨界攜液流速和臨界壓力梯度的平均相對誤差分別為1.19%和3.02%。同時，臨界攜液流速隨傾角的增加先增大后減小，最大臨界攜液流速對應(yīng)的管道傾角隨管徑的增加而增加，同時也受流動介質(zhì)物性的影響。本文模型模擬計算發(fā)現(xiàn)，標(biāo)況下天然氣攜礦化水時，最大臨界攜液流速對應(yīng)的管道傾角分別為 30°（D=50.67 mm）、37°（D=62.00 mm）、40°（D=75.9 mm）和 45°（D=100.53 mm）。

圖9 本文模型預(yù)測結(jié)果與實驗結(jié)果對比

本文模型應(yīng)用于高壓系統(tǒng)時需要對預(yù)測結(jié)果進行轉(zhuǎn)換。本文從能量角度出發(fā)，認為高、低壓管道系統(tǒng)內(nèi)，臨界條件下單位體積氣液兩相間傳遞的機械能應(yīng)相等，即：

為便于比較，將臨界攜液流速轉(zhuǎn)換為標(biāo)況下的臨界攜液流量：

當(dāng)天然氣的實際流速大于傾斜氣井中臨界攜液流速最大值時，氣井才能連續(xù)攜液。氣田現(xiàn)場需結(jié)合方位測井資料來確定傾斜氣井的最大臨界攜液流速。選取中國普光氣田和延長氣田42口氣井的生產(chǎn)數(shù)據(jù)，并與實際產(chǎn)氣量繪制交會圖（見圖10），可以看到本文模型僅誤判1口氣井，誤判率2.38%，說明本文模型的預(yù)測結(jié)果具有較高的可信度，可對傾斜氣井積液進行有效判斷。

圖10 本文模型預(yù)測結(jié)果與對比

4 結(jié)論

低持液率流動條件下，相界面形狀更容易受到壁面潤濕性和表面張力的影響呈現(xiàn)彎曲，管徑越小、持液率越小、井斜角越小、氣體流速越大、氣相密度越大時，相界面彎曲越明顯。

氣井臨界攜液流速隨井斜角的增加呈現(xiàn)非線性先增后減的趨勢；最大臨界攜液流速對應(yīng)的井斜角隨井筒直徑的增加而增大，同時也受氣液兩相物性的影響。

本文模型預(yù)測臨界攜液流速與臨界壓力梯度的平均相對誤差分別為1.19%和3.02%，現(xiàn)場傾斜氣井積液誤判率2.38%，可對傾斜氣井積液進行有效判斷。

符號注釋：

an，bn，A，B——常系數(shù)，無因次；AT——油管截面面積，m2；Ag——氣相流動通道面積，m2；ΔAgl——單位管長氣液相接觸面面積的變化，m2/m；Al——液相流動通道面積，m2；ΔAlw——單位管長液壁接觸面面積的變化，m2/m；Bo——邦德數(shù)，無因次；Cg，m——與流態(tài)有關(guān)的系數(shù)，無因次；D——油管內(nèi)徑，m；DH,g——氣相水力直徑，m；Dl（φ0P）——不同液相分布角下的等效層流直徑，m；fg——氣相摩擦因子，無因次；fi——相界面摩擦因子，無因次；fl——液相摩擦因子，無因次；f（α）——相界面摩擦因子修正參數(shù)，無因次；f（φ0P）——與液相分布角有關(guān)的函數(shù)，無因次；F——中間變量，無因次；g——重力加速度，m/s2；k——相界面絕對粗糙度，m；n——常系數(shù)編號；p——壓力，Pa；Qcr——臨界攜液流量，m3/d；r——距離井筒軸心的徑向位置，m；R——油管半徑，m；Reg——氣相雷諾數(shù)，無因次；R1——虛擬圓半徑，m；Sg——氣相濕周，m；Si——相界面濕周，m；si——無因次相界面濕周，無因次；Sl——液相濕周，m；T——熱力學(xué)溫度，K；ug——氣相真實流速，m/s；g——平均氣相真實速度，m/s；ug,max——最大氣相真實速度，m/s；ui——相界面真實流速，m/s；ul——液相真實流速，m/s；usg——氣相表觀流速，m/s；Δusg——氣相表觀流速迭代步長，m/s；usg,cr——臨界攜液氣相表觀流速，m/s；usg,cr,hp——高壓系統(tǒng)下的臨界攜液氣相表觀流速，m/s；Wesg——氣相表觀韋伯?dāng)?shù)；x——流動方向坐標(biāo)，m；YGPC——彎曲相界面的重心位置，m；YGP——水平相界面的重心位置，m；Z——壓縮因子；α——管道傾角，（°）；β——井斜角，（°）；ΔEk,g——單位管長氣相的動量改變量，J/m；ΔEp——單位管長氣液兩相系統(tǒng)的重力勢能改變量，J/m；ΔEs——單位管長氣液兩相系統(tǒng)的表面自由能改變量，J/m；ΔEtot——單位管長氣液兩相系統(tǒng)的總能量改變量，J/m；δl——液膜厚度，m；δl（φ）——液膜厚度分布，m；ε——持液率；ε0——持液率迭代初值；Δε——持液率迭代步長；θ——潤濕角，rad；μa——空氣動力黏度，Pa·s；μg——氣相動力黏度，Pa·s；μgw——水、甘油混合溶液動力黏度，Pa·s；μl——液相動力黏度，Pa·s；μw——水動力黏度，Pa·s；ρa——空氣密度，kg/m3；ρg——氣相密度，kg/m3；ρg,hp——高壓系統(tǒng)下的氣相密度，kg/m3；ρgw——水、甘油混合溶液密度，kg/m3；ρl——液相密度，kg/m3；ρw——水密度，kg/m3；σaw——空氣、水界面的表面張力，N/m；σagw——空氣與水、甘油混合溶液的表面張力，N/m；σgl——氣液界面的表面張力，N/m；τi——相界面剪切應(yīng)力，Pa；——平均相界面切應(yīng)力，Pa；τig——相界面反作用于氣相的相界面切應(yīng)力，Pa；τgi——氣相作用于相界面的相界面切應(yīng)力，Pa；τwg——氣壁面剪切應(yīng)力，Pa；τwl——液壁面剪切應(yīng)力，Pa；τwl——平均液壁面切應(yīng)力，Pa；φ——任意液相分布角，rad；φ0——彎曲相界面的液相分布角，rad；φ0P——水平相界面的液相分布角，rad；φPA——相界面彎曲角，rad。下標(biāo)：g—氣相；l—液相。