李元生 ，李相方，藤賽男，徐大融，和向楠

1.中國石油大學(xué)石油工程教育部重點(diǎn)實(shí)驗(yàn)室，北京 昌平102249；2.中海石油（中國）有限公司上海分公司研究院，上海 徐匯200030 3.中國石化上海海洋油氣分公司研究院，上海 浦東200120；4.陜西延長石油（集團(tuán)）有限責(zé)任公司研究院，陜西 西安710000

引 言

常規(guī)的二項(xiàng)式產(chǎn)能方程中，將氣井的泄氣半徑當(dāng)成氣井高速非達(dá)西外邊界[1-2]，實(shí)際上高速非達(dá)西外邊界是隨著產(chǎn)量變化的，并且只存在于近井地帶比較小的范圍內(nèi)，從而會對氣井的產(chǎn)能預(yù)測造成很大影響。在研究高速非達(dá)西時(shí)一般用雷諾數(shù)表征，比較合理的公式為前蘇聯(lián)學(xué)者卡佳霍夫提出的表達(dá)式[3-4]，但是目前在應(yīng)用該公式到氣井產(chǎn)能方程中時(shí)，一方面將地層平均壓力[5-8]看作等于非達(dá)西動邊界處的壓力，而實(shí)際上非達(dá)西動邊界處的壓力與儲層的物性以及氣井工作制度有關(guān)不能進(jìn)行簡單的處理；另一方面方程的推導(dǎo)也并不完整，可以對氣體的密度進(jìn)一步簡化[3-7]。

1 高速非達(dá)西流動邊界半徑rh 的確定

對于二維徑向流，按滲流規(guī)律可將滲流區(qū)域分為兩個(gè)區(qū)域[8]：（1）達(dá)西流動區(qū)域，邊界為（rh，re）。該區(qū)域由于氣體滲流速度小，服從達(dá)西流動規(guī)律。（2）高速非達(dá)西流動區(qū)域，邊界為（rw，rh）。該區(qū)域由于氣體過流面積減小，加上氣體膨脹，氣體不服從達(dá)西滲流規(guī)律。 rh稱為非達(dá)西動邊界。目前的產(chǎn)能方程中通常認(rèn)為rh=re即整個(gè)滲流區(qū)域都是高速非達(dá)西流動區(qū)域。而實(shí)際整個(gè)區(qū)域可以用如圖1 與式（1）所示。

圖1 地層滲流示意圖Fig.1 Schematic of flow

由于在地層中的滲流和管流的相似性，可以用雷諾數(shù)值來判斷流體在地層中是否服從達(dá)西定律[4]。雷諾數(shù)有許多不同求法，目前公認(rèn)比較合理的是前蘇聯(lián)學(xué)者卡佳霍夫提出的表達(dá)式[3-4]

通過實(shí)驗(yàn)獲得臨界雷諾數(shù)為0.2～0.3[3-4]，雷諾數(shù)小于臨界雷諾數(shù)時(shí)流體流動服從達(dá)西定律，雷諾數(shù)大于臨界雷諾數(shù)時(shí)流體流動不服從達(dá)西定律。

假設(shè)非達(dá)西流動邊界rh處氣體的密度、體積系數(shù)和氣體的速度分別為

將臨界雷諾數(shù)以及式（3）～式（5）代入式（2）得到不同流量下的非達(dá)西動邊界半徑，如式（6）所示

從式（6）中可以看出非達(dá)西外邊界與氣井的產(chǎn)量有關(guān)，由于氣井產(chǎn)能方程為

結(jié)合式（6）與式（7）可以建立非達(dá)西動邊界與氣井產(chǎn)能方程的關(guān)系，進(jìn)行迭代求解可以得到高速非達(dá)西外邊界的大小，如圖2 所示。從圖中可以看出，隨著儲層滲透率的增大，高速非達(dá)西區(qū)域逐漸增大；隨著井底流壓的增大，生產(chǎn)壓差逐漸減小，高速非達(dá)西區(qū)域減小。并且滲透率在0.8～5.0 mD 變化時(shí)，非達(dá)西動邊界為0.07～1.20，根據(jù)式（7）可以發(fā)現(xiàn)高速非達(dá)西動邊界對方程右邊第二項(xiàng)系數(shù)的影響很大。

2 高速非達(dá)西臨界滲透率與臨界產(chǎn)量

2.1 非達(dá)西臨界滲透率

由式（6）可知，當(dāng)?shù)貙訁?shù)滲透率、孔隙度、氣體黏度和地層厚度都不變時(shí)，非達(dá)西動邊界只與氣體的產(chǎn)量有關(guān)。當(dāng)氣體的產(chǎn)量為絕對無阻流量時(shí)即qg=qAOF，最大的非達(dá)西動邊界為

當(dāng)井底流壓為大氣壓時(shí)，氣井的產(chǎn)能方程為

由式（8）、式（9）可知，只有當(dāng)最大的非達(dá)西外半徑要大于井筒的井眼即：rmaxh＞rw時(shí)，地層中才會有高速達(dá)西出現(xiàn)。故存在臨界滲透率使得當(dāng)氣井的產(chǎn)能為其無阻流量時(shí)，最大高速非達(dá)西動邊界等于井眼半徑即：rmaxh=rw。

結(jié)合式（8）、（9）可得臨界滲透率為

當(dāng)?shù)貙訚B透率大于臨界滲透率時(shí)，地層中才有可能出現(xiàn)高速非達(dá)西滲流。當(dāng)?shù)貙訚B透率小于臨界滲透率時(shí)，無論氣井以多少產(chǎn)量生產(chǎn)，高速非達(dá)滲流都不會在地層中產(chǎn)生。如圖3 所示，隨著地層壓力的增大，臨界滲透率逐漸減小。

圖3 臨界滲透率與地層壓力的關(guān)系曲線Fig.3 Curves between formation pressure and critical permeability

2.2 高速非達(dá)西臨界產(chǎn)量

當(dāng)滲透率大于KD時(shí)，氣井工作制度的變化影響了地層中的滲流形態(tài)。當(dāng)氣井產(chǎn)量很小時(shí)，高速非達(dá)西動邊界仍可能比rw小，地層中的滲流符合達(dá)西定律；當(dāng)氣井產(chǎn)量很大時(shí)，高速非達(dá)西動邊界比rw大，地層中線性流破壞。故存在使得非達(dá)西動邊界rh等于rw時(shí)的產(chǎn)量，定義該氣井產(chǎn)量為氣井高速非達(dá)西臨界產(chǎn)量

3 近井地帶高速非達(dá)西產(chǎn)能方程的建立

由于二項(xiàng)式產(chǎn)能方程中非達(dá)西項(xiàng)為[8]

非達(dá)西系數(shù)為

將式（6）與式（13）代入式（12），得

將式（15）代入氣井產(chǎn)能方程中，則產(chǎn)能方程為

其中

得到無阻流量

從式（15）和AD的表達(dá)式可以看出，非達(dá)西邊界的影響可以看成為高速非達(dá)西表皮

其中β 是與滲透率、孔隙度相關(guān)的函數(shù)，其綜合表達(dá)式可以寫成

許多學(xué)者通過對單相流體的實(shí)驗(yàn)研究得到不同的β 函數(shù)關(guān)系式，對于有裂縫的鹽水或者油氣體系，系數(shù)c 和d 可以為0，β 僅是與滲透率有關(guān)的函數(shù)，表1 為不同的學(xué)者通過回歸得到的不同表達(dá)式。

表1 僅與滲透率有關(guān)的β 表達(dá)式數(shù)據(jù)表Tab.1 Data of β expressions concerned only with permeability

不同的學(xué)者根據(jù)實(shí)驗(yàn)假設(shè)條件下也得到β 與孔隙度、迂曲度等參數(shù)的關(guān)系，如表2 所示。

表2 與滲透率、孔隙度以及迂曲度有關(guān)的β 表達(dá)式數(shù)據(jù)表Tab.2 Data of β expressions concerned with permeability，porosity and capillary tortuosity

從表1 和表2 可以發(fā)現(xiàn)，眾多學(xué)者回歸的公式中滲透率都處于分母位置，系數(shù)基本大于0.5；孔隙度既有處于分子也有處于分母位置，當(dāng)其處于分母時(shí)系數(shù)基本小于1.5；迂曲度處于分子位置。故將系數(shù)β 代入式（17）中可以得到表皮的一般表達(dá)式

其中，f =-5.48×10-9×bg=1.5-ce=a-0.5。

另外，由式（15）可知，高速非達(dá)西產(chǎn)生的條件為

將不同的β 表達(dá)式代入式（17）中可以得到不同的e，d，g 系數(shù)，現(xiàn)將常用的式（14）中β 表達(dá)式（1）代入式（17），計(jì)算并取臨界雷諾數(shù)NRel=0.3。則非達(dá)西的邊界表皮為

SD是地層滲透率和孔隙度的函數(shù)。其主要反映了非達(dá)西區(qū)域半徑的影響：孔隙度越大，高速非達(dá)西半徑也越小，高速非達(dá)西表皮越小；滲透率越大，高速非達(dá)西半徑也越大，高速非達(dá)西表皮越小。其中SD是該產(chǎn)能方程區(qū)別于常規(guī)二項(xiàng)式產(chǎn)能方程的地方。

4 氣井產(chǎn)能比較

結(jié)合式（16）和式（22）繪出常規(guī)二項(xiàng)式產(chǎn)能方程、達(dá)西線性流動產(chǎn)能方程以及本文推到產(chǎn)能方程在不同滲透率下的IPR 曲線，如圖4 所示。

從圖4 中可以看出，當(dāng)滲透率小時(shí)，氣井臨界產(chǎn)量大于氣井的無阻流量時(shí)，高速非達(dá)西滲流在地層不存在，地層中的滲流符合達(dá)西定律，故IPR 曲線完全與線性流動時(shí)的IPR 曲線重合如圖4a 前半段所示。當(dāng)氣井無阻流量大于臨界產(chǎn)量如圖4a～圖4c 所示，隨著滲透率的增大，本文產(chǎn)能方程的IPR曲線不再與線性流動產(chǎn)能方程的IPR 曲線重合，越來越向常規(guī)二項(xiàng)式產(chǎn)能方程靠近。當(dāng)滲透率足夠大時(shí)，本文推導(dǎo)的產(chǎn)能方程的IPR 曲線與常規(guī)二項(xiàng)式產(chǎn)能方程IPR 曲線在氣井產(chǎn)量大于臨界產(chǎn)量時(shí)完全重合。所以本文的產(chǎn)能方程是符合實(shí)際情況的。

圖4 不同滲透率下的IPR 曲線Fig.4 IPR curves at different permeability

5 結(jié) 論

（1）低滲氣藏中氣井高速非達(dá)西動邊界半徑要遠(yuǎn)小于氣井的泄氣半徑，只存在于井筒附近。并且滲透率越大高速非達(dá)西動邊界越大；產(chǎn)量越大非達(dá)西動邊界半徑越大。并且當(dāng)滲透率很低時(shí)，高速非達(dá)西動邊界對系數(shù)B 的影響很大，隨著動邊界的增大系數(shù)B 增大，但動邊界大于10 倍的井眼半徑時(shí)，系數(shù)B 與常規(guī)二項(xiàng)式產(chǎn)能方程的系數(shù)B 基本相等。

（2）當(dāng)氣井產(chǎn)量為其無阻流量時(shí)，存在使得高速非達(dá)西動邊界等于井眼半徑時(shí)的臨界滲透率，當(dāng)儲層滲透率大于該滲透率時(shí)，儲層的高速非達(dá)西才可能存在。并且臨界滲透率隨著地層壓力的變化較大，隨著地層壓力的增大，臨界滲透率變小。

（3）近井地帶高速非達(dá)西的產(chǎn)能方程與常規(guī)產(chǎn)能方程相比多了近井地帶高速非達(dá)西表皮，該表皮只與地層孔隙度、滲透率有關(guān)的參數(shù)。滲透率越大表皮越小，孔隙度越大表皮也越大。

（4）與常規(guī)的二項(xiàng)式產(chǎn)能方程相比，由于本文推導(dǎo)的產(chǎn)能方程是考慮近井地帶高速非達(dá)西動邊界的影響，更符合實(shí)際的生產(chǎn)情況，用二項(xiàng)式產(chǎn)能方程預(yù)測低滲氣藏的產(chǎn)能誤差很大，所以本文推導(dǎo)的產(chǎn)能方程對預(yù)測氣藏的產(chǎn)能方程更為準(zhǔn)確。

符號說明

p—壓力，MPa；

r—波及半徑，m；

μ—?dú)怏w黏度，mPa·s；

K—滲透率，mD；

υ—滲流速度，m/s；

β—速度系數(shù)，m-1；

ρ—?dú)怏w的密度，kg/m3；

re—泄氣半徑，m；

rh—高速非達(dá)西動邊界半徑，m；

rw—井眼半徑，m；

Re—雷諾數(shù)，無因次；

?—孔隙度，無因次；

ρh—非達(dá)西動邊界處的密度，kg/m3；

γg—?dú)怏w的相對密度，無因次；

ph—非達(dá)西動邊界處的壓力，MPa；

psc—標(biāo)況下的壓力，MPa；

T—地層溫度，K；

Zh—非達(dá)西動邊界處壓力對應(yīng)的偏差因子，無因次；

Bgh—非達(dá)西動邊界處壓力對應(yīng)的體積系數(shù)，無因次；

Tsc—標(biāo)況下的溫度，K；

υh—非達(dá)西動邊界處的速度，m/s；

qg—產(chǎn)量，m3/d；

h—儲層厚度，m；

Rel—臨界雷諾數(shù)，無因次；

ˉpR—地層平均壓力，MPa；

pwf—井底流壓，MPa；

rmaxh—最大高速非達(dá)西動邊界，m；

qAOF—無阻流量，m3/d；

KD—高速非達(dá)西臨界滲透率，mD；

Z—偏差因子，無因次；

qminh—高速非達(dá)西臨界產(chǎn)量，m3/d；

Δp—壓差，MPa；

D—紊流系數(shù)，（m3/d）2；

SD—高速非達(dá)西邊界表皮，無因次；

a，b，c，d—速度系數(shù)表達(dá)式的系數(shù)，無因次；

τ—迂曲度，無因次。

[1] 李士倫.天然氣工程[M].北京：石油工業(yè)出版社，2008.

[2] 羅伯特，沃特恩伯格.氣藏工程[M].王玉普，譯.北京：石油工業(yè)出版社，2007.

[3] 楊勝來，魏俊之.油藏物理學(xué)[M].北京：石油工業(yè)出版社，2004.

[4] 卡佳霍夫.油層物理基礎(chǔ)[M].北京：石油工業(yè)出版社，1985.

[5] 喻西崇，趙金洲，鄔亞玲，等.深部凝析氣井流入動態(tài)分析研究[J].石油勘探與開發(fā)，2002，29（3）：71-73.Yu Xichong，Zhao Jinzhou，Wu Yaling，et al.A study on deep condensate gas well inflow performance equation[J].Petroleum Exploration and Development，2002，29（3）：71-73.

[6] 康曉東，李相方，郝偉.氣井高速非達(dá)西流動附加壓降計(jì)算公式的修正[J].油氣井測試，2004，13（5）：4-5.Kang Xiaodong，Li Xiangfang，Hao Wei. Pressure drop calculation of non-darcy flow at high velocity for gas well[J].Well Testing，2004，13（5）：4-5.

[7] 康曉東，李相方，程時(shí)清. 考慮流動邊界影響的氣井非達(dá)西效應(yīng)評價(jià)[J].中國石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2006，30（1）：82-85.Kang Xiaodong，Li Xiangfang，ChenShiqing.Evaluation of non-Darcy flow effect in gas well considering influence of flow boundary[J].Journal of China University of Petroleum，2006，30（1）：82-85.

[8] 楊筱璧，李祖友，魯敏蘅，等.高速非達(dá)西流氣井產(chǎn)能方程的新形式[J]. 特種油氣藏，2008，15（5）：74-75，83.Yang Xiaobi，Li Zuyou，Lu Minheng，et al.New deliverability equation for gas wells with high velocity non-Darcy flow[J]. Special Oil & Gas Reservoirs，2008，15（5）：74-75，83.

[9] 陳元千，董寧宇.確定氣井湍流系數(shù)和湍流表皮系數(shù)的新方法[J].斷塊油氣田，2001，8（1）：20-23.Chen Yuanqian，Dong Ningyu.New method of determining turbulence factor and turbulence skin factor[J].Fault-Block Oil&Gas Field，2001，8（1）：20-23.

[10] Noman R，Shrimanker N，Archer J S. Estimation of the coefficient of inertial resistance in high rate gas wells[C].SPE 14207，1985.

[11] Tessem R. High velocity coefficient dependence of rock properties[M]. Norway：Thesis Pet.Inst，NTH，1980：26-30.

[12] Firoozabodi A，Katz D L.An analysis of high velocity gas flow through porous media[C].SPE 6827，1979.

[13] 陳元千，郭二鵬，張楓.確定氣井高速湍流系數(shù)方法的應(yīng)用與對比[J].斷塊油氣田，2008，15（5）：53-55，90.Chen Yuanqian，Guo Erpeng，Zhang Feng. Application and comparison of method for determining high velocity turbulence flow coefficient of gas well[J].Fault-Block Oil&Gas Field，2005，15（5）：53-55，90.

[14] Jones S C.Using the inertial coefficient β to characterize heterogeneity in reservoir rock[C].SPE 16949，1987.

[15] Katz D L，Cornell D. Handbook of natural gas engineering[M].NewYork：McGraw-Hill Book Company Inc.1959：49-50.

[16] Irmay S. On the theoretical derivation of darcy and forchheimer formulas[J]. Transactions，1958，39（4）：702-707.

[17] Geertsma J. Estimating the coefficient of inertial resistance in fluid flow through porous media[J]. SPE 4706，1974.

[18] Liu X，Civan F，Evans R D.Correlation of the non-darcy flow coefficient[C].SPE 95-10-05，1995.

[19] Thauvin F，Mohanty K K. Network modeling of non-Darcy flow through porous media[J].Transport in Porous Media，1998，31：19-37.

[20] Coles M E，Hartman K J.Non-darcy measurements in dry core and the effect of immobile liquid[C]. SPE 39977，1998.

[21] Cooper J W，Wang Xiuli，Mohanty K K.Non-Darcy flowstudies in anisotropic porous media[C].SPE 57755，1999.

[22] Li Dacun，Svec R K，Engler T W，et al. Modelingand simulation of the wafer non-darcy flow experiments[C].SPE 68822，2001.

[23] Janicek，John Daniel，Donald La Verne Katz.Applications of unsteady state gas flow calculations[C]//Flow of Natrual Gas from Reservoirs，Michigan，1955.

[24] Li Dacun，Engler T W.Literature review on correlations of the non-darcy coefficient[C].SPE 70015，2001.