近井地带高速非达西动边界对产能预测的影响

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

1.中国石油大学石油工程教育部重点实验室，北京 昌平102249；2.中海石油（中国）有限公司上海分公司研究院，上海 徐汇200030 3.中国石化上海海洋油气分公司研究院，上海 浦东200120；4.陕西延长石油（集团）有限责任公司研究院，陕西 西安710000

引 言

常规的二项式产能方程中，将气井的泄气半径当成气井高速非达西外边界[1-2]，实际上高速非达西外边界是随着产量变化的，并且只存在于近井地带比较小的范围内，从而会对气井的产能预测造成很大影响。在研究高速非达西时一般用雷诺数表征，比较合理的公式为前苏联学者卡佳霍夫提出的表达式[3-4]，但是目前在应用该公式到气井产能方程中时，一方面将地层平均压力[5-8]看作等于非达西动边界处的压力，而实际上非达西动边界处的压力与储层的物性以及气井工作制度有关不能进行简单的处理；另一方面方程的推导也并不完整，可以对气体的密度进一步简化[3-7]。

1 高速非达西流动边界半径rh 的确定

对于二维径向流，按渗流规律可将渗流区域分为两个区域[8]：（1）达西流动区域，边界为（rh，re）。该区域由于气体渗流速度小，服从达西流动规律。（2）高速非达西流动区域，边界为（rw，rh）。该区域由于气体过流面积减小，加上气体膨胀，气体不服从达西渗流规律。 rh称为非达西动边界。目前的产能方程中通常认为rh=re即整个渗流区域都是高速非达西流动区域。而实际整个区域可以用如图1 与式（1）所示。

图1 地层渗流示意图Fig.1 Schematic of flow

由于在地层中的渗流和管流的相似性，可以用雷诺数值来判断流体在地层中是否服从达西定律[4]。雷诺数有许多不同求法，目前公认比较合理的是前苏联学者卡佳霍夫提出的表达式[3-4]

通过实验获得临界雷诺数为0.2～0.3[3-4]，雷诺数小于临界雷诺数时流体流动服从达西定律，雷诺数大于临界雷诺数时流体流动不服从达西定律。

假设非达西流动边界rh处气体的密度、体积系数和气体的速度分别为

将临界雷诺数以及式（3）～式（5）代入式（2）得到不同流量下的非达西动边界半径，如式（6）所示

从式（6）中可以看出非达西外边界与气井的产量有关，由于气井产能方程为

结合式（6）与式（7）可以建立非达西动边界与气井产能方程的关系，进行迭代求解可以得到高速非达西外边界的大小，如图2 所示。从图中可以看出，随着储层渗透率的增大，高速非达西区域逐渐增大；随着井底流压的增大，生产压差逐渐减小，高速非达西区域减小。并且渗透率在0.8～5.0 mD 变化时，非达西动边界为0.07～1.20，根据式（7）可以发现高速非达西动边界对方程右边第二项系数的影响很大。

2 高速非达西临界渗透率与临界产量

2.1 非达西临界渗透率

由式（6）可知，当地层参数渗透率、孔隙度、气体黏度和地层厚度都不变时，非达西动边界只与气体的产量有关。当气体的产量为绝对无阻流量时即qg=qAOF，最大的非达西动边界为

当井底流压为大气压时，气井的产能方程为

由式（8）、式（9）可知，只有当最大的非达西外半径要大于井筒的井眼即：rmaxh＞rw时，地层中才会有高速达西出现。故存在临界渗透率使得当气井的产能为其无阻流量时，最大高速非达西动边界等于井眼半径即：rmaxh=rw。

结合式（8）、（9）可得临界渗透率为

当地层渗透率大于临界渗透率时，地层中才有可能出现高速非达西渗流。当地层渗透率小于临界渗透率时，无论气井以多少产量生产，高速非达渗流都不会在地层中产生。如图3 所示，随着地层压力的增大，临界渗透率逐渐减小。

图3 临界渗透率与地层压力的关系曲线Fig.3 Curves between formation pressure and critical permeability

2.2 高速非达西临界产量

当渗透率大于KD时，气井工作制度的变化影响了地层中的渗流形态。当气井产量很小时，高速非达西动边界仍可能比rw小，地层中的渗流符合达西定律；当气井产量很大时，高速非达西动边界比rw大，地层中线性流破坏。故存在使得非达西动边界rh等于rw时的产量，定义该气井产量为气井高速非达西临界产量

3 近井地带高速非达西产能方程的建立

由于二项式产能方程中非达西项为[8]

非达西系数为

将式（6）与式（13）代入式（12），得

将式（15）代入气井产能方程中，则产能方程为

其中

得到无阻流量

从式（15）和AD的表达式可以看出，非达西边界的影响可以看成为高速非达西表皮

其中β 是与渗透率、孔隙度相关的函数，其综合表达式可以写成

许多学者通过对单相流体的实验研究得到不同的β 函数关系式，对于有裂缝的盐水或者油气体系，系数c 和d 可以为0，β 仅是与渗透率有关的函数，表1 为不同的学者通过回归得到的不同表达式。

表1 仅与渗透率有关的β 表达式数据表Tab.1 Data of β expressions concerned only with permeability

不同的学者根据实验假设条件下也得到β 与孔隙度、迂曲度等参数的关系，如表2 所示。

表2 与渗透率、孔隙度以及迂曲度有关的β 表达式数据表Tab.2 Data of β expressions concerned with permeability，porosity and capillary tortuosity

从表1 和表2 可以发现，众多学者回归的公式中渗透率都处于分母位置，系数基本大于0.5；孔隙度既有处于分子也有处于分母位置，当其处于分母时系数基本小于1.5；迂曲度处于分子位置。故将系数β 代入式（17）中可以得到表皮的一般表达式

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

另外，由式（15）可知，高速非达西产生的条件为

将不同的β 表达式代入式（17）中可以得到不同的e，d，g 系数，现将常用的式（14）中β 表达式（1）代入式（17），计算并取临界雷诺数NRel=0.3。则非达西的边界表皮为

SD是地层渗透率和孔隙度的函数。其主要反映了非达西区域半径的影响：孔隙度越大，高速非达西半径也越小，高速非达西表皮越小；渗透率越大，高速非达西半径也越大，高速非达西表皮越小。其中SD是该产能方程区别于常规二项式产能方程的地方。

4 气井产能比较

结合式（16）和式（22）绘出常规二项式产能方程、达西线性流动产能方程以及本文推到产能方程在不同渗透率下的IPR 曲线，如图4 所示。

从图4 中可以看出，当渗透率小时，气井临界产量大于气井的无阻流量时，高速非达西渗流在地层不存在，地层中的渗流符合达西定律，故IPR 曲线完全与线性流动时的IPR 曲线重合如图4a 前半段所示。当气井无阻流量大于临界产量如图4a～图4c 所示，随着渗透率的增大，本文产能方程的IPR曲线不再与线性流动产能方程的IPR 曲线重合，越来越向常规二项式产能方程靠近。当渗透率足够大时，本文推导的产能方程的IPR 曲线与常规二项式产能方程IPR 曲线在气井产量大于临界产量时完全重合。所以本文的产能方程是符合实际情况的。

图4 不同渗透率下的IPR 曲线Fig.4 IPR curves at different permeability

5 结 论

（1）低渗气藏中气井高速非达西动边界半径要远小于气井的泄气半径，只存在于井筒附近。并且渗透率越大高速非达西动边界越大；产量越大非达西动边界半径越大。并且当渗透率很低时，高速非达西动边界对系数B 的影响很大，随着动边界的增大系数B 增大，但动边界大于10 倍的井眼半径时，系数B 与常规二项式产能方程的系数B 基本相等。

（2）当气井产量为其无阻流量时，存在使得高速非达西动边界等于井眼半径时的临界渗透率，当储层渗透率大于该渗透率时，储层的高速非达西才可能存在。并且临界渗透率随着地层压力的变化较大，随着地层压力的增大，临界渗透率变小。

（3）近井地带高速非达西的产能方程与常规产能方程相比多了近井地带高速非达西表皮，该表皮只与地层孔隙度、渗透率有关的参数。渗透率越大表皮越小，孔隙度越大表皮也越大。

（4）与常规的二项式产能方程相比，由于本文推导的产能方程是考虑近井地带高速非达西动边界的影响，更符合实际的生产情况，用二项式产能方程预测低渗气藏的产能误差很大，所以本文推导的产能方程对预测气藏的产能方程更为准确。

符号说明

p—压力，MPa；

r—波及半径，m；

µ—气体黏度，mPa·s；

K—渗透率，mD；

υ—渗流速度，m/s；

β—速度系数，m-1；

ρ—气体的密度，kg/m3；

re—泄气半径，m；

rh—高速非达西动边界半径，m；

rw—井眼半径，m；

Re—雷诺数，无因次；

ϕ—孔隙度，无因次；

ρh—非达西动边界处的密度，kg/m3；

γg—气体的相对密度，无因次；

ph—非达西动边界处的压力，MPa；

psc—标况下的压力，MPa；

T—地层温度，K；

Zh—非达西动边界处压力对应的偏差因子，无因次；

Bgh—非达西动边界处压力对应的体积系数，无因次；

Tsc—标况下的温度，K；

υh—非达西动边界处的速度，m/s；

qg—产量，m3/d；

h—储层厚度，m；

Rel—临界雷诺数，无因次；

¯pR—地层平均压力，MPa；

pwf—井底流压，MPa；

rmaxh—最大高速非达西动边界，m；

qAOF—无阻流量，m3/d；

KD—高速非达西临界渗透率，mD；

Z—偏差因子，无因次；

qminh—高速非达西临界产量，m3/d；

Δp—压差，MPa；

D—紊流系数，（m3/d）2；

SD—高速非达西边界表皮，无因次；

a，b，c，d—速度系数表达式的系数，无因次；

τ—迂曲度，无因次。

[1] 李士伦.天然气工程[M].北京：石油工业出版社，2008.

[2] 罗伯特，沃特恩伯格.气藏工程[M].王玉普，译.北京：石油工业出版社，2007.

[3] 杨胜来，魏俊之.油藏物理学[M].北京：石油工业出版社，2004.

[4] 卡佳霍夫.油层物理基础[M].北京：石油工业出版社，1985.

[5] 喻西崇，赵金洲，邬亚玲，等.深部凝析气井流入动态分析研究[J].石油勘探与开发，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] 康晓东，李相方，郝伟.气井高速非达西流动附加压降计算公式的修正[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] 康晓东，李相方，程时清. 考虑流动边界影响的气井非达西效应评价[J].中国石油大学学报：自然科学版，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] 杨筱璧，李祖友，鲁敏蘅，等.高速非达西流气井产能方程的新形式[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] 陈元千，董宁宇.确定气井湍流系数和湍流表皮系数的新方法[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] 陈元千，郭二鹏，张枫.确定气井高速湍流系数方法的应用与对比[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.