引言

页岩气藏开发实践表明，页岩气藏存在多重复杂的运移机理^[1-2]，并且脆性矿物含量较高、天然裂缝发育^[3-4]，加之体积压裂技术的广泛应用使得页岩气藏裂缝形态更加复杂。运移机理的多样性及次生裂缝的复杂性加大了页岩气藏流动特征分析的难度。尽管考虑页岩气藏多重运移机理的耦合模型^[5-8]已相对成熟，但对于次生裂缝形态的考虑较少，大多基于常规的等效物理模型，因而无法体现次生裂缝对气井压力响应特征的影响。目前主要的流动物理模型有：多重介质等效模型^[9-10]、三/五线性流模型^[11]、双翼平板缝模型^[12]、多区复合模型^[13-14]、树状分形模型^[15]、离散裂缝模型^[16-18]。其中前两种模型多采用解析的方法进行求解，计算速度快但缺乏对裂缝形态和次生裂缝的描述；而平板缝模型和多区复合模型多采用半解析的方法进行求解，对气藏流动解析处理，对裂缝或者SRV区进行离散化数值求解，该类模型能有效区分改造区和未改造区，但由于交汇裂缝流量分配问题及边界条件复杂的原因，少见有考虑次生裂缝交汇的情况；离散裂缝模型多采用数值手段进行求解，虽然能详细刻画裂缝，但计算速度慢。随后，贾品等人^[19-20]基于半解析方法，采用星三角变换来处理交汇裂缝的流量分配问题，但该方面的理论并未在页岩气藏进行过应用研究。为此，基于前人的研究^[19-20]，建立考虑次生裂缝形态特征的页岩气藏压裂水平井有限导流模型，并结合Laplace变换、局部坐标转换、叠加原理等数学手段对模型进行半解析求解，以期快速获取压力响应曲线，为页岩气藏压裂水平井压力动态特征分析做指导。

1 物理模型

考虑有限导流次生裂缝的影响对页岩气藏压裂水平井物理模型做如下假设：

(1) 页岩气藏为双重介质气藏，气体在基质中以拟稳态形式流动、天然裂缝中以渗流形式流动，气体解吸满足Langmuir单分子层等温吸附定律^[21]。(2)人工主裂缝与井筒正交等间距分布，多组次生裂缝与主缝成一定角度相交，如图 1所示。主缝与次缝均为有限导流裂缝。(3)气井以恒定产量进行生产，井筒中不存在压降，具有无限导流能力。(4)不考虑重力和页岩储层的压缩性。

2 数学模型

页岩气从气藏到井筒的不稳定流动分为两部分：双重介质页岩气藏中的流动及压裂裂缝中的流动。

2.1 页岩气藏裂缝渗流-基质拟稳态扩散的双重介质模型

气藏中的流动又可分为基质中的拟稳态扩散和天然裂缝中的渗流。

联立基质扩散与裂缝渗流微分方程，并结合Langmuir等温吸附定律，可得到球坐标下的气藏流动微分方程

$ \dfrac{1}{{r_{\rm{D}}^2}}\dfrac{\partial }{{\partial {r_{\rm{D}}}}}\left( {r_{\rm{D}}^2\dfrac{{\partial {{\overline \psi }_{{\rm{fD}}}}}}{{\partial {r_{\rm{D}}}}}} \right) = f\left( u \right){\overline \psi _{{\rm{fD}}}} $

(1)

式中：

$r_{\rm{D}}$—球形基质径向距离，无因次；

$\overline \psi _{{\rm{fD}}}$— Laplace空间拟压力，无因次；

$f(u)$—流动方程系数，无因次；

$f\left(u \right) = \omega u + \dfrac{{\alpha u\lambda \left({1 - \omega } \right)}}{{u + \lambda }}$；

$\omega$—储容比，无因次；

$u$— Laplace变量；

$\lambda$—窜流系数，无因次；

$\alpha$—吸附解吸常数，无因次；

$\alpha = \dfrac{{{V_{\rm{L}}}{p_{\rm{L}}}}}{{\left({{p_{\rm{L}}} + {p_{\rm{f}}}} \right)\left({{p_{\rm{L}}} + {p_{\rm{i}}}} \right)}}\dfrac{{{p_{{\rm{sc}}}}{q_{{\rm{sc}}}}T}}{{{\rm{ \mathsf{ π} }}K_{\rm{f}}h{T_{{\rm{sc}}}}}}\dfrac{{{\mu _{\rm{i}}}{Z_{\rm{i}}}}}{{2{p_{\rm{i}}}}}$；

$V$—基质中气体平均浓度，sm$^3$/m$^3$；

$p$—压力，MPa；

$q$—页岩气产量，m$^3$/s；

$T$—温度，K；

$\mu$—气体黏度，Pa$\cdot$s；

$K_{\rm{f}}$—裂缝渗透率，m$^2$；

$h$—气藏厚度，m；

$Z$—气体偏差因子，无因次。

下标

f—天然裂缝；

D—无因次量；

sc—地面状态；

i—初始状态；

L— langmuir常量。

上述微分方程式(1)通过引入吸附解吸常数$\alpha$来反映页岩气吸附解吸的特性，其形式与Ozkan提出的点源函数^[10]基本方程相似，仅$f(u)$的表达式不同。由此基于Laplace空间源函数可求得顶底封闭侧向无限大页岩气藏的连续线源解

$ {\rm{\Delta }}{\overline \psi _{\rm{f}}} = \dfrac{{{p_{{\rm{sc}}}}T}}{{{T_{{\rm{sc}}}}}}\dfrac{{\overline q (u)}}{{{\rm{ \mathsf{ π} }}{K_{\rm{f}}}L{h_{\rm{D}}}}}{{{K}}_0}\left( {{R_{\rm{D}}}\sqrt {f\left( u \right)} } \right) $

(2)

式中：

$\overline q (u)$— Laplace空间下点源流量，m$^3$/s；

$L$—参考长度，m；

$K_0$()—零阶贝塞尔函数；

$R_{\rm{D}}$—源点到场点的距离，无因次。

气藏流体向人工裂缝中流动时，相当于多个足够小的微元段同时生产，因此，各微元段对气藏中某一点的压力叠加响应可通过对式(2)积分获得。另外，对于存在倾角的裂缝，采用局部坐标变换，将斜线积分转化为一维积分

$ \overline \psi _{{\rm{fD}}}^{}\left( {{x_{{\rm{D}}p, q}}, {y_{{\rm{D}}p, q}}} \right) = \sum\limits_i {\sum\limits_j {} } {\overline q _{{\rm{D}}i, j}}\left( u \right)\int\limits_{} {{{{K}}_0}} \\ \left[ {\sqrt {f\left( u \right)} \sqrt {{{\left( {{x_{{\rm{D}}p, q}} - {x_{{\rm{wD}}}}} \right)}^2} + {{\left( {{y_{{\rm{D}}p, q}} - {x_{{\rm{wD}}}}\tan \theta } \right)}^2}} } \right]\sqrt {1 + {{\tan }^2}\theta } {\rm{d}}{x_{{\rm{wD}}}} $

(3)

式中：

$\left({{x_{{\rm{D}}p, q}}, {y_{{\rm{D}}p, q}}} \right)$—局部坐标系下场点离散段($p, q$)坐标；

$x_{{\rm{wD}}}$—无因次源点$x$坐标；

$p$，$q$—场点主缝编号和次缝编号；

$i$，$j$—源点主缝编号和次缝编号；

$\theta$—主缝与次缝夹角，(°)。

2.2 有限导流裂缝流动模型

相比于无限导流裂缝，有限导流裂缝中存在压降，根据质量守恒得到

$ \dfrac{{{\partial ^2}{{\overline \psi }_{{\rm{flD}}}}}}{{\partial l_{\rm{D}}^2}} - \dfrac{{2{\rm{ \mathsf{ π} }}}}{{{R_{{\rm{fD}}}}}}{\overline q _{{\rm{flD}}}}\left( {{l_{\rm{D}}}, u} \right) = \dfrac{u}{{{\eta _{\rm{D}}}}}{\overline \psi _{{\rm{flD}}}} $

(4)

式中：

$l_{\rm{D}}$—裂缝长度，无因次；

$R_{\rm{fD}}$—裂缝导流能力，无因次；

$\eta_{\rm{D}}$—导压系数，无因次。

下标

fl—有限导流裂缝。

对式(4)进行差分格式转换

$ {A_N}{\overline \psi _{{\rm{WD}}\left( N \right)}} + {B_N}{\overline \psi _{{\rm{flD}}\left( N \right)}} + {C_N}{\overline \psi _{{\rm{flD}}\left( {N + 1} \right)}} + {D_N}{\overline \psi _{{\rm{flD}}\left( {N - 1} \right)}} - {E_N}{\overline q _{{\rm{flD}}\left( N \right)}} = 0 $

(5)

式中：

$A_N$，$B_N$，$C_N$，$D_N$，$E_N$—单元$N$的方程系数，与单元位置有关；

$N$—离散段编号；

${\overline \psi _{{\rm{wD}}}}$— Laplace空间下的井底拟压力，无因次。

由于受到缝网内外边界条件及流量交汇分配的影响，如图 2所示。在不同位置处各离散段的差分方程格式不同，即对应系数$A_N$，$B_N$，$C_N$，$D_N$，$E_N$不同^[22]。

(1) 内部裂缝单元

$ \left( { - \dfrac{{{R_{{\rm{fD}}(N)}}}}{{{\rm{ \mathsf{ π} }}\Delta l_{\rm{D}}^2}} - \dfrac{{{R_{{\rm{fD}}(N)}}u}}{{2{\rm{ \mathsf{ π} }}{\eta _{\rm{D}}}}}} \right){\overline \psi _{{\rm{flD}}\left( N \right)}} \\ + \dfrac{{{R_{{\rm{fD}}(N)}}}}{{2{\rm{ \mathsf{ π} }}\Delta l_{\rm{D}}^2}}{\overline \psi _{{\rm{flD}}\left( {N + 1} \right)}} + \dfrac{{{R_{{\rm{fD}}(N)}}}}{{2{\rm{ \mathsf{ π} }}\Delta l_{\rm{D}}^2}}{\overline \psi _{{\rm{flD}}\left( {N - 1} \right)}} - {\overline q _{{\rm{flD}}\left( N \right)}} = 0 $

(6)

(2) 与井筒相邻单元(如图 2中离散段1)

$ \left( { - \dfrac{{3{R_{{\rm{fD}}(N)}}}}{{2{\rm{ \mathsf{ π} }}\Delta l_{\rm{D}}^2}} - \dfrac{{{R_{{\rm{fD}}(N)}}u}}{{2{\rm{ \mathsf{ π} }}{\eta _{\rm{D}}}}}} \right){\overline \psi _{{\rm{flD}}\left( N \right)}} + \\ \dfrac{{{R_{{\rm{fD}}(N)}}}}{{2{\rm{ \mathsf{ π} }}\Delta l_{\rm{D}}^2}}{\overline \psi _{{\rm{flD}}\left( {N + 1} \right)}} + \dfrac{{{R_{{\rm{fD}}(N)}}}}{{{\rm{ \mathsf{ π} }}\Delta l_{\rm{D}}^2}}{\overline \psi _{{\rm{wD}}}} - {\overline q _{{\rm{flD(}}N)}} = 0 $

(7)

(3) 封闭端裂缝单元(如图 2中离散段4、5、8)

$ \left( { - \dfrac{{{R_{{\rm{fD}}(N)}}}}{{2{\rm{ \mathsf{ π} }}\Delta l_{\rm{D}}^2}} - \dfrac{{{R_{{\rm{fD}}(N)}}u}}{{2{\rm{ \mathsf{ π} }}{\eta _{\rm{D}}}}}} \right){\overline \psi _{{\rm{flD}}\left( N \right)}} + \dfrac{{{R_{{\rm{fD}}(N)}}}}{{2{\rm{ \mathsf{ π} }}\Delta l_{\rm{D}}^2}}{\overline \psi _{{\rm{flD}}\left( {N - 1} \right)}} - {\overline q _{{\rm{flD}}\left( N \right)}} = 0 $

(8)

(4) 交汇点附近单元(如图 2中离散段2、3、6、7)

$ \dfrac{{{R_{{\rm{fD}}(N)}}}}{{2{\rm{ \mathsf{ π} }}}}\sum\limits_{i = 1}^N {\dfrac{{{T_{{\rm{D}}i, 0}}{T_{{\rm{D}}N, 0}}}}{{\sum\limits_{k = 1}^N {{T_{{\rm{D}}k, 0}}} }}} {\overline \psi _{{\rm{flD}}(i)}}\!-\! \dfrac{{{R_{{\rm{fD}}(N)}}}}{{2{\rm{ \mathsf{ π} }}}} \\ \left( {\sum\limits_{i = 1}^N {\dfrac{{{T_{{\rm{D}}i, 0}}{T_{{\rm{D}}N, 0}}}}{{\sum\limits_{k = 1}^N {{T_{{\rm{D}}k, 0}}} }}}\!+\!\dfrac{u}{{{\eta _{\rm{D}}}}}\!+\!{T_{{\rm{D}}N + 1, N}}} \right) \\ {\overline \psi _{{\rm{flD}}(N)}}\!+\! \dfrac{{{R_{{\rm{fD}}(N)}}}}{{2{\rm{ \mathsf{ π} }}}}{T_{{\rm{D}}N\!+\! 1, N}}{\overline \psi _{{\rm{flD}}(N + 1)}}\!-\!{\overline q _{{\rm{flD}}\left( N \right)}}\!=\!0 $

(9)

式(9)中，${T_{{\rm{D}}N + 1, N}}$为离散段$N$与$N+1$之间的传导率，${T_{{\rm{D}}N, 0}}$为与交点相邻节点的传导率

$ {T_{{\rm{D}}N + 1, N}} = \dfrac{1}{{\Delta {l_{{\rm{D}}(N)}}\left( {\dfrac{{{l_{{\rm{D}}(N + 1)}}}}{2} + \dfrac{{{l_{{\rm{D}}(N)}}}}{2}} \right)}} = \dfrac{1}{2}{T_{{\rm{D}}N, 0}} $

(10)

综上，式(6)$\sim$式(9)分别为人工裂缝中不同位置处的流动差分控制方程。

2.3 模型求解

综合气藏流动方程和裂缝流动方程，若将所有裂缝离散为$n$段，则未知数共有$2n+1$个，包括$n$个位置处的线密度流量${\overline q _{{\rm{flD}}\left(N \right)}}$和拟压力${\overline \psi _{{\rm{flD}}\left(N \right)}}$以及1个未知井底拟压力${\overline \psi _{{\rm{wD}}}}$。

由气藏流动方程(3)可得到$n$个方程，由有限导流裂缝流动方程式(6)$\sim$式(9)可得到$n$个方程，因此还需要一个方程。考虑在Laplace空间下各压裂段流量之和满足归一化方程

$ \sum\limits_{i = 1}^{{N_{\rm{w}}}} { - {{\overline \psi }_{{\rm{flD}}\left( i \right)}} + {N_{\rm{w}}}} {\overline \psi _{{\rm{wD}}}} = \dfrac{{{\rm{ \mathsf{ π} }}\Delta {l_{\rm{D}}}}}{{2u{R_{{\rm{fD}}}}}} $

(11)

式中：

$N_{\rm{w}}$—与井筒相连的离散单元的个数，无因次。

联立以上$2n+1$个方程即可线性求解。最后根据Duhamel原理和Stehfest数值反演方法，可得到实空间下考虑井筒储集效应和表皮系数的井底压力解，由此可绘制出典型试井曲线。

3 典型曲线及压力动态特征分析

图 3为计算得到的考虑次生裂缝影响的页岩气藏压裂水平井典型图版，根据曲线特征主要划分为以下10个流动阶段。

Ⅰ井储阶段：该阶段在形态上与常规曲线一致，表现为重合的斜率为1的两条曲线。Ⅱ井储后的过渡流阶段：表现为与常规曲线类似的“拱形”。Ⅲ双线性流动阶段：该阶段拟压力导数曲线表现为斜率为“1/4”的直线。Ⅳ线性流阶段：双对数曲线表现为两条斜率为1/2的平行线，此时各裂缝段间互不干扰。Ⅴ裂缝干扰流动阶段：各压裂段内，主缝与次生相互干扰。Ⅵ早期第一径向流阶段：各压裂段间互不干扰，段内主缝与其相沟通的次缝形成一个整体进行协调生产。表现出拟压力导数曲线值为1/2$m$的特点($m$—压裂段数)。Ⅶ椭圆流阶段：该阶段相邻两压裂段之间开始干扰。Ⅷ天然裂缝系统的拟径向流阶段：改造区域裂缝对曲线的影响已经结束。表现出拟压力导数曲线值为1/2的特点。Ⅸ基质向天然裂缝的窜流段：基岩中的气体解吸脱附并以拟稳态的形式向天然裂缝进行扩散补给，在双对数曲线上表现出双重介质的“凹子”特征。Ⅹ晚期总系统的拟径向流段：在双对数曲线上表现出拟压力导数曲线值为1/2的水平直线特征。

本文仅对反映页岩气藏特征和裂缝形态的参数进行影响因素分析。

3.1 吸附解吸常数$\alpha$

模型引入了表征页岩气吸附解吸能力的特征变量$\alpha$，图 4为该特征变量对曲线形态的影响。从图 4中可以看出，吸附解吸常数$\alpha$越小，阶段Ⅸ的“凹子”就越浅，反映出基质的吸附能力越强，相同条件下，该阶段获得补给的能力就越弱，因而拟压力曲线更高，持续时间越短。

3.2 主缝与次缝夹角$\theta$

分别绘制出主缝与次缝夹角$\theta$为10°、30°、50°、70°及90°时的双对数曲线，如图 5所示。从图 5可知，主缝与次缝夹角$\theta$仅对早期线性流和裂缝干扰流阶段有影响。在次生裂缝长度相同的情况下，主缝与次缝夹角越接近90°，表明次生裂缝的横向延伸范围越大，改造区域越大，因而双对数曲线在该对应阶段内的曲线越靠下。主缝与次缝夹角越接近0°，主缝与次缝的干扰发生越快，因而双对数曲线越靠上。

3.3 次生裂缝组数$m_{\rm{f}}$

次生裂缝组数$m_{\rm{f}}$，反映了沿主缝延伸方向上次生裂缝的密集程度。对比次生裂缝组数$m_{\rm{f}}$取1、3、5时的双对数曲线(图 6)，次生裂缝组数$m_{\rm{f}}$主要影响井储过渡流到裂缝干扰流动段的曲线形态，次生裂缝组数$m_{\rm{f}}$越少，该阶段内的双对数曲线越靠上，但其对第一拟径向流之后的曲线形态的影响不大。

3.4 裂缝导流能力$R_{\rm{fD}}$

裂缝导流能力$R_{\rm{fD}}$主要影响过渡流、双线性流、线性流及裂缝干扰流阶段的曲线形态，如图 7所示。

裂缝导流能力$R_{\rm{fD}}$越大，过渡段出现越早，双线性流阶段持续时间越短，而线性流越早出现。当裂缝导流能力$R_{\rm{fD}}$大于一定值时双线性流阶段消失，与无限导流情况下的曲线形态一致。

3.5 无因次导压系数$\eta_{\rm{fD }}$

$\eta_{\rm{fD}}$对双对数曲线形态的影响如图 8所示。从图 8中可以看出，当$\eta_{\rm{fD}}$很小时，拟压力导数曲线在过渡段之后出现了一个“凹子”，且$\eta_{\rm{fD}}$越小“凹子”越深，$\eta_{\rm{fD}}$反映出主裂缝和次生裂缝中流体的压缩性与流动能力的大小，当$\eta_{\rm{fD}}$较大时，即裂缝内流体的压缩性可忽略不计，主裂缝与次生裂缝渗透率远大于天然裂缝渗透率时，“凹子”消失。由于页岩气系统中的主裂缝与次生裂缝体积相比于整个页岩气藏体积很小，故压裂裂缝中由于气体弹性膨胀所引起的流量可忽略不计，此时$\eta_{\rm{fD}}$无穷大。

4 结论

(1) 建立了考虑次生裂缝影响的有限导流裂缝流动模型，利用半解析的方法，有效解决了有限导流裂缝交点的流量分配问题。

(2) 计算得到了井底压力的解，并绘制出典型曲线，共可分为10个流动阶段。

(3) 通过引入吸附解吸常数$a$来表征页岩气吸附解吸能力，$a$越大，阶段Ⅸ的“凹子”越深。

(4) 引入了表征次生裂缝的相关参数：次生裂缝组数$m_{\rm{f}}$越多，角度$\theta$越接近90°，拟压力导数曲线越低。裂缝导流能力$R_{\rm{fD}}$越大，过渡段出现越早，双线性流阶段持续时间越短，而线性流越早出现。