引言

根据内部孔隙结构的分布情况，碳酸盐岩储层可分为裂缝型、缝洞型、基质孔隙型、溶洞型等多种类型^[1-2]。其中，缝洞型储层最为复杂，一般具有岩石基质、裂缝和溶洞3种孔隙类型，3者同时共存于储层中^[3-4]。其中，岩石基质孔隙度和渗透率非常低^[5]，可归类为低储低渗单元，基本不参与流体流动；裂缝开度较大，渗透率较高，但相对于溶洞来说储集能力较小，可归类为低储高渗单元，是流体流动的主要通道；溶洞体积巨大，可归类为高储高渗单元，是流体的主要储存场所^[6]。

国内外学者针对此类复杂储层引入了三重孔隙的概念，建立了三重孔隙模型^[7-12]和三孔双渗模型^[13]，并在多重介质的基础上进行了试井理论研究^[14-16]。此外，部分学者还以管流^[17-18]、渗-管流耦合^[19-20]等理论为基础进行了相关研究。但是，三重介质理论属于微观非均质性理论，无法表征储层内部的大尺度宏观缝洞储集体^[21-23]；且储层中大开度裂缝内的流体流动特征与管流特征并不相符^[24]，而是更接近于平板流特征。因此，建立以裂缝平板流与宏观非均质性理论^[25-26]为基础的储层模型，更加切合缝洞型储层的客观特性，对于串珠状缝洞型油气藏的动态分析与试井解释，深入认识缝洞连通情况和缝洞参数有一定的理论意义。

1 数学模型建立

缝洞型碳酸盐岩储层内部缝洞连通形式各异，溶洞与裂缝之间串并联关系复杂多变。其中，“串珠状”缝洞型储层是最为常见的连通形式之一^[27-30]。所以，本次研究中主要针对串珠状缝洞型储层进行数学模型的建立。

串珠状缝洞型储层数学模型如图 1所示，图中串珠状缝洞储集体共由若干个溶洞与若干条裂缝串联而成，生产井位于裂缝之上。缝洞系统编号原则：对于各溶洞来说，将井筒编为1号溶洞，与井筒直接相连的溶洞为2号溶洞，依次串联的溶洞分别编为3、4、5、……、N号溶洞；对于裂缝来说，连接于x号溶洞与y号溶洞之间的裂缝编为x_y号，如连接于3号与2号溶洞之间的裂缝编号为3_2。

为简化研究，对模型作出如下假设条件：

(1) 模型内流体为单相非理想可压缩气体，流动为等温过程，满足流体力学原理。

(2) 由于岩石基质相对于缝洞储集体的低渗低储特性，故忽略流其内部的流体流动。

(3) 边部裂缝中有一口生产气井，并以定产量的生产制度进行产气。

(4) 由于溶洞内部完全由气体充填，不符合渗流规律，故将溶洞视为一个等势体^{[6, 21, 31-32]}，即其内部压力处处相等。

(5) 考虑井筒储集效应，忽略缝洞储集体的压缩性、井筒压降和井底表皮效应。

1.1 裂缝内气体流动方程

由于缝洞型碳酸盐岩储层中的大型裂缝开度较大^[21]，可达毫米级别。故普通的达西渗流公式^[33]与立方定律^[34]不再适用，需要通过流体力学连续性方程、能量方程以及动量方程^[35]，推导出大开度的平板流公式。

假设气体流过某一裂缝，裂缝长度为L，开度为d，高度为H，裂缝横截面面积为$A=Hd$，湿周为$C=2(H+d)$。裂缝前后两端压力分别为$p_1$和$p_2$，如图 2所示。

沿着裂缝方向的气体受力应该满足

${p_1}A - {p_2}A - \rho {\rm{g}}LA\sin \alpha - {\bar \tau _0}CL = 0$

(1)

式中： ${p_1}，{p_2}$—裂缝两端压力，Pa；

A—裂缝截面积，m²；

$\rho$—气体密度，kg/m³；

g—重力加速度，g=9.8 m/s²；

L—裂缝长度，m；

$\alpha$—裂缝倾斜角度，(°)；

C—裂缝过流截面湿周，m；

${\bar \tau _0}$—裂缝壁面对气体的平均切向力，Pa。

裂缝壁面对气体的平均切向力计算公式为

${\bar \tau _0} = \dfrac{{\int_0^C {\tau {\rm{d}}C} }}{C}$

(2)

式中： $\tau$—裂缝壁面对气体的切向力，Pa。

由图 2可知，$\sin \alpha = \left( {{z_2} - {z_1}} \right)/L$，代入式(1)后，两端同时除以$\rho {\rm g}A$，得

$\dfrac{{{p_1}}}{{\rho {\rm{g}}}} - \dfrac{{{p_2}}}{{\rho {\rm{g}}}} - {z_2} + {z_1} = {\bar \tau _0}\dfrac{{CL}}{{\rho {\rm{g}}A}}$

(3)

式中： $z_1{，}z_2$—裂缝两端高度，m。

由伯努利原理，该流动的摩阻水头为

${h_{\rm{f}}} = \left( {{z_1} + \dfrac{{{p_1}}}{{\rho {\rm g}}}} \right) - \left( {{z_2} + \dfrac{{{p_2}}}{{\rho {\rm g}}}} \right)$

(4)

式中： ${h_{\rm{f}}}$—摩阻水头，m。

结合式(3)，可得到

${h_{\rm{f}}} = {\bar \tau _0}\dfrac{{CL}}{{\rho {\rm g}A}}$

(5)

通过量纲分析^[35]，平均切向应力可表示为

${\bar \tau _0} = {C_{\rm{f}}}\rho \dfrac{{{v^2}}}{2}$

(6)

式中： v—裂缝内气体流速，m/s；

${C_{\rm{f}}}$—气体摩擦系数，无因次。

将式(6)代入式(5)，有

${h_{\rm{f}}} = {C_{\rm{f}}}\dfrac{{CL}}{A}\dfrac{{{v^2}}}{{2{\rm g}}}$

(7)

忽略气体在流动中由于摩擦作用导致的温度上升，将气体在裂缝中的流动视为等温过程，根据流体力学中的能量方程^[35]

$\dfrac{{{\rm{d}}p}}{\rho } + {\rm{gd}}z + v{\rm{d}}v + \dfrac{{\tau C}}{{\rho A}}{\rm{d}}s = 0$

(8)

式中： p—气体压力，Pa；

s—移动距离，m。

结合式(2)、式(5)，有

$\dfrac{{\tau C}}{{\rho A}} = {\rm{g}}\dfrac{{{h_{\rm{f}}}}}{L}$

(9)

将式(7)代入式(9)，得

$\dfrac{{\tau C}}{{\rho A}} = \dfrac{{{C_{\rm{f}}}}}{{A/C}}\dfrac{{{v^2}}}{2}$

(10)

式(10)代入式(8)，得

$\dfrac{{{\rm{d}}p}}{\rho } + {\rm{gd}}z + v{\rm{d}}v + \dfrac{{{C_{\rm{f}}}}}{{A/C}}\dfrac{{{v^2}}}{2}{\rm{d}}s = 0$

(11)

假设该裂缝中气体的质量流量为$\dot m$，根据连续性方程

$\dot m = \rho Av$

(12)

式中： $\dot m$—气体质量流量，kg/s。

同时，气体状态方程为

$\rho = p/Z{\rm R'}T$

(13)

式中： Z—气体体积偏差系数，无因次；

$\rm R'$—气体常数，Pa·m³/(kg·K)；

T—温度，K。

结合气体状态方程及式(12)、式(13)，有

$v = \dfrac{{\dot mZ{\rm R'}T}}{{pA}}$

(14)

将式(14)代入式(11)，假设裂缝倾角不大，并且气体密度较小，忽略式中的重力项，整理，得

$ - \left( {\dfrac{{2{A^2}}}{{{{\dot m}^2}Z{\rm{R'}}T}}} \right)p{\rm{d}}p = \dfrac{{{C_{\rm{f}}}}}{{A/C}}{\rm{d}}s + \dfrac{{2{\rm{d}}v}}{v}$

(15)

将式(15)沿裂缝长度方向积分，整理，得到气体在裂缝中流动的质量流量方程

$\dot m = A\sqrt {\dfrac{{ {p_1^2 - p_2^2} }}{{Z{\rm{R'}}T\left( {{C_{\rm{f}}}\dfrac{L}{{A/C}} + 2\ln \dfrac{{{p_1}}}{{{p_2}}}} \right)}}}$

(16)

代入裂缝截面积A和湿周C的表达式，整理，可得

$\dot m = Hd\sqrt {\dfrac{{\left( {p_1^2 - p_2^2} \right)}}{{Z{\rm{R'}}T\left[{{C_{\rm{f}}}\dfrac{{2L(H + d)}}{{Hd}} + 2\ln \dfrac{{{p_1}}}{{{p_2}}}} \right]}}}$

(17)

式中： H—裂缝高度，m；

d—裂缝开度，m。

换算成标况下气体在裂缝中的体积流量公式为

${q_{{\rm{sc}}}} = \dfrac{{Hd}}{{{\rho _{{\rm{sc}}}}}}\sqrt {\dfrac{{\left( {p_1^2 - p_2^2} \right)}}{{Z{\rm{R'}}T\left[{{C_{\rm{f}}}\dfrac{{2L(H + d)}}{{Hd}} + 2\ln \dfrac{{{p_1}}}{{{p_2}}}} \right]}}}$

(18)

式中： ${q_{{\rm{sc}}}} $—标况下气体表体积流量，m³/s；

${\rho_{{\rm{sc}}}} $—标况下气体密度，kg/m³。

与流体力学中的管流类似，气体摩擦系数可表示为^[35]

${C_{\rm{f}}} = f/4$

(19)

式中：f—气体摩阻系数，无因次。

气体摩阻系数可通过迭代法求解，具体算法为

${\mathop{ Re}\nolimits} = \dfrac{{\rho vd}}{\mu }$

(20)

式中： $\mu$—黏度，Pa·s。

$Re$—雷诺数，无因次。

$v = \dfrac{{202000ZT{q_{{\rm{sc}}}}}}{{273\left( {{p_1} + {p_2}} \right)A}}$

(21)

当${Re}\leqslant 2000$时

$f = \dfrac{{64}}{{{\mathop{ Re}\nolimits} }}$

(22)

当${Re} > 2000$时

$\dfrac{1}{{\sqrt f }} = 1.14 - 2\lg \left (\dfrac{E}{d} + \dfrac{{21.25}}{{{{{\mathop{Re}\nolimits} }^{0.9}}}}\right )$

(23)

式中： E—壁面绝对粗糙度(一般取0.000 016 m)，m。

1.2 溶洞内气体物质守恒方程

气体物质守恒方程即为气体状态方程^[36]，广泛应用于天然气工程中。它描述了一定物质量气体的压力、温度和体积之间的关系。

若假设气体分子间无作用力、分子体积与总体积相比小到可忽略不计，以及分子间、分子与容器间碰撞完全是弹性碰撞而无内能损耗，那么可以得到理想气体状态方程^[37]，它主要用于低压条件。实际应用中主要涉及高压条件下的实际气体，描述实际气体的状态方程很多。而本次研究选用最常用的Clapeyron方程^[38]

$pV = nZ{\rm R}T$

(24)

式中： V—气体体积，m³；

n—气体物质的量，mol；

R—普适气体常数，R = 8.314 J/(k·mol)。

1.3 串珠状缝洞模型气体流动方程

根据图 1所示数学模型与缝洞单元编号规则，并考虑到裂缝高度远远大于裂缝开度($H \gg d$) ，联立裂缝流动方程与质量守恒方程，得出串珠状缝洞模型在定产量产气时流动方程

$q_{1\_0} = {常数}$

(25)

$\begin{align} & {{q}_{{{i}\_{(}}i-1)}}=\frac{{{H}_{{{i}\_{(}}i-1)}}{{d}_{{{i}\_{(}}i-1)}}}{{{\rho }_{\text{sc}}}}\sqrt{\frac{p_{i}^{2}-p_{i-1}^{2}}{{{{\bar{Z}}}_{{{i}\_{(}}i-1)}}{R}'T\left[ {{f}_{{{i}\_{(}}i-1)}}\frac{{{L}_{{{i}\_{(}}i-1)}}}{2{{d}_{{{i}\_{(}}i-1)}}}+2\ln \frac{{{p}_{i}}}{{{p}_{i-1}}} \right]}} \\ & \left( i=2,3,\cdots ,N \right) \\ \end{align}$

(26)

式中： $q_{1\_0}$—气体产量(标况)，m³/s；

${q_{{{i\_(i}} - {{1)}}}}$—${{{i\_(i}} - {{1)}}}$号裂缝内的气体流量(标况)，m³/s；

${H_{{{i\_(i}} - {{1)}}}}$—${{{i\_(i}} - {{1)}}}$号裂缝的高度，m； $p_{i}$—i号溶洞内压力，Pa；

${{\overline Z }_{{{i\_(i}} - {{1)}}}}$—${{{i\_(i}} - {{1)}}}$号裂缝内气体的偏差因子，无因次；

${f_{{{i\_(i}} - {{1)}}}}$—${{{i\_(i}} - {{1)}}}$号裂缝内气体的摩阻系数，无因次；

${L_{{{i\_(i}} - {{1)}}}}$—${{{i\_(i}} - {{1)}}}$号裂缝的长度，m；

${d_{{{i\_(i}} - {{1)}}}}$—${{{i\_(i}} - {{1)}}}$号裂缝的开度，m；

N—溶洞数。

根据气体状态方程，有

$ {p_{i}}{V_{i}} = {n_{i}}{Z_{i}}{\rm R}T{，}{\kern 10pt} \left (i = 2,3,\cdots,N \right)$

(27)

式中： $V_{i}$—i号溶洞体积，m³；

$n_{i}$—i号溶洞内气体物质的量数，mol；

$Z_{i}$—i号溶洞内气体偏差因子，无因次。

在给定时间dt内，各溶洞内气体物质摩尔量的变化为

${\rm{d}}{n_i} = \dfrac{{ - {q_{i\_(i - 1)}}{\rm{d}}t}}{{22.4/1000}}{，}{\kern 10pt} \left (i = 2,3,\cdots,N \right)$

(28)

联立式(25)~式(28)，即得到了描述串珠状缝洞模型气体流动的方程组。

2 数值计算

假设某一串珠状缝洞型碳酸盐岩储层中充满天然气(甲烷)，初始压力、温度、缝洞参数、产量、开井时间等基本参数如表 1所示。

根据表 1数据，利用牛顿迭代法与差分法^[39]求解珠状缝洞模型气体流动方程，可计算出各溶洞内压力变化与各裂缝中气体流量变化。

3 结果分析

以4溶洞模型为例，代入表 1中各参数，利用数值计算方法得出各溶洞压力数据，绘制出各个溶洞压力导数与时间的半对数曲线(图 3)。由图 3可见，各溶洞压力导数半对数曲线呈现出先上升后下降的趋势，其中距离井筒越近的溶洞压力导数的上升和下降都更快，转折的时间点也更早。这是由于靠近井筒的溶洞初期的产气速度最大，导致压力下降很快；随着其压力的不断下降，后部串联的溶洞开始大量向其补充气体，促使压力下降变缓，压力导数下降。随着时间的推移，晚期各溶洞流量趋于平衡，压力导数趋于一致，说明压力下降速度也趋于一致。

利用数值计算方法得出该模型中的井底压力数据，并绘制出井底压力双对数曲线，如图 4所示。

从图 4可以看出，压力导数曲线可以分为4个阶段，第一段与常规试井相同，为井筒储集反应段；第二段为裂缝反应阶段，由于裂缝相对于溶洞来说导流能力有限，类似于封闭边界的情况，导数曲线上升；第三段为溶洞反应阶段，当溶洞产量随着井底压力的减小迅速升高时，井筒得到大量气体补充，压力下降减缓，导数出现下降；第四段为边界反应阶段，由于忽略基质渗流，故可视为封闭边界，导数曲线上升。

随后，改变模型中不同的缝洞参数，计算不同模型定产量生产时的井底压力数据，并绘制出井底压力双对数曲线，进行对比分析。

3.1 溶洞数量的影响

仅改变模型溶洞数量，利用数值计算方法得出各溶洞压力数据，绘制出不同溶洞数量时井底压力双对数曲线，如图 5所示。

从图 5中可以看出，溶洞个数对井底压力导数曲线有明显影响：当溶洞个数越多时，压降曲线晚期与压力导数曲线第四阶段，即边界反应阶段的上升速度越小，说明溶洞数量越多，气体的补充能力越强，封闭边界反应越弱。

3.2 溶洞体积与排列顺序的影响

溶洞个数设定为3，仅改变模型溶洞排列顺序和体积，绘制出不同溶洞排列顺序和不同溶洞体积时井底压力双对数曲线，如图 6、图 7所示。

由图 6可以看出，溶洞体积越大，压降曲线晚期上升速度越小；导数曲线第三阶段，即溶洞反应阶段下降速度与幅度越大。这说明溶洞体积越大，对井筒的气体补充作用越明显。

由图 7可以看出，若靠近井筒的溶洞体积越大，压降曲线中期上升速度越小，晚期上升速度有增大趋势；导数曲线第三阶段下降特征越明显，若靠近井筒的溶洞体积过小，有可能出现导数曲线下降阶段缺失的现象。两个模型的导数曲线在第四阶段晚期有趋向一致的趋势。

3.3 裂缝开度与排列顺序的影响

溶洞个数设定为3，仅改变模型裂缝开度和排列顺序，绘制出不同裂缝开度和裂缝排列顺序时井底压力双对数曲线，如图 8、图 9所示。

裂缝开度增大，压降曲线整体下移；压力导数曲线第二阶段，即裂缝反应阶段上升速度与幅度越小。这说明裂缝开度越大，溶洞对井筒的气体补充作用越明显，封闭边界效应越弱(图 8)。

靠近井筒的裂缝开度越大，压降曲线整体越偏小，中期上升速度越小，晚期上升速度有增大趋势；压力导数曲线第二阶段上升幅度越小，第三、第四阶段开始时间越早，并且两者曲线在晚期有趋向一致的趋势(图 9)。

4 结论

(1) 建立了串珠状缝洞型储层中气体流动的数学模型，并通过数值方法求解出了定产量生产时储层内压力变化规律。

(2) 模型中各溶洞压力导数半对数曲线呈现出先上升后下降的趋势；初期阶段，由于靠近井筒的溶洞对产气量贡献更大，所以其压力导数的上升和下降时间都更早。末期阶段，各溶洞流量趋于平衡，压力导数趋于一致。

(3) 压力导数双对数曲线可以分为4段：井筒储集反应段、裂缝反应阶段、溶洞反应阶段、边界反应阶段。

(4) 模型中缝洞参数的变化对井底压力双对数曲线的对应阶段特征都有比较明显的影响。所以，理论上可以通过对比井底压力双曲线变化特征，判断模型的参数情况。