周朝 何祖清 付道明 羅璇 劉歡樂 孫志揚

1. 中國石化石油工程技術研究院；2. 頁巖油氣富集機理與有效開發(fā)國家重點實驗室；3. 中國石化江漢油田石油工程技術研究院

0 引言

多段壓裂頁巖氣水平井在生產(chǎn)早期有大量的壓裂液返排到井筒，造成頁巖氣井長期帶液生產(chǎn)，井筒積液問題不容忽視［1］。臨界攜液流量法是現(xiàn)場普遍采用的積液預測方法，然而，以Turner 模型、Coleman 模型為代表的球狀液滴模型對于國內(nèi)氣井的積液預測精度并不理想，考慮液滴變形的李閩模型對于垂直井筒的積液預測精度較高，但是并未考慮井斜變化和產(chǎn)液量對井筒積液的影響［2］。之后，王志彬等［3］、潘杰等［4］考慮液滴變形和產(chǎn)液量，建立了垂直井筒臨界攜液流量模型，Belfroid等［5］、李麗等［6］、Chen 等［7］考慮井斜角變化，建立了傾斜(水平) 井筒臨界攜液流量模型。Wang 等［8］根據(jù)室內(nèi)實驗數(shù)據(jù)修正Belfroid 模型，得到考慮井斜和液相流速的連續(xù)油管臨界攜液相關式，但應用條件受限。綜合分析，目前的臨界攜液流量模型均未完整反映頁巖氣井的復雜井身結(jié)構(gòu)和返排液量變化特征，為此，以液滴模型為基礎，綜合考慮井筒產(chǎn)液量、液滴變形和造斜率變化，根據(jù)變形液滴能量平衡關系，建立了頁巖氣井全井筒臨界攜液流量模型。

1 模型建立

1.1 物理模型

以橢球形變形液滴為基礎(圖1)，推導建立了臨界攜液流量新模型。液滴從球形變形為橢球形，以單個變形液滴為研究對象，主要進行如下假設：

圖 1 液滴變形示意圖Fig. 1 Schematic droplet deformation

(1)球形液滴以球心所在平面發(fā)生對稱變形；

(2)球形液滴變形后為光滑的標準橢球體，橢球體的兩個長軸長度相等；

(3)氣流方向垂直于橢球形液滴兩個長軸所在平面；

(4)忽略液滴與氣流之間的質(zhì)量交換和熱量交換。

1.2 臨界攜液流量模型

根據(jù)液滴動力學理論，當液滴在氣流中保持滯止狀態(tài)時，作用在液滴上的作用力相互平衡，在垂直井筒中，液滴所受重力等于液滴所受浮力與氣體曳力之和，即

代入各個力的表達式，有

式中，F(xiàn)G為液滴重力，N；Fb為液滴所受浮力，N；Fd為氣體曳力，N；ρl為液相密度，kg/m3；Vd為液滴體積，m3；g為重力加速度，m/s2；ρg為氣相密度，kg/m3；Cd為曳力系數(shù)，無因次；Ad為液滴迎風面積，m2；ug為氣相流速，m/s。

液滴從球形變形為橢球形后，式(2)可改寫為

式中，de為橢球形液滴長軸長度，m；h為橢球形液滴短軸長度，m。

根據(jù)液滴攜液理論，如果氣流能將井筒中的最大液滴攜帶到井口，則氣井能夠連續(xù)攜液，此時的氣體流速稱為臨界攜液流速［9］。定義橢球形液滴的軸比α為h/de，并考慮臨界攜液時的最大穩(wěn)定橢球形液滴長軸長度為dec，則有

式中，ucv為垂直井筒臨界攜液流速，m/s；α為橢球形液滴的軸比，小數(shù)；dec為最大穩(wěn)定橢球形液滴長軸長度，m。

根據(jù)橢球形液滴軸比實驗擬合得到的相關式為［10］

將式(5)代入式(4)整理，得到垂直井筒臨界攜液流速

在井筒傾斜段和水平段，由于造斜率和井身結(jié)構(gòu)的變化，會影響井筒流體的流動。在造斜率發(fā)生變化的地方，液滴與井筒碰撞并反彈，氣體及其夾帶的液滴的流動方向發(fā)生改變，以匹配新的造斜率和井身結(jié)構(gòu)?？紤]傾斜段和水平段中液滴與管壁碰撞和反彈造成的能量損失，對垂直井筒臨界攜液流速進行修正［11］，得到全井筒臨界攜液流速

式中，ucw為全井筒臨界攜液流速，m/s；R為最大造斜率(最大全角變化率)，(°)/30 m。

全井筒臨界攜液流量為

式中，qc為臨界攜液流量(標況)，m3/d；A為管道橫截面積，m2；p為壓力，MPa；Z為氣相偏差系數(shù)，無因次；T為溫度，K。

為求取臨界攜液流量，需首先確定關鍵參數(shù)dec、Cd和表面張力σ，并選取適合的井筒兩相流模型。根據(jù)式(8)計算出頁巖氣井井筒各個深度處的臨界攜液流量后，取各臨界攜液流量的最大值作為井筒積液判斷標準。

2 關鍵參數(shù)確定和模型求解

2.1 最大穩(wěn)定橢球形液滴尺寸

在環(huán)霧流中，液滴在氣相湍流作用下向上運動，或者下降形成積液。湍流動能表征了湍流的劇烈程度［12］，氣相湍流動能與液滴表面自由能間的平衡關系決定了最大穩(wěn)定橢球形液滴尺寸。最大穩(wěn)定橢球形液滴的總表面自由能最小，達到臨界積液狀態(tài)。

2.1.1 最大穩(wěn)定橢球形液滴總表面自由能

單個最大穩(wěn)定橢球形液滴的表面積為

式中，S為最大穩(wěn)定橢球形液滴的表面積，m2。

表面自由能等于產(chǎn)生表面面積所做的功［13］，則最大穩(wěn)定橢球形液滴的總表面自由能為

式中，ESmin為最大穩(wěn)定橢球形液滴總表面自由能，W；σ為氣液表面張力，N/m；vsl為液相表觀流速，m/s。

2.1.2 氣相湍流動能

井筒中單位體積氣相湍流動能為

式中，eT為單位體積氣相湍流動能，W/m3；v′r為氣相湍流平均徑向速度脈動，m/s；v′θ為氣相湍流平均切向速度脈動，m/s；v′z為氣相湍流平均軸向速度脈動，m/s。

由于平均徑向速度脈動的均方根近似等于摩擦速度，則有［14］

式中，vsg為氣相表觀流速，m/s；fsg為以氣相表觀流速流經(jīng)光滑管道的氣相摩阻系數(shù)，無因次。

在環(huán)霧流中通常為充分發(fā)展的氣相湍流流動，可假設氣相湍流為各向同性，有v′r=v′θ=v′z，將式(12)代入式(11)，則井筒中氣相總湍流動能為［15］

式中，ET為氣相總湍流動能，W。

2.1.3 最大穩(wěn)定橢球形液滴長軸長度

當井筒處于積液的臨界狀態(tài)時，最大穩(wěn)定橢球形液滴總表面自由能與氣相湍流動能相互平衡，并且成正比［15］

其中

式中，Ce為比例系數(shù)，無因次；θd為井斜角，°。

將式(10)和式(13)代入式(14)，可得最大穩(wěn)定橢球形液滴長軸長度

氣相摩阻系數(shù)可用Blasius 公式計算［16］，對于充分發(fā)展的氣相湍流流動，雷諾數(shù)Re＞2 100，則有

式中，μg為氣相黏度，Pa·s；D為管道內(nèi)徑，m。

對于垂直井筒，有Ce=0.75，并考慮到達到最大穩(wěn)定橢球形液滴時有vsg=ucv，將式(17)代入式(16)整理得

式中，Ql為液相流量，m3/s。

2.2 曳力系數(shù)

曳力系數(shù)的計算比較復雜，對于球形液滴的曳力系數(shù)，可采用Stokes 定律、Newton 公式和標準阻力曲線近似計算［17］。但是，對于橢球形變形液滴來說，其迎風面積與球形液滴相比發(fā)生變化，液滴顯著變形的情況下，變形液滴的曳力能達到相同條件下剛性球體曳力的兩倍以上［18］，因此，球形液滴的曳力系數(shù)計算方法對于橢球形液滴不再適用［19］。在井筒高雷諾數(shù)條件下，考慮液滴內(nèi)部流動，對剛性橢球體的曳力系數(shù)進行修正，得到橢球形液滴的曳力系數(shù)［3］

2.3 表面張力

目前的臨界攜液流量模型中，常取表面張力為常數(shù)，但是，表面張力的取值受到液相含量、溫度和壓力的影響，沿井深是差異分布的，如果取表面張力為常數(shù)，會造成臨界攜液流量計算結(jié)果出現(xiàn)偏差，最大偏差可達40%［20］。頁巖氣井中的液相為返排壓裂液，Sutton 公式［21］適用于天然氣-鹽水體系，能夠反映頁巖氣井產(chǎn)出頁巖氣和返排壓裂液的生產(chǎn)特點，因此選用Sutton 公式計算氣-液表面張力

式中，Cs為礦化度，mg/L。

2.4 臨界攜液流量模型求解

式(7)、式(8)構(gòu)成了頁巖氣井全井筒臨界攜液流量模型，模型中的關鍵參數(shù)dec、Cd和σ分別由式(18)、式(19)、式(20)確定，其中dec和Cd反映了井筒產(chǎn)液量與液滴變形的影響。由于臨界攜液流速ucw、橢球形液滴軸比α和曳力系數(shù)Cd均與最大穩(wěn)定橢球形液滴長軸長度dec相關，可消去ucw，采用Newton-Raphson 迭代法求解

其中

根據(jù)式(21)求解出dec后，即可求解出ucw和qc。

2.5 井筒兩相流模型優(yōu)選

臨界攜液流量受到井筒溫壓分布的影響，為了準確計算臨界攜液流量，需要優(yōu)選適合的兩相流模型。由于頁巖氣水平井井斜角變化較大，因此選取可以計算水平和傾斜井筒壓降的兩相流模型：Beggs-Brill 模 型［22］、Baker Jardine 模 型［23］、Dulker-Flanigan 模型［24-25］和Mukherjee-Brill 模型［26］。根據(jù)涪陵現(xiàn)場39 口頁巖氣井的實測油(套)壓數(shù)據(jù)和生產(chǎn)數(shù)據(jù)，分別利用4 種兩相流模型計算井底流壓，并與實測井底流壓結(jié)果比較。4 種兩相流模型的計算結(jié)果如圖2 所示，誤差分析如表1 所示?？梢钥闯觯琈ukherjee-Brill 模型的平均相對誤差的絕對值和均方根誤差均最小，壓力計算精度最高，可用于計算頁巖氣水平井的井筒壓力分布，滿足現(xiàn)場壓力計算精度要求。

表 1 兩相流模型誤差分析Table 1 Error analysis of two-phase flow models

圖 2 兩相流模型壓力計算結(jié)果Fig. 2 Pressure calculation result of two-phase flow model

3 應用實例

3.1 現(xiàn)場參數(shù)

根據(jù)涪陵現(xiàn)場39 口頁巖氣井的井筒溫壓測試數(shù)據(jù)(以最大井筒流壓梯度0.45 MPa/100 m 作為積液界限參考［27］)和生產(chǎn)動態(tài)數(shù)據(jù)，得到氣井實際積液狀態(tài)，并選用Belfroid 模型［5］、修正李閩模型［2］、李麗模型［6］、Chen 模型［7］、Wang 模型［8］和本文新模型分別計算現(xiàn)場氣井的臨界攜液流量，根據(jù)積液預測精度對新模型進行驗證。其中的修正李閩模型是指應用Fiedler 形狀函數(shù)修正［28］，使模型可應用于傾斜(水平)井筒的積液預測。

Belfroid 模型和修正李閩模型具有類似的表達式

式中，uc為臨界攜液流速，m/s；C為系數(shù)，Belfroid 模型中C取6.6，修正李閩模型中C取2.5；fw為壁面摩阻系數(shù)，無因次；NB為Bond 數(shù)，無因次。

39 口頁巖氣井生產(chǎn)數(shù)據(jù)和實際積液狀態(tài)如表2所示。天然氣相對密度0.56~0.57，返排液密度1.01~1.05 g/cm3，油管內(nèi)徑為62 mm 或50.6 mm，套管內(nèi)徑為115 mm。采用Yarborough-Hall 方法［29］計算氣相偏差系數(shù)。

表 2 頁巖氣井生產(chǎn)數(shù)據(jù)和積液狀態(tài)Table 2 Field parameters of shale gas wells

3.2 積液預測結(jié)果分析

根據(jù)表2 的頁巖氣井現(xiàn)場參數(shù)，以及各井的井筒溫壓分布數(shù)據(jù)，分別用6 種臨界攜液流量模型計算全井筒臨界攜液流量，并取所有臨界攜液流量中的最大值與實際產(chǎn)氣量進行比較，判斷井筒積液狀態(tài)。6 種模型的積液預測精度和預測結(jié)果如表3 和圖3 所示，圖3 中的對角線左上方區(qū)域表示氣井不積液，右下方區(qū)域表示氣井積液，而處于接近積液狀態(tài)的氣井，其數(shù)據(jù)點應位于對角線±15% 偏差范圍內(nèi)［4］(兩條虛線之間的區(qū)域)。將6 種臨界攜液流量模型的積液預測結(jié)果與井筒實際積液狀態(tài)比較，得到各模型的積液預測精度。

由表3 和圖3 可知，對于頁巖氣水平井的積液預測，Belfroid 模型的預測結(jié)果偏差最大，Belfroid 模型和Wang 模型均高估了臨界攜液流量，因此這兩種模型可以準確預測積液井，但是對不積液井的預測誤差較大；修正李閩模型低估了臨界攜液流量，因此可以準確預測不積液井，但是對積液井的預測誤差較大；李麗模型和Chen 模型對于各種積液狀態(tài)井的預測結(jié)果均有偏差；新模型的積液預測結(jié)果與實際積液情況的符合率最高，總體積液預測精度達92.3%，新模型可以準確預測積液井和接近積液井，只是對3 口不積液井的預測結(jié)果出現(xiàn)偏差，盡管如此，新模型對不積液井的預測精度仍然能夠滿足現(xiàn)場應用需求。綜合來看，新模型對于頁巖氣水平井的積液預測精度最高，可以有效指導現(xiàn)場積液判斷與排采工藝選擇。

4 結(jié)論

(1)依據(jù)液滴動力學分析和能量分析，綜合考慮井筒產(chǎn)液量、液滴變形和造斜率變化引起的液滴能量損失，根據(jù)變形液滴能量平衡關系，建立了頁巖氣井全井筒臨界攜液流量模型，彌補了現(xiàn)有模型沒有同時考慮頁巖氣水平井復雜井身結(jié)構(gòu)和返排液量變化的缺陷。

(2)根據(jù)變形液滴能量平衡關系，確定了最大穩(wěn)定變形液滴長軸長度；選取了適用于頁巖氣水平井的曳力系數(shù)公式和表面張力公式；優(yōu)選Mukherjee-Brill 兩相流模型用于計算頁巖氣水平井的井筒壓力分布，均方根誤差為8.96%，能滿足現(xiàn)場壓力計算精度要求。

(3) 39 口頁巖氣井的實例分析表明，與現(xiàn)有臨界攜液流量模型相比，新模型對于頁巖氣水平井的積液預測符合率最高，積液預測精度達92.3%，新模型可以有效指導頁巖氣水平井積液判斷。