引言

长井段水平井+缝网压裂是实现页岩气藏商业开发的关键技术。然而，目前对于页岩气藏压裂水平井渗流模型的研究仍十分欠缺，成为制约页岩气高效开发的瓶颈之一。建立页岩气的渗流模型比常规、致密气藏更为困难，其原因主要在于：(1)页岩气的赋存状态特殊(见图 1)，大量的页岩气以吸附态存在于基质孔隙表面^[1-2]，需考虑吸附气解吸规律的定量表征。(2)页岩气藏孔隙类型多样且尺度跨度很大：发育有纳米级基质孔隙和微米级天然裂缝，还存在压裂所形成的微米至毫米级的次生裂缝和具有更大尺度的水力主裂缝。继而，需考虑页岩气多尺度渗流的特点：气体在微裂缝中主要为黏性流，但在基质孔隙中的运移主要受滑脱、扩散等微观渗流机理的控制。另外，储层压裂后产生复杂缝网使得各尺度空间内的气体流动规律十分复杂。

调研发现，线性耦合方法在很大程度上降低了渗流模型建立和求解的难度，已有很多学者基于该方法建立了页岩气藏多级压裂水平井渗流模型，但这些模型还存在以下几点不足：(1)采用双重介质模型处理页岩基质的方法过于理想^[3-7]，研究表明，原始状态下极为致密的储层发育有大量离散分布的孤立不连通孔隙-微裂缝^[8]，当储层压裂后，这些结构被重新“激活”扩张，这类天然裂缝对页岩气藏渗流能力的提升不能像常规、致密气藏那样忽略^[9-11]。(2)基于现有“基质孔隙+双重裂缝”形式的三孔隙页岩基质模型建立的渗流方程非线性程度高^[12-13]，难以用其进行渗流机理的拓展研究。(3)当前模型大多应用Fick第一定律描述解吸气的扩散过程^[14-15]，而非稳态扩散更符合实际情形。(4)气体在纳米级基质孔隙中低速流动时需克服一定程度的启动压力梯度，但现有模型对该效应缺乏考虑。

针对以上不足，本文提出一种新的圆柱状三孔隙页岩基质模型，与五线性流模型结合，建立考虑基质启动压力梯度的页岩气藏多级压裂水平井不稳定渗流模型；借助Green函数、Laplace变换、Stehfest数值反演得到了模型的解，并绘制无因次拟压力及其导数曲线；根据曲线和渗流特征，划分模型的各渗流阶段，并进行相关参数的敏感性分析。

1 页岩气藏压裂水平井物理模型 1.1 页岩三重孔隙基质模型

Apaydin等人曾提出基于De Swaan模型^[16]的改进三重孔隙模型^[17]，如图 2所示。该模型将球形颗粒划分为均质基质内核与外部一层基质内嵌裂缝的双重孔隙同心球壳，各基质颗粒之间的区域为微裂缝系统，该方法得到广泛采用^[18-20]。

然而，当考虑启动压力梯度时，基于Apaydin模型所建立的球坐标偏微分方程难以求解。为解决该问题，本文借鉴Apaydin模型，建立一种柱状三重孔隙页岩基质模型，模型由柱状均质基质内核与薄筒状的双重介质表层构成，各圆柱基质间的空间为微裂缝部分(图 3)，基质内核中的气体沿径向流至表层外壳的天然裂缝。显然，基于本文模型所建立的页岩基质渗流方程为柱坐标偏微分方程，只需借助Green函数即可求解由于启动压力梯度在渗流方程中产生的非齐次项。

1.2 页岩气藏压裂水平井五线性流模型

目前，页岩气藏压裂水平井五线性流模型应用较为广泛^[21-23]，该模型的各流动区域位置关系及流体流动方向如图 4所示。

将本文提出的柱状页岩三重孔隙基质模型与五线性流模型结合如图 5所示，各流动区具有以下特点与关系：(1)区域一为压裂改造区，由圆柱形三重孔隙基质与微裂缝系统构成，内部流体由微裂缝流向水力裂缝区；(2)区域二为均质单一孔隙基质，流体流向区域一的微裂缝系统，流线方向垂直于水力裂缝延伸方向；(3)区域三和区域四视为均质的单一孔隙介质，其内部流体分别流向区域一和区域二。

2 页岩基质系统渗流模型

为了便于研究，对模型作如下基本假设：(1)页岩气藏为矩形，厚度为$h$，外边界封闭。(2)仅存在单相气体渗流，且不考虑温度、重力对渗流的影响。(3)基质孔隙中气体渗流考虑Knudsen扩散以及启动压力梯度的影响；裂缝系统的渗流符合Darcy定律；吸附气吸附解吸满足Langmuir等温吸附定律。(4)模型中各介质间窜流均为非稳态。

2.1 页岩基质内核的渗流控制方程

考虑基质孔隙内气体输运过程中的滑脱、Knudsen扩散效应，其运动方程表示为

$ {v_{{\rm{mc}}}} = {v_{\rm{D}}} + {v_{{\rm{Kn}}}} = - \left( {F\dfrac{{{K_{\rm{m}}}}}{{{\mu _{\rm{g}}}}} + {D_{\rm{K}}}{c_{\rm{g}}}} \right)\dfrac{{\partial {p_{{\rm{mc}}}}}}{{{r_{{\rm{mc}}}}}} $

(1)

式中：

$v_{{\rm{mc}}}$—基质中气体渗流速度，m/s；

$v_{{\rm{D}}}$—受黏性流控制的速度分量，m/s；

$v_{{\rm{Kn}}}$—受Knudsen扩散效应控制的速度分量，m/s；

$F$—滑脱因子，无因次；

$K_{{\rm{m}}}$—基质渗透率，mD；

$\mu_{{\rm{g}}}$—气体黏度，Pa$\cdot$s；

$D_{{\rm{K}}}$—Knudsen扩散系数，m$^2$/s；

$c_{{\rm{g}}}$—气体压缩系数，Pa$^{-1}$；

$p_{{\rm{mc}}}$—基质内核孔隙压力，Pa；

$r_{{\rm{mc}}}$—基质核中任意点至圆柱基质轴线的距离，m。

根据Javadpour提出的滑脱因子计算方法^[24]以及Igwe提出的Knudsen扩散系数方法^[25]，并引入表观渗透率的概念，运动方程可进一步表示为

$ {v_{{\rm{mc}}}} = - \dfrac{{{K_{\rm{m}}}}}{{{\mu _{\rm{g}}}}}\left\{ {\left[ {1 + \sqrt {\dfrac{{8\mathsf{π}{\rm R}T}}{{{M_{\rm{g}}}}}} \dfrac{{{\mu _{\rm{g}}}}}{{{r_{\rm{p}}}{p_{{\rm{avg}}}}}}\left( {\dfrac{2}{\alpha } - 1} \right)} \right] + \dfrac{{2{r_{\rm{p}}}}}{3}\sqrt {\dfrac{{8{\rm R}T}}{{\mathsf{π}{M_{\rm{g}}}}}} \dfrac{{{c_{\rm{g}}}{\mu _{\rm{g}}}}}{{{K_{\rm{m}}}}}} \right\}\\\dfrac{{\partial {p_{{\rm{mc}}}}}}{r_{\rm{mc}}} = - \dfrac{{{K_{{\rm{m, a}}}}}}{{{\mu _{\rm{g}}}}}\dfrac{{\partial {p_{{\rm{mc}}}}}}{{{r_{{\rm{mc}}}}}} $

(2)

式中：

R—气体常数，R=8.314 J/(mol$\cdot$K)；

$T$—体系温度，K；

$M_{{\rm{g}}}$—气体摩尔质量，kg/mol；

$r_{{\rm{p}}}$—孔隙半径，m；

$p_{{\rm{avg}}}$—体系平均压力，Pa；

$\alpha$—角动量调节因子，无因次；

$K_{{\rm{m, a}}}$—基质表观渗透率，mD。

基于Prada和Civan提出的启动压力梯度表示方法^[26]结合式(2)，可得连续性方程

$ \nabla \left[ {{\rho _{\rm{g}}}\dfrac{{{K_{{\rm{mc, a}}}}}}{{{\mu _{\rm{g}}}}}\left( {\nabla {p_{{\rm{mc}}}} - {\lambda _{\rm{T}}}} \right)} \right] + {\widetilde q_{{\rm{mc}}}} = \dfrac{\partial }{{\partial t}}\left( {{\phi _{{\rm{mc}}}}{\rho _{\rm{g}}}} \right) $

(3)

式中：$\rho_{{\rm{g}}}$—气体密度，kg/m$^3$；

$\lambda_{{\rm{T}}}$—启动压力梯度，Pa/m；

$\widetilde {q}_{{\rm{mc}}}$—基质孔隙表面吸附气体解吸的质量流量强度，kg/(m$^3$$\cdot$s)；

$\phi_{\rm{mc}}$—基质内核孔隙度，无因次；

$K_{\rm{mc, a}}$—基质内核表观渗透率，mD。

基于式(3)，结合Langmuir等温吸附原理并引入拟时间概念^[27]，得到圆柱基质内核渗流方程

$ \dfrac{1}{{{r_{{\rm{mc}}}}}}\dfrac{\partial }{{\partial {r_{{\rm{mc}}}}}}\left[ {{r_{{\rm{mc}}}}\left( {\dfrac{{\partial {m_{{\rm{mc}}}}}}{{\partial {r_{{\rm{mc}}}}}} - {\lambda _{{\rm{Tm}}}}} \right)} \right] = \\ {\kern 40pt} \dfrac{{{\phi _{{\rm{mc}}}}{\mu _{{\rm{gi}}}}{{\overline c }_{{\rm{tmci}}}}}}{{{K_{{\rm{mc, a}}}}}}\dfrac{{\partial {m_{{\rm{mc}}}}}}{{\partial {t_{\rm{a}}}}} $

(4)

其中

$ m = 2\int\limits_{\rm{0}}^p {\dfrac{p}{{\mu z}}{\rm{d}}p}, \\{\overline c _{{\rm{tmc}}}} = \dfrac{{2{p_{{\rm{sc}}}}T}}{{{T_{{\rm{sc}}}}}}\dfrac{{{V_{\rm{L}}}{p_{\rm{L}}}z}}{{{{({p_{\rm{L}}} + {p_{\rm{m}}})}^2}{p_{\rm{m}}}}}\dfrac{{(1 - {\phi _{{\rm{mc}}}})}}{{{\phi _{\rm{m}}}}}{\rm{ + }}{c_{{\rm{tmc}}}}, \\{\lambda _{{\rm{Tm}}}}{\rm{ = }}\dfrac{{2{p_{{\rm{mc}}}}}}{{{\mu _{\rm{g}}}z}}{\lambda _{\rm{T}}}, \\{t_{\rm{a}}} = \int\limits_{{t_0}}^t {\dfrac{{{\mu _{{\rm{gi}}}}{c_{{\rm{gi}}}}}}{{\mu (t){c_{\rm{g}}}(t)}}} {\rm{d}}t $

式中：

$m_{{\rm{mc}}}$—基质内核孔隙拟压力，Pa/s；

$\lambda_{{\rm{Tm}}}$—拟启动压力梯度，Pa/(m$\cdot$s)；

$\mu_{{\rm{gi}}}$—初始压力下的气体黏度，Pa$\cdot$s；

${c}_{{\rm{tmci}}}$—原始压力下的基质总压缩系数，Pa$^{-1}$；

$t_{{\rm{a}}}$—拟时间，s；

$z$—天然气偏差因子，无因次；

$V_{{\rm{L}}}$—Langmuir吸附体积，m$^3$/m$^3$；

$p_{{\rm{L}}}$—Langmuir压力，Pa；

${c}_{{\rm{tmc}}}$—基质总压缩系数，Pa$^{-1}$；

$c_{{\rm{gi}}}$—初始压力下的气体压缩系数，Pa$^{-1}$。

2.2 页岩基质表层的渗流控制方程

对于基质表层，可将圆柱“薄筒”近似展开成多个以栅状平行排列且尺寸相同的平板状天然裂缝和基质，如图 6所示。此时页岩气由天然裂缝向微裂缝的流动、表层基质向天然裂缝的流动简化为一维直线流。

由此可得页岩表层基质、天然裂缝的渗流方程

$ \dfrac{{{\partial ^2}{m_{{\rm{ms}}}}}}{{\partial {\chi ^2}}} = \dfrac{{{\phi _{{\rm{ms}}}}{\mu _{{\rm{gi}}}}{{\overline c }_{{\rm{tmsi}}}}}}{{{K_{{\rm{ms, a}}}}}}\dfrac{{\partial {m_{{\rm{ms}}}}}}{{\partial {t_{\rm{a}}}}} $

(5)

$ \dfrac{{{\partial ^2}{m_{{\rm{nf}}}}}}{{\partial {r^2}}} + \dfrac{2}{{{K_{{\rm{nf}}}}}}\dfrac{{{\rm R}T}}{{{M_{\rm{g}}}}}{\left. {{{\widetilde q}_{{\rm{msf}}}}} \right|_{{h_{{\rm{mm}}}}/2}} = \dfrac{{{\phi _{{\rm{nf}}}}{\mu _{\rm{gi}}}{{\overline c }_{{\rm{tnfi}}}}}}{{{K_{{\rm{nf}}}}}}\dfrac{{\partial {m_{{\rm{nf}}}}}}{{\partial {t_{\rm{a}}}}} $

(6)

式中：

$m_{\rm{ms}}$—表层基质孔隙拟压力，Pa/s；

$\chi$—表层基质中任意点至中心面的距离，m；

$\phi_{\rm{ms}}$—表层基质孔隙度，无因次；

${c}_{{\rm{tmsi}}}$—初始地层压力下的表层基质总压缩系数，Pa$^{-1}$；

$K_{\rm{ms, a}}$—表层基质表观渗透率，mD；

$m_{\rm{nf}}$—天然裂缝孔隙拟压力，Pa/s；

$K_{{\rm{nf}}}$—天然裂缝渗透率，mD；

$\widetilde{q}_{{\rm{msf}}}$—流体由表层基质向天然裂缝的质量流量强度kg/(m$^3$$\cdot$s)；

$h_{\rm{mm}}$—表层双重介质展开后的平板状基质厚度，m；

$\phi_{{\rm{nf}}}$—天然裂缝孔隙度，无因次；

${c}_{{\rm{tnfi}}}$—初始地层压力下的天然裂缝总压缩系数，Pa$^{-1}$。

式(6)中，表层基质流体流向天然裂缝的质量流量强度表达式如下

$ {\widetilde q_{{\rm{msf}}}}({R_{\rm{m}}}, t) = \dfrac{{{q_{{\rm{ms}}}}({h_{{\rm{mm}}}}/2, t)}}{{{V_{{\rm{nf}}}}/2}}=\\ {\kern 40pt} 2{\rho _{\rm{g}}}\dfrac{{N{A_{{\rm{ms}}}}}}{{N{A_{{\rm{ms}}}}{h_{{\rm{nf}}}}}}{V_{{\rm{ms}}}}({{{h_{{\rm{mm}}}}}}/{2}, t) = \\ {\kern 80pt} - \dfrac{{{M_{\rm{g}}}}}{{{\rm{R}}T}}\dfrac{{{K_{{\rm{ms, a}}}}}}{{{h_{{\rm{nf}}}}}}{\left. {\dfrac{{\partial {m_{{\rm{ms}}}}}}{{\partial \chi }}} \right|_{\chi = {h_{{\rm{mm}}}}/2}} $

(7)

式中：$ q_{\rm{ms}}$—表层基质窜流至天然裂缝的质量流量，kg/s；

$V_{{\rm{nf}}}$—天然裂缝体积，m$^3$；

$N$—表层介质展开后的天然裂缝条数，常数；

$A_{{\rm{ms}}}$—基质与天然裂缝接触面积，m$^2$；

$h_{{\rm{nf}}}$—天然裂缝宽度，m；

$V_{{\rm{ms}}}$—渗流速度，m/s。

3 页岩气藏压裂水平井渗流模型求解 3.1 参数处理

为简化渗流方程，定义导压系数如下

$ {\eta _\xi } = \dfrac{{{K_\xi }}}{{{\mu _{{\rm{gi}}}}{{\left( {\phi {{\overline c }_{\rm{t}}}} \right)}_\xi }}} $

(8)

式中：

$\eta_\xi$—导压系数，m$^2$/s；

下标$\xi$取mc、ms、2m、3m、4m、nf、f、F，其中，mc—基质核心；ms—表层基质；2m—二区基质；3m—三区基质；4m—四区基质；nf—天然裂缝；f—微裂缝；F—水力裂缝。

同时，按照表 1定义的无因次变量，对渗流方程进行无因次化处理。

3.2 模型求解

首先对各区渗流方程进行Laplace变换，完成各区域渗流方程的求解，具体过程如下。

区域四的渗流方程为

$ \dfrac{{{\partial ^2}{{\overline m }_{{\rm{4D}}}}}}{{\partial x_{\rm{D}}^{\rm{2}}}} - \dfrac{s}{{{\eta _{{\rm{4mD}}}}}}{\overline m _{{\rm{4D}}}} = 0 $

(9)

式中：

$s$—Laplace变量。

满足边界条件

$ {\left. {\dfrac{{\partial {{\overline m }_{{\rm{4D}}}}}}{{\partial {x_{\rm{D}}}}}} \right|_{{x_{\rm{D}}} = {x_{{\rm{eD}}}}}}{\rm{ = }}0, \\ {\left. {{{\overline m }_{{\rm{4D}}}}} \right|_{{x_{\rm{D}}} = 1}} = {\left. {{{\overline m }_{{\rm{2D}}}}} \right|_{{x_{\rm{D}}} = 1}} $

由此得到区域四的渗流方程解

$ {\overline m _{{\rm{4D}}}} = {\left. {{{\overline m }_{{\rm{2D}}}}} \right|_{{x_{\rm{D}}} = 1}}\dfrac{{\cosh \left[({x_{\rm{D}}} - {x_{{\rm{eD}}}})\sqrt {{u_4}} \right]}}{{\cosh \left[(1 - {x_{{\rm{eD}}}})\sqrt {{u_4}}\right]}} $

(10)

式中：

${u_4} = s/{\eta _{{\rm{4mD}}}}$。

区域三与区域四的渗流方程形式上基本一致、解法相同。区域三渗流方程解为

$ {\overline m _{{\rm{3D}}}} = {\left. {{{\overline m }_{{\rm{fD}}}}} \right|_{{x_{\rm{D}}} = 1}}\dfrac{{\cosh \left[({x_{\rm{D}}} - {x_{{\rm{eD}}}})\sqrt {{u_3}} \right]}}{{\cosh \left[(1 - {x_{{\rm{eD}}}})\sqrt {{u_3}} \right]}} $

(11)

式中：

${u_3} = s/{\eta _{{\rm{3mD}}}}$。

区域二渗流方程为

$ \dfrac{{{\partial ^2}{{\overline m }_{{\rm{2D}}}}}}{{\partial y_{\rm{D}}^2}}{\rm{ + }}\dfrac{{{K_{{\rm{m4}}}}}}{{{K_{{\rm{m2}}}}}}{\left. {\dfrac{{\partial {{\overline m }_{{\rm{4D}}}}}}{{\partial {x_{\rm{D}}}}}} \right|_{{x_{\rm{D}}} = 1}} = \dfrac{s}{{{\eta _{{\rm{2mD}}}}}}{\overline m _{{\rm{2D}}}} $

(12)

该区域内边界条件

$ {\left. {\dfrac{{\partial {{\overline m }_{{\rm{2D}}}}}}{{\partial {y_{\rm{D}}}}}} \right|_{{y_{\rm{D}}} = {y_{{\rm{eD}}}}}} = 0, $

外边界条件

$ {\left. {{m_{{\rm{2D}}}}} \right|_{{y_{\rm{D}}} = {y_{{\rm{1D}}}}}}{\rm{ = }}{\left. {{m_{{\rm{fD}}}}} \right|_{{y_{\rm{D}}} = {y_{{\rm{1D}}}}}} $

由此，得到式(12)解为

$ {\overline m _{{\rm{2D}}}} = {\left. {{{\overline m }_{{\rm{fD}}}}} \right|_{{y_{\rm{D}}} = {y_{1{\rm{D}}}}}}\dfrac{{\cosh \left[({y_{\rm{D}}} - {y_{{\rm{eD}}}})\sqrt {{u_2}(s)} \right]}}{{\cosh \left[({y_{{\rm{1D}}}} - {y_{{\rm{eD}}}})\sqrt {{u_2}(s)} \right]}} $

(13)

式中：

${u_2} = \dfrac{s}{{{\eta _{{\rm{2mD}}}}}} - \dfrac{{{K_{{\rm{m4}}}}}}{{{K_{{\rm{m2}}}}}}\sqrt {{u_4}(s)} \tanh [(1 - {x_{{\rm{eD}}}})\sqrt {{u_4}(s)} ]$。

区域一基质内核渗流方程为

$ \dfrac{{{\partial ^2}{{\overline m }_{{\rm{mcD}}}}}}{{\partial r_{{\rm{mcD}}}^2}} + \dfrac{1}{{{r_{{\rm{mcD}}}}}}\dfrac{{\partial {{\overline m }_{{\rm{mcD}}}}}}{{\partial {r_{{\rm{mcD}}}}}} + \dfrac{{{\lambda _{{\rm{TmD}}}}}}{{s{r_{{\rm{mcD}}}}}} -\\ {\kern 40pt} \dfrac{s}{{{\eta _{{\rm{mcD}}}}}}{\overline m _{{\rm{mcD}}}} = 0 $

(14)

根据基质内核渗流方程满足的边界条件

$ \mathop {\lim }\limits_{{r_{\rm{D}}} \to 0} {\dfrac{{\partial {{\overline m }_{{\rm{mcD}}}}}}{{\partial {r_{\rm{D}}}}}} = 0\\ {\left. {{{\overline m }_{{\rm{mcD}}}}} \right|_{{r_{\rm{D}}} = {r_{{\rm{mcD}}}}}} = {\left. {{{\overline m }_{{\rm{nfD}}}}} \right|_{{r_{\rm{D}}} = {r_{{\rm{mcD}}}}}} $

得到方程解为

$ {\overline m _{{\rm{mcD}}}} = \left[ {{{\left. {{{\overline m }_{{\rm{mfD}}}}} \right|}_{{r_{\rm{D}}} = {R_{{\rm{mcD}}}}}} - \int\limits_0^{ + \infty } {{\mathop{\rm G}\nolimits} ({r_{{\rm{mcD}}}}, \tau )} {\rm{d}}\tau } \right]\cdot \\ {\kern 40pt} \dfrac{{{I_0}\left( {{r_{{\rm{mcD}}}}\sqrt {{u_{{\rm{mc}}}}} } \right)}} {{{I_0}\left( {{R_{{\rm{mcD}}}}\sqrt {{u_{{\rm{mc}}}}} } \right)}} + \int\limits_0^{ + \infty } {{\mathop{\rm G}\nolimits} ({r_{{\rm{mcD}}}}, \tau )} {\rm{d}}\tau $

(15)

式中：

${u_{{\rm{mc}}}} = \dfrac{s}{{{\eta _{{\rm{mcD}}}}}}$；

G($r_{\rm{mcD}}$, $\tau$)为Green函数，表达式为

$ {\mathop{\rm G}\nolimits} ({r_{{\rm{mcD}}}}, s) =\\ {\kern 40pt} \dfrac{{{\lambda _{{\rm{TmD}}}}}}{s}\left\{ \begin{array}{l} {I_0}(\tau \sqrt {{u_{{\rm{mc}}}}} ){K_0}({r_{{\rm{mcD}}}}\sqrt {{u_{{\rm{mc}}}}} ) , \\ (0 < \tau \leqslant {r_{{\rm{mcD}}}})\\ {I_0}({r_{{\rm{mcD}}}}\sqrt {{u_{{\rm{mc}}}}} ){K_0}(\tau \sqrt {{u_{{\rm{mc}}}}} ) , \\ ({r_{{\rm{mcD}}}} < \tau \leqslant + \infty ) \end{array} \right. $

(16)

表层基质渗流数学方程为

$ \dfrac{{{\partial ^2}{{\overline m }_{{\rm{msD}}}}}}{{\partial \chi _{\rm{D}}^2}} - \dfrac{s}{{{\eta _{{\rm{msD}}}}}}{\overline m _{{\rm{msD}}}} = 0 $

(17)

式(17)满足边界条件

$ {\left. {{{\overline m }_{{\rm{msD}}}}} \right|_{\chi = 1}} = {\left. {{{\overline m }_{{\rm{mfD}}}}} \right|_{{\chi _{\rm{D}}} = 1}}, \\{\left. {\dfrac{{\partial {{\overline m }_{{\rm{msD}}}}}}{{\partial {\chi _{\rm{D}}}}}} \right|_{\chi = 0}} = 0 $

由此解得

$ {\overline m _{{\rm{msD}}}} = {\left. {{{\overline m }_{{\rm{nfD}}}}} \right|_{{\chi _{\rm{D}}} = 1}}\dfrac{{\cosh \left( {{\chi _{\rm{D}}}\sqrt {{u_{{\rm{ms}}}}} } \right)}}{{\cosh {\sqrt {{u_{{\rm{ms}}}}} } }} $

(18)

式中：

${u_{{\rm{ms}}}} = s/{\eta _{{\rm{msD}}}}$。

天然裂缝渗流数学方程为

$ \dfrac{{{\partial ^2}{{\overline m }_{{\rm{nfD}}}}}}{{\partial r_{\rm{D}}^2}} - \dfrac{1}{{{h_{{\rm{mmD}}}}{h_{{\rm{nfD}}}}}}\dfrac{{{K_{\rm{ms, a}}}}}{{{K_{{\rm{nf}}}}}}{\left. {\dfrac{{\partial {{\overline m }_{{\rm{msD}}}}}}{{\partial {\chi _{\rm{D}}}}}} \right|_{{\chi _{\rm{D}}} = 1}} -\\ {\kern 40pt} \dfrac{s}{{{\eta _{{\rm{msD}}}}}}{\overline m _{{\rm{msD}}}} = 0 $

(19)

根据边界条件

$ {\left. {\dfrac{{\partial {{\overline m }_{{\rm{nfD}}}}}}{{\partial {r_{\rm{D}}}}}} \right|_{{r_{\rm{D}}} = {R_{{\rm{mcD}}}}}} = \dfrac{{{k_{{\rm{mc, a}}}}}}{{{k_{{\rm{nf}}}}}}\dfrac{{{h_{{\rm{mm}}}}}}{{{h_{{\rm{nf}}}}}}{\left. {\dfrac{{\partial {{\overline m }_{{\rm{mcD}}}}}}{{\partial {r_{{\rm{mcD}}}}}}} \right|_{{r_{\rm{D}}} = {R_{{\rm{mcD}}}}}}, \\ {\left. {{{\overline m }_{{\rm{nfD}}}}} \right|_{{r_{\rm{D}}} = {R_{{\rm{msD}}}}}} = {\left. {{{\overline m }_{{\rm{fD}}}}} \right|_{{r_{\rm{D}}} = {R_{{\rm{msD}}}}}} $

天然裂缝渗流方程解为

$ {\overline m _{{\rm{nfD}}}} = {\overline m _{{\rm{fD}}}}\left\{ \begin{array}{l} \cosh \left[ {\left( {{r_{\rm{D}}} - {R_{{\rm{mcD}}}}} \right)\sqrt {{f_{{\rm{mf}}}}} } \right] + \dfrac{{{\theta _1}}}{{{\theta _3}}}\dfrac{{\sinh \left[ {\left( {{r_{\rm{D}}} - {R_{{\rm{mcD}}}}} \right)\sqrt {{f_{{\rm{mf}}}}} } \right]}}{{\sqrt {{f_{{\rm{mf}}}}} {{\cosh }^2}\left[\left( {{R_{{\rm{msD}}}} - {R_{{\rm{mcD}}}}} \right)\sqrt {{f_{{\rm{mf}}}}} \right]}} - \\ \tanh \left[\left( {{R_{{\rm{msD}}}} - {R_{{\rm{mcD}}}}} \right)\sqrt {{u_{{\rm{mf}}}}} \right]\sinh \left[ {\left( {{r_{\rm{D}}} - {R_{{\rm{mcD}}}}} \right)\sqrt {{f_{{\rm{mf}}}}} } \right] \end{array} \right\} - \\ {\kern 40pt} \left( {\dfrac{{{\theta _1}{\theta _4}}}{{{\theta _3}}} + {\theta _2}} \right)\dfrac{{\sinh \left[ {\left( {{r_{\rm{D}}} - {R_{{\rm{mcD}}}}} \right)\sqrt {{f_{{\rm{mf}}}}} } \right]}}{{\sqrt {{f_{{\rm{mf}}}}} \cosh \left[\left( {{R_{{\rm{msD}}}} - {R_{{\rm{mcD}}}}} \right)\sqrt {{f_{{\rm{mf}}}}} \right]}} $

(20)

式中：

$ {\theta _1} = \dfrac{{{K_{{\rm{mc, a}}}}{h_{{\rm{mm}}}}}}{{{K_{{\rm{mf}}}}{h_{{\rm{nf}}}}}}\dfrac{{{I_1}({R_{{\rm{mcD}}}}\sqrt {{u_{{\rm{mc}}}}} )}}{{{I_0}({R_{{\rm{mcD}}}}\sqrt {{u_{{\rm{mc}}}}} )}}\sqrt {{u_{{\rm{mc}}}}} $

$ {\theta _2} = {\theta _1}\int\limits_0^{ + \infty } {G({R_{{\rm{mcD}}}}, \tau ){\rm{d}}\tau } - \\ {\kern 40pt} {\left. {\dfrac{{\rm{d}}}{{{\rm{d}}{r_{{\rm{mcD}}}}}}\left( {\int\limits_0^{ + \infty } {G({r_{{\rm{mcD}}}}, \tau ){\rm{d}}\tau } } \right)} \right|_{{r_{{\rm{mcD}}}} = {R_{{\rm{mcD}}}}}} $

$ {\theta _3} = 1 - \dfrac{{{\theta _1}}}{{\sqrt {{f_{{\rm{mf}}}}} }}\tanh \left[ {\left( {{R_{{\rm{msD}}}} - {R_{{\rm{mcD}}}}} \right)\sqrt {{f_{{\rm{mf}}}}} } \right] \\ {\theta _4} = \dfrac{{{\theta _2}}}{{\sqrt {{f_{{\rm{mf}}}}} }}\tanh \left[ {\left( {{R_{{\rm{msD}}}} - {R_{{\rm{mcD}}}}} \right)\sqrt {{f_{{\rm{mf}}}}} } \right]\\ {f_{{\rm{nf}}}} = \dfrac{1}{{{h_{{\rm{mmD}}}}{h_{{\rm{nfD}}}}}}\dfrac{{{K_{{\rm{mc, a}}}}}}{{{K_{{\rm{nf}}}}}}\sqrt {\dfrac{s}{{{\eta _{{\rm{msD}}}}}}} \tanh {\sqrt {\dfrac{s}{{{\eta _{{\rm{msD}}}}}}} } + \dfrac{s}{{{\eta _{{\rm{nfD}}}}}} $

微裂缝渗流方程为

$ \dfrac{{{\partial ^{\rm{2}}}{{\overline m }_{{\rm{fD}}}}}}{{\partial y_{\rm{D}}^2}} - \dfrac{{2{K_{\rm{f}}}{h_{\rm{nfD}}}}}{{{K_{{\rm{nf}}}}{h_{\rm{fD}}}{h_{{\rm{mmD}}}}}}{\left. {\dfrac{{\partial {{\overline m }_{{\rm{nfD}}}}}}{{\partial {r_{\rm{D}}}}}} \right|_{{r_{\rm{D}}} = {R_{{\rm{msD}}}}}} + \\ {\kern 40pt} \dfrac{{{K_{{\rm{m3}}}}}}{{{K_{{\rm{fi}}}}}}{\left. {\dfrac{{\partial {{\overline m }_{{\rm{m3}}}}}}{{\partial {x_{\rm{D}}}}}} \right|_{{x_{\rm{D}}} = 1}} - \dfrac{s}{{{\eta _{{\rm{fD}}}}}}{\overline m _{{\rm{fD}}}} = 0 $

(21)

满足边界条件

$ {\left. {\dfrac{{\partial {{\overline m }_{{\rm{fD}}}}}}{{\partial {y_{\rm{D}}}}}} \right|_{{y_{\rm{D}}} = {y_{1{\rm{D}}}}}} = \dfrac{{{K_{{\rm{m2}}}}}}{{{K_{\rm{f}}}}}{\left. {\dfrac{{\partial {{\overline m }_{{\rm{2D}}}}}}{{\partial {y_{{\rm{1D}}}}}}} \right|_{{y_{\rm{D}}} = {y_{{\rm{1D}}}}}}, \\ {\left. {{{\overline m }_{{\rm{fD}}}}} \right|_{{w_{{\rm{FD}}}}/2}} = {\left. {{{\overline m }_{{\rm{FD}}}}} \right|_{{w_{{\rm{FD}}}}/2}} $

得到微裂缝渗流方程解

$ {\overline m _{{\rm{fD}}}} = {\overline m _{{\rm{FD}}}}{\zeta _2} + {\zeta _3} $

(22)

式中：

$ {\zeta _3} = \dfrac{{{a_{\rm{f}}}}}{{\sqrt {{f_{\rm{f}}}} }} \dfrac{1}{{1 + {\zeta _1}\tanh \left[ {\left( {{w_{{\rm{FD}}}}/2 - {y_{{\rm{1D}}}}} \right)\sqrt {{f_{\rm{f}}}} } \right]}}\dfrac{{\rm{1}}}{{\cosh \left[\left( {{w_{{\rm{FD}}}}/2 - {y_{{\rm{1D}}}}} \right)\sqrt {{f_{\rm{f}}}} \right]}} - \\ {\kern 40pt} \dfrac{{{a_{\rm{f}}}}}{{\sqrt {{f_{\rm{f}}}} }} \tanh \left[ {\left( {{w_{{\rm{FD}}}}/2 - {y_{{\rm{1D}}}}} \right)\sqrt {{f_{\rm{f}}}} } \right]\left\{ {1 - \dfrac{{{\zeta _1}}}{{1 + {\zeta _1}\tanh \left[ {\left( {{w_{{\rm{FD}}}}/2 - {y_{{\rm{1D}}}}} \right)\sqrt {{f_{\rm{f}}}} } \right]}}} \right\} $

$ {f_{\rm{f}}} = - \dfrac{{2{K_{\rm{f}}}{h_{{\rm{nfD}}}}}}{{{K_{{\rm{nf}}}}{h_{{\rm{fD}}}}{h_{{\rm{mmD}}}}}}\left\{ {\sqrt {{u_{{\rm{nf}}}}} \tanh \left[ {\left( {{R_{{\rm{msD}}}} - {R_{{\rm{mcD}}}}} \right)\sqrt {{u_{{\rm{nf}}}}} } \right] + \dfrac{{{\theta _1}{\theta _4}}}{{{\theta _3}}}} \right\}\\ + \dfrac{s}{{{\eta _{{\rm{fD}}}}}} - \dfrac{{{K_{{\rm{3m}}}}}}{{{K_{\rm{f}}}}}\sqrt {{u_3}} \tanh \left[ {\left( {1 - {x_{{\rm{eD}}}}} \right)\sqrt {{u_3}} } \right] $

$ {\zeta _1} = \dfrac{{{K_{{\rm{m2}}}}}}{{{K_{\rm{f}}}}}\sqrt {\dfrac{{{u_2}}}{{{f_{\rm{f}}}}}} \tanh \left[ {\left( {{y_{{\rm{1D}}}} - {y_{{\rm{eD}}}}} \right)\sqrt {{u_2}} } \right] $

$ {\zeta _2} = \sqrt {{f_{\rm{f}}}} \dfrac{{\tanh \left[ {\left( {{w_{{\rm{FD}}}}/2 - {y_{{\rm{1D}}}}} \right)\sqrt {{f_{\rm{f}}}} } \right]{\rm{ + }}{\zeta _1}}}{{1 + {\zeta _1}\tanh \left[ {\left( {{w_{{\rm{FD}}}}/2 - {y_{{\rm{1D}}}}} \right)\sqrt {{f_{\rm{f}}}} } \right]}} $

$ {a_{\rm{f}}} = \dfrac{{2{K_{\rm{f}}}{h_{{\rm{nfD}}}}}}{{{K_{{\rm{nf}}}}{h_{{\rm{fD}}}}{h_{{\rm{mmD}}}}}}\left( {\dfrac{{{\theta _1}{\theta _4}}}{{{\theta _3}}} + {\theta _2}} \right) $

水力裂缝渗流方程为

$ \dfrac{{{\partial ^2}{{\overline m }_{{\rm{FD}}}}}}{{\partial x_{\rm{D}}^{\rm{2}}}} + \dfrac{2}{{{F_{{\rm{CD}}}}}}{\left. {\dfrac{{\partial {{\overline m }_{{\rm{fD}}}}}}{{\partial {y_{\rm{D}}}}}} \right|_{{y_{\rm{D}}} = \dfrac{{{w_{F{\rm{D}}}}}}{2}}} - \dfrac{s}{{{\eta _{{\rm{FD}}}}}}{\overline m _{{\rm{FD}}}} = 0 $

(23)

根据水力裂缝渗流方程满足的边界条件

$ {\left. {\dfrac{{\partial {{\overline m }_{{\rm{FD}}}}}}{{\partial {x_{\rm{D}}}}}} \right|_{{x_{\rm{D}}} = 1}} = 0, \\ {\left. {\dfrac{{\partial {{\overline m }_{{\rm{FD}}}}}}{{\partial {x_{\rm{D}}}}}} \right|_{{x_{\rm{D}}} = 1}} = - \dfrac{{\mathsf{π}}}{{s {F_{{\rm{CD}}}}}} $

解得井底无因次拟压力

$ {\overline m _{{\rm{FD}}}} = \dfrac{\mathsf{π} }{{s{F_{{\rm{CD}}}}\sqrt {{f_{\rm{F}}}} \tanh \left( {\sqrt {{f_{\rm{F}}}} } \right)}} + \dfrac{{{a_{\rm{F}}}}}{{{f_{\rm{F}}}}} $

(24)

式中：

$ {f_{\rm{F}}} = \sqrt {\dfrac{s}{{{\eta _{{\rm{FD}}}}}}} - \dfrac{{2{\zeta _2}}}{{{F_{{\rm{CD}}}}}}\\ {a_{\rm{F}}} = \dfrac{{2{\zeta _3}}}{{{F_{{\rm{CD}}}}}} $

4 结果分析

利用Stehfest数值反演^[28]，可以将模型在Laplace空间下的解反演到真实空间。模型计算所用的基本参数如表 2所示。通过改变参数取值，对考虑各种复杂渗流机理的页岩气藏多段压裂水平井压裂响应曲线动态特征进行影响因素分析。

图 7为模型计算的井底无因次拟压力响应曲线，可划分7个流动阶段：(1)阶段Ⅰ为纯井筒储集阶段，拟压力及导数曲线呈斜率为1的直线；(2)阶段Ⅱ为过渡段；(3)阶段Ⅲ为水力裂缝、微裂缝双线性流阶段，无因次拟压力导数曲线呈斜率为0.25的直线；(4)阶段Ⅳ为天然裂缝至微裂缝的流体补充阶段；(5)阶段Ⅴ为页岩基质内部孔隙至天然裂缝的流体补充阶段；(6)阶段Ⅵ为为二、三、四区过渡流动阶段；(7)阶段Ⅶ为边界控制流动阶段，无因次拟压力及其导数曲线呈斜率为1的直线。

启动压力梯度对井底无因次拟压力的影响如图 8所示。启动压力梯度主要影响阶段Ⅴ及以后的流动阶段，启动压力梯度越大，渗流阻力越大、气井定产生产时所需的压差越大、无因次拟压力及其导数值越大。

图 9反映了基质渗透率对无因次拟压力曲线的影响。基质渗透率越低，气体在基质系统中的流动阻力越大，阶段Ⅴ和阶段Ⅵ的无因次拟压力及其导数值越大、曲线位置越高且斜率越大，阶段Ⅵ经历的时间越长，边界控制流动阶段越晚出现。

页岩基质的双重介质表层厚度$h_{\rm s}$对无因次拟压力曲线的影响如图 10所示，该参数主要影响了阶段Ⅳ和阶段Ⅴ。当表层基质厚度越大时，水力裂缝、微裂缝双线性流阶段经历时间越短，阶段Ⅳ和阶段Ⅴ内的拟压力导数值越大，并且曲线的斜率越大。

图 11反映了区域一的长度$y_1$值对无因次拟压力曲线的影响。若$y_1$值越小，则水力裂缝、微裂缝双线性流阶段经历的时间越短，阶段Ⅳ和阶段Ⅴ的拟压力及导数值越大、曲线位置越高。

5 结论

（1）本文基于柱状三重孔隙页岩基质模型和五线性流模型，建立并求解了考虑启动压力梯度的页岩气藏压裂水平井渗流模型，绘制了典型试井曲线并分析了相关参数的敏感性。

（2）启动压力梯度主要影响基质内部孔隙至天然裂缝的流体补充阶段及其以后的流动阶段。启动压力梯度越大，渗流阻力越大，气井定产生产时所需的压差越大，无因次拟压力及其导数值越大。

（3）基质渗透率主要影响页岩气渗流的中后期。基质渗透率越小，基质孔隙流体向微裂缝补充阶段无因次拟压力导数曲线位置越高且斜率越大、边界流动阶段开始时间越晚。

（4）页岩基质表层双重介质厚度越大，则双线性流阶段经历时间越短、天然裂缝至微裂缝流体补充阶段和页岩基质内部孔隙至天然裂缝流体补充阶段的无因次拟压力导数值越大、曲线斜率也越大。