引言

裂缝性气藏储层多为致密储层，裂缝是该类储层中重要的气体储存空间与传输通道。裂缝系统形态与体积决定了裂缝性气藏气井的产能与稳产能力^[1]。对于该类型气藏的试井模型研究方法主要分为4类：第1类研究方法为使用连续双重介质描述裂缝性气藏非均质储层，结合渗流力学理论，建立裂缝性气藏试井模型，模拟气体在裂缝性油气藏内的流动，该类方法无法描述裂缝系统形态与试井压力动态曲线的关系^[2-4]；第2类研究方法是基于数值试井理论方法，利用离散裂缝网络(DFN)描述裂缝性油气藏的大尺度裂缝，模拟其试井时的压力变化规律，该类方法的缺点在于仅能模拟大尺度裂缝，受限于计算机处理能力与内存限制，无法描述小尺寸分支裂缝及其形成的网络系统^[5-7]；第3类研究方法是使用分形连续双重介质描述不同尺度下储层的非均质性，该类方法同样存在裂缝形态描述上的问题^[8-10]；第4类方法使用分形参数控制不同长度的裂缝，形成平行或相交网络，嵌入储层基质中，形成具有多尺度分形性质的天然裂缝性油气藏模型^[11-13]，该方法能够很好地表征裂缝网络的形态与分布，但该类裂缝网络以直线和相交的直线为主，并不能完全模拟趋向井底的径向流动^[14-15]。

本文利用树状分形网络对趋向井底的径向流动模拟的优越性，使用树状分形网络模拟气藏裂缝系统，并将树状分形网络嵌入到基质中，提出基于树状分形网络的裂缝性气藏试井模型，分析基于树状分形网络的裂缝性气藏压力动态特征的影响因素，为裂缝性气藏的动态描述及试井分析提供了一种新方法。

1 物理模型

树状分形网络由一系列分叉结构组合而成，考虑到网络的对称性，可以在计算时只分析网络的一个单元，如图 1所示。生成网络之前，需要生成多个不同级(总分叉级数$M$ )的同心圆环，各个圆环以原点$O$ 为圆心，使树状分形网络的各级分叉末端处于同一级圆环上，即可生成树状分形网络^[16]。树状分形网络由$N$ 条初始裂缝分形得到，初始裂缝的初始长度为$l_{{\rm 0}}$ 、直径为$d_{{\rm 0}}$ 。树状分形网络采用“Y”型结构^{[17, 18]}，树状分形的网络结构参数采用两个尺度因子，分别为长度比$\alpha$ 和直径比$\beta$ ，“Y”型结构分叉角度为$\theta$ ($\theta < \pi/2$ )。

本文研究的是裂缝性气藏圆形地层，井位于圆心$O$ 处，储层厚度$h$ ，井半径$r_{{\rm w}}$ 。使用树状分形网络模拟地层裂缝系统，并将树状分形网络嵌入到水平圆形地层内。模型基于如下假设条件：地层可被树状分形网络划分为$M$ 个环状区域，每一区裂缝系统性质不同，而基质系统及流体性质相同；裂缝渗透率远大于基质渗透率，流体只能由裂缝系统流向井筒，储层内气体流动为等温渗流；流体在各区的基质系统、裂缝系统内的流动满足线性渗流规律并忽略重力影响。

2 数学模型

气藏储层由多个不同级(总分叉级数$M$ )以原点$O$ 为圆心的同心圆环组成(图 1)，裂缝系统参数直接由树状分形网络生成。$k$ 级裂缝即为区域$k$ ($k = 0, 1, \cdots, M$ )内的裂缝，单条$k$ 级裂缝的长度

$ {{l}_{k}}={{\alpha }^{k}}{{l}_{0}} $

(1)

式中：$ l_{{k}}$ —$k$ 级裂缝长度，m；

$\alpha$ —裂缝长度比，无因次；

$k$ —裂缝分叉级数，无因次；

$l_{{\rm 0}}$ —裂缝初始长度，m。

单条$k$ 级裂缝的直径为

$ {{d}_{k}}={{\beta }^{k}}{{d}_{0}} $

(2)

式中：$ d_{{ k}}$ —$k$ 级裂缝直径，m；

$\beta$ —裂缝直径比，无因次；

$d_{{\rm 0}}$ —裂缝初始直径，m。

径向距离为井到各区域界面的距离，由式(1)，可得

$ {{r}_{k}}=\sum\limits_{i=0}^{k}{{{l}_{i}}\cos \theta }={{l}_{0}}\left[1+\dfrac{\alpha \left( 1-{{\alpha }^{k}} \right)\cos \theta }{1-\alpha } \right] $

(3)

式中：$ r_{{k}}$ —$k$ 级裂缝径向距离，m；

$\theta$ —裂缝分叉角度，(°)。

Xu等将单条$k$ 级裂缝渗透率表示为^[19]

$ {{K}_{k}}=10^{-12}\dfrac{d_{k}^{2}}{32}\dfrac{1}{{{T}_{k}}} $

(4)

式中：$K_{{k}}$ —$k$ 级裂缝渗透率，D；

$T_{{k}}$ —区域$k$ 的曲迂度，无因次。

区域$k$ 的曲迂度$T_{k}$ 在模型中表示为

$ {{T}_{k}}=\dfrac{{{l}_{k}}}{{{r}_{k}}-{{r}_{k-1}}}=\left\{ \begin{array}{*{35}{l}}1, {\kern 32pt} k=0 \\ \sec \theta, {\kern 18pt} k>0 \\\end{array} \right. $

(5)

联立式(2)~式(5)，得到区域$k$ 的裂缝系统渗透率

$ {{K}_{{\rm f}k}}=\left\{ \begin{array}{*{35}{l}}\dfrac{10^{-12}Nd_{0}^{2}}{32}, {\kern 52pt} k=0 \\[8pt] \dfrac{10^{-12}N{{n}^{k}}{{\beta }^{2k}}d_{0}^{2}\cos \theta }{32}, {\kern 10pt} k>0 \\\end{array} \right. $

(6)

式中：${{K}_{{\rm f}k}}$ —区域$k$ 的裂缝系统渗透率，D；

$N$ —初始裂缝数，条；

$n$ —分形网络分叉数，本文$n=2$ 。

区域$k$ 的总系统体积为

$ {{V}_{{\rm t}k}}=\left\{ \begin{array}{*{35}{l}} \pi hr_{0}^{2}, {\kern 48pt} k=0 \\ \pi h\left( r_{k}^{2}-r_{k-1}^{2} \right), {\kern 10pt} k>0 \\ \end{array} \right. $

(7)

式中：${{V}_{{\rm t}k}}$ —区域$k$ 的总系统体积，m³；

$h$ —储层厚度，m；

$r_{{\rm 0}}$ —初始裂缝径向距离，m。

区域$k$ 的裂缝系统孔隙体积为

$ {{V}_{{\rm ft}k}}=N{{n}^{k}}\dfrac{\pi {{l}_{k}}d_{k}^{2}}{4}=\dfrac{N\pi {{n}^{k}}{{\alpha }^{k}}{{\beta }^{2k}}{{l}_{0}}d_{0}^{2}}{4} $

(8)

式中：${{V}_{{\rm ft}k}}$ —区域$k$ 的裂缝系统孔隙体积，m³。

区域$k$ 的基质系统孔隙体积为

$ {{V}_{{\rm mt}k}}={{\phi }_{{\rm mp}}}\left( {{V}_{{\rm t}k}}-{{V}_{{\rm tf}k}} \right) $

(9)

式中：${{V}_{{\rm mt}k}}$ —区域$k$ 的基质系统孔隙体积，m³；

${{\phi }_{{\rm mp}}}$ —基质孔隙度，%。

联立式(7)、式(8)，区域$k$ 的裂缝系统孔隙度表示为

$ {{\phi }_{{\rm f}k}}=\dfrac{{{V}_{{\rm ft}k}}}{{{V}_{{\rm t}k}}}=\left\{ \begin{array}{*{35}{l}}\dfrac{Nd_{0}^{2}}{4h{{l}_{0}}}, {\kern 65pt} k=0 \\[6pt] \dfrac{N{{n}^{k}}{{\alpha }^{k}}{{\beta }^{2k}}{{l}_{0}}d_{0}^{2}}{4h\left( r_{k}^{2}-r_{k-1}^{2} \right)} k=0, {\kern 5pt} k>0 \\\end{array} \right. $

(10)

式中：${{\phi }_{{\rm f}k}}$ —区域$k$ 的裂缝系统孔隙度，%。

联立式(7)~式(10)，区域$k$ 的基质系统孔隙度表示为

$ {{\phi }_{{\rm m}k}}=\dfrac{{{V}_{{\rm mt}k}}}{{{V}_{{\rm t}k}}}={{\phi }_{{\rm mp}}}\left( 1-{{\phi }_{{\rm f}k}} \right) $

(11)

式中：${{\phi }_{{\rm m}k}}$ —区域$k$ 的基质系统孔隙度，%。

在双分叉结构下，裂缝性气藏各区域裂缝系统渗透率、裂缝体积相等的条件为$n=2$ ，$\alpha =1$ ，$\beta =0.707$ 。

当$\beta < 0.707$ 时，裂缝性气藏各区域裂缝系统渗透率随半径增大而减小；当$\beta>0.707$ 时，裂缝性气藏各区域裂缝系统渗透率随半径增大而增大。根据渗流力学理论^[20]，基于树状分形网络的裂缝性气藏的渗流数学模型为^[14]

$ \left\{ \begin{array}{*{35}{l}} {{K}_{{\rm f}k}}\left( \dfrac{{{\partial }^{2}}{{\psi }_{{\rm f}k}}}{\partial {{r}^{2}}}+\dfrac{1}{r}\dfrac{\partial {{\psi }_{{\rm f}k}}}{\partial r} \right)+a{{K}_{{\rm m}}}\left( {{\psi }_{{\rm m}k}}-{{\psi }_{{\rm f}k}} \right)= \dfrac{\mu {{\phi }_{{\rm f}k}}{{C}_{{\rm tf}}}}{3.6}\dfrac{\partial {{\psi }_{{\rm f}k}}}{\partial t}, {\kern 28pt} \\ {{r}_{k-1}}\leqslant r \leqslant {{r}_{k}} \\ -a{{K}_{{\rm m}}}\left( {{\psi }_{{\rm m}k}}-{{\psi }_{{\rm f}k}} \right)=\dfrac{\mu {{\phi }_{{\rm m}k}}{{C}_{{\rm tm}}}}{3.6}\dfrac{\partial {{\psi }_{{\rm m}k}}}{\partial t}, {\kern 112pt} \\ {{r}_{k-1}}\leqslant r\leqslant {{r}_{k}} \\ {{\psi }_{{\rm f}k}}(r, 0)={{\psi }_{{\rm m}k}}(r, 0)={{\psi }_{{\rm i}}}, {\kern 158pt} \\ k=0, 1, \cdots, M \\ {{\left. \left( r\dfrac{\partial {{\psi }_{{\rm f}k}}}{\partial r} \right) \right|}_{{{r}_{{\rm e}}}={{r}_{{\rm w}}, k =1}}}{\rm =}-\left( 1-C\dfrac{{\rm d}{{\psi }_{{\rm wf}}}}{{\rm d}t} \right), {\kern 120pt} \\ {{\psi }_{{\rm wf}}}={{\left. {{\psi }_{{\rm f}k}} \right|}_{k=1, r={{r}_{{\rm w}}}{{{\rm e}}^{-S}}}} \\ {{\left. {{\psi }_{{\rm f}k}} \right|}_{r={{r}_{k}}}}={{\left. {{\psi }_{{\rm f}\left( k+1 \right)}} \right|}_{r={{r}_{k}}}}, {\kern 175pt} \\ k=0, 1, \cdots, M-1 \\[6pt] {{\left. \dfrac{\partial {{\psi }_{{\rm f}k}}}{\partial r} \right|}_{r={{r}_{k}}}}={{\left. \dfrac{{{K}_{{\rm f}\left( k+1 \right)}}}{{{K}_{{\rm f}k}}}\dfrac{\partial {{\psi }_{{\rm f}\left( k+1 \right)}}}{\partial r} \right|}_{r={{r}_{k}}}}, {\kern 131pt} \\ k=0, 1, \cdots, M-1 \\ {{\psi }_{{\rm f}k}}\left( t \right)={{\psi }_{{\rm i}}}, {\kern 217pt} \\ r\to \infty, k= M \\ \end{array} \right. $

(12)

式中：${{\psi }_{{\rm f}k}}$ —区域$k$ 裂缝系统拟压力，MPa²/(mPa·s)；

$r$ —半径，m；

$a$ —形状因子，无因次；

$K_{{\rm m}}$ —基质系统渗透率，D；

${{\psi }_{{\rm m}k}}$ —区域$k$ 基质系统拟压力，MPa²/(mPa·s)；

$\mu$ —黏度，Pa·s；

${{C}_{{\rm tf}}}$ —裂缝系统综合压缩系数，MPa^-1；

$t$ —时间，s；

${{C}_{{\rm tm}}}$ —基质系统综合压缩系数，MPa^-1；

${{\psi }_{{\rm i}}}$ —拟原始地层压力，MPa²/(mPa·s)；

$M$ —环状区域数；

$C$ —系数，无因次；

$r_{\rm e}$ —泄油半径，m；

$r_{\rm w}$ —井眼半径，m；

$S$ —表皮系数，无因次。

裂缝系统和基质系统的拟压力可以分别表示为

$ {\psi _{{\text{f}}k}} = \int_{{p_0}}^{{p_{{\text{f}}k}}} {\frac{{2p}}{{\mu Z}}{\text{d}}p}, \;\;\;\;\;\;{r_{k-1}} \leqslant r \leqslant {r_k} $

(13)

$ {\psi _{{\text{m}}k}} = \int_{{p_0}}^{{p_{{\text{m}}k}}} {\frac{{2p}}{{\mu Z}}{\text{d}}p}, \;\;\;\;\;\;{r_{k-1}} \leqslant r \leqslant {r_k} $

(14)

式中：$p_{{\rm 0}}$ —参考压力，MPa，为方便运算可取为0；

$p$ —压力，MPa；

${{{p}_{{\rm f}k}}}$ —区域$k$ 裂缝系统的压力，MPa;

$Z$ —压缩因子，无因次；

${{{p}_{{\rm m}k}}}$ —区域$k$ 基质系统的压力，MPa。

3 数学模型的求解

为简化数学模型及求解方便，定义下列无因次变量^[21]。

区域$k$ 无因次裂缝系统拟压力

$ {{\psi }_{{\rm Df}k}}=\dfrac{78.489{{K}_{{\rm f}k}}h}{T{{q}_{{\rm sc}}}}\left( {{\psi }_{{\rm i}}}-{{\psi }_{{\rm f}k}} \right) $

(15)

式中：$\psi _{{\rm Df}k}$ —区域$k$ 内的无因次裂缝系统拟压力；

$T$ —温度，K；

$q_{\rm sc}$ —产量，m³/d。

区域$k$ 无因次基质系统拟压力为

$ {{\psi }_{{\rm Dm}k}}=\dfrac{78.489{{K}_{{\rm m}k}}h}{T{{q}_{{\rm sc}}}}\left( {{\psi }_{{\rm i}}}-{{\psi }_{{\rm m}k}} \right) $

(16)

式中：$\psi _{{\rm Dm}k}$ —区域$k$ 内无因次基质系统拟压力；

${K}_{{\rm m}k}$ —区域$k$ 基质系统渗透率，D。

无因次有效半径为

$ {{r}_{{\rm De}}}=\dfrac{r}{{{r}_{{\rm w}}}}{{{\rm e}}^{S}} $

(17)

式中：

$r_{\rm De}$ —无因次有效半径。

无因次有效界面半径为

$ {{r}_{{\rm De}k}}=\dfrac{{{r}_{k}}}{{{r}_{{\rm w}}}}{{{\rm e}}^{S}} $

(18)

式中：

$r_{{\rm De}k}$ —无因次有效界面半径。

无因次有效时间为

$ {{t}_{{\rm De}}}={\dfrac{{3.6{K_{{\rm{f}}0}}t{{\rm{e}}^{2S}}}}{{{\phi _{{\rm{fm0}}}}{C_{{\rm{tfm0}}}}\mu r_{\rm{w}}^2}}} $

(19)

式中：$t_{\rm De}$ —无因次有效时间；

$K_{\rm f0}$ —区域0内基质系统渗透率，D；

${\phi _{{\rm{fm0}}}}$ —区域0孔隙度，$\%$ ；

${C _{{\rm{tfm0}}}}$ —区域0综合压缩系数，MPa^-1。

无因次有效井筒储积系数为

$ {{C}_{{\rm De}}}=\dfrac{C{{{\rm e}}^{2S}}}{2\pi h{\phi _{{\rm{fm0}}}}{C_{{\rm{tfm0}}}}r_{{\rm w}}^{2}} $

(20)

式中：

$C_{\rm De}$ —无因次有效井筒储积系数。

区域$k$ 的弹性储容比为

$ {{\omega }_{k}}=\dfrac{{{\phi }_{{\rm f}k}}{{C}_{{\rm ft}}}}{{{\phi }_{{\rm m}k}}{{C}_{{\rm tm}}}+{{\phi }_{{\rm f}k}}{{C}_{{\rm ft}}}} $

(21)

式中：

$\omega _k$ —区域$k$ 的弹性储容比，无因次。

区域$k$ 的窜流系数为

$ {{\lambda }_{k}}=a\dfrac{{{K}_{\rm m}}}{{{K}_{{\rm f}k}}}r_{\rm w}^{2} $

(22)

式中：

$\lambda _k$ —区域$k$ 的窜流系数，无因次。

无因次的渗流数学模型为

$ \left\{ \begin{array}{*{35}{l}} \dfrac{{{\partial }^{2}}{{\psi }_{{\rm Df}k}}}{\partial r_{{\rm De}}^{2}}+\dfrac{1}{{{r}_{{\rm De}}}}\dfrac{\partial {{\psi }_{{\rm Df}k}}}{\partial {{r}_{{\rm De}}}}+{{\lambda }_{k}}{{{\rm e}}^{-2S}}\left( \dfrac{{{K}_{{\rm f}k}}}{{{K}_{{\rm m}}}}{{\psi }_{{\rm Dm}k}}-{{\psi }_{{\rm Df}k}} \right)=\\ \dfrac{{{\omega }_{k}}}{{{C}_{{\rm De}}}}\dfrac{{{K}_{{\rm f}0}}}{{{K}_{{\rm f}k}}}\dfrac{\partial {{\psi }_{{\rm Df}k}}}{\partial \left( {{{t}_{{\rm De}}}}/{{{C}_{{\rm De}}}}\; \right)}, {\kern 10pt} {{r}_{{\rm De}\left( k-1 \right)}}\leqslant {{r}_{{\rm De}}}\leqslant {{r}_{{\rm De}k}} \\ -{{\lambda }_{k}}{{{\rm e}}^{-2S}}\left( {{\psi }_{{\rm Dm}k}}-\dfrac{{{K}_{{\rm m}}}}{{{K}_{{\rm f}k}}}{{\psi }_{{\rm Df}k}} \right)=\dfrac{{\rm 1}-{{\omega }_{k}}}{{{C}_{{\rm De}}}}\dfrac{{{K}_{{\rm f}0}}}{{{K}_{{\rm f}k}}}\dfrac{\partial {{\psi }_{{\rm Dm}k}}}{\partial \left( {{{t}_{{\rm De}}}}/{{{C}_{{\rm De}}}}\; \right)}, {\kern 86pt} \\ {{r}_{{\rm De}\left( k-1 \right)}}\leqslant {{r}_{{\rm De}}}\leqslant {{r}_{{\rm De}k}} \\ {{\psi }_{{\rm Df}k}}({{r}_{{\rm De}}}, 0)={{\psi }_{{\rm Dm}k}}({{r}_{{\rm De}}}, 0)=0, {\kern 178pt}\\ k=0, 1, \cdots, M \\ \dfrac{{\rm d}{{\psi }_{{\rm wD}}}}{{\rm d}\left( {{{t}_{{\rm De}}}}/{{{C}_{{\rm De}}}}\; \right)}-{{\left. \left( {{r}_{{\rm De}}}\dfrac{\partial {{\psi }_{{\rm Df}0}}}{\partial {{r}_{{\rm De}}}} \right) \right|}_{{{r}_{{\rm De}}}=1}}=1, {\kern 156pt} \\ {{\psi }_{{\rm wD}}}={{\left. {{\psi }_{{\rm Df0}}} \right|}_{{{r}_{{\rm De}}}=1}} \\[6pt] {{\left. {{\psi }_{{\rm Df}k}} \right|}_{{{r}_{{\rm De}}}={{r}_{{\rm De}k}}}}={{\left. {{\psi }_{{\rm Df}\left( k+1 \right)}} \right|}_{{{r}_{{\rm De}}}={{r}_{{\rm De}k}}}}, {\kern 185pt}\\ k=0, 1, \cdots, M-1 \\[8pt] {{\left. \dfrac{\partial {{\psi }_{{\rm Df}k}}}{\partial {{r}_{{\rm D}}}} \right|}_{{{r}_{{\rm De}}}={{r}_{{\rm De}k}}}}=\dfrac{{{K}_{{\rm f}\left( k+1 \right)}}}{{{K}_{{\rm f}k}}}{{\left. \dfrac{\partial {{\psi }_{{\rm Df}\left( k+1 \right)}}}{\partial {{r}_{{\rm De}}}} \right|}_{{{r}_{{\rm De}}}={{r}_{{\rm De}k}}}}, {\kern 140pt} \\ k=0, 1, \cdots, M-1 \\[8pt] {{\psi }_{{\rm Df}k}}\left( {{t}_{{\rm De}}} \right)=0, {\kern 251pt} \\ {{r}_{{\rm De}}}\to \infty; k = M \\\end{array} \right. $

(23)

将式(23)进行基于$t_{{\rm De}}$ /$C_{{\rm De}}$ 的拉普拉斯变换，得到该模型在拉式空间内的表达形式

$ \left\{ \begin{array}{*{35}{l}}\dfrac{{{\rm d}^{2}}{{{\bar{\psi }}}_{{\rm Df}k}}}{{\rm d}r_{{\rm De}}^{2}}+\dfrac{1}{{{r}_{{\rm De}}}}\dfrac{{\rm d}{{{\bar{\psi }}}_{{\rm Df}k}}}{{\rm d}{{r}_{{\rm De}}}}+{{\lambda }_{k}}{{{\rm e}}^{-2S}}\left( \dfrac{{{K}_{{\rm f}k}}}{{{K}_{{\rm m}}}}{{{\bar{\psi }}}_{{\rm Dm}k}}-{{{\bar{\psi }}}_{{\rm Df}k}} \right)=\dfrac{{{\omega }_{k}}}{{{C}_{{\rm De}}}}\dfrac{{{K}_{{\rm f}0}}}{{{K}_{{\rm f}k}}}z{{{\bar{\psi }}}_{{\rm Df}k}}, {\kern 35pt} \\ {{r}_{{\rm De}\left( k-1 \right)}}\leqslant {{r}_{{\rm De}}}\leqslant {{r}_{{\rm De}k}} \\ -{{\lambda }_{k}}{{{\rm e}}^{-2S}}\left( {{{\bar{\psi }}}_{{\rm Dm}k}}-\dfrac{{{K}_{{\rm m}}}}{{{K}_{{\rm f}k}}}{{{\bar{\psi }}}_{{\rm Df}k}} \right)=\dfrac{1-{{\omega }_{k}}}{{{C}_{{\rm De}}}}\dfrac{{{K}_{{\rm f}0}}}{{{K}_{{\rm f}k}}}z{{{\bar{\psi }}}_{{\rm Dm}k}}, {\kern 108pt} \\ {{r}_{{\rm De}\left( k-1 \right)}}\leqslant {{r}_{{\rm De}}}\leqslant {{r}_{{\rm De}k}} \\ {{{\bar{\psi }}}_{{\rm Df}k}}({{r}_{{\rm De}}}, 0)={{{\bar{\psi }}}_{{\rm Dm}k}}({{r}_{{\rm De}}}, 0)=0, {\kern 173pt}\\ k=0, 1, \cdots, M \\ z{{{\bar{\psi }}}_{{\rm wD}}}-{{\left. \left( {{r}_{{\rm De}}}\dfrac{d{{{\bar{\psi }}}_{{\rm Df}0}}}{d{{r}_{{\rm De}}}} \right) \right|}_{{{r}_{{\rm De}}}=1}}=\dfrac{1}{z}, {\kern 180pt} \\ {{{\bar{\psi }}}_{{\rm wD}}}={{\left. {{{\bar{\psi }}}_{{\rm Df0}}} \right|}_{{{r}_{{\rm De}}}=1}} \\ {{\left. {{{\bar{\psi }}}_{{\rm Df}k}} \right|}_{{{r}_{{\rm De}}}={{r}_{{\rm De}k}}}}={{\left. {{{\bar{\psi }}}_{{\rm Df}\left( k+1 \right)}} \right|}_{{{r}_{{\rm De}}}={{r}_{{\rm De}k}}}}, {\kern 180pt} \\ k=0, 1, \cdots, M-1 \\[8pt] {{\left. \dfrac{{\rm d}{{{\bar{\psi }}}_{{\rm Df}k}}}{{\rm d}{{r}_{{\rm De}}}} \right|}_{{{r}_{{\rm De}}}={{r}_{{\rm De}k}}}}=\dfrac{{{K}_{{\rm f}\left( k+1 \right)}}}{{{K}_{{\rm f}k}}}{{\left. \dfrac{{\rm d}{{{\bar{\psi }}}_{{\rm Df}\left( k+1 \right)}}}{{\rm d}{{r}_{{\rm De}}}} \right|}_{{{r}_{{\rm De}}}={{r}_{{\rm De}k}}}}, {\kern 135pt} \\ k=0, 1, \cdots, M-1 \\[8pt] {{{\bar{\psi }}}_{{\rm Df}k}}\left( z \right)=0, {\kern 254pt} \\ {{r}_{{\rm De}}}\to \infty ; k = M \\ \end{array} \right. $

(24)

式中：$\bar \psi _{{\rm Df}k}$ —拉氏空间内区域$k$ 无因次裂缝系统拟压力；$\bar \psi _{{\rm Dm}k}$ —拉氏空间内区域$k$ 无因次基质系统拟压力。

式(24)中，${{\bar{\psi }}_{{\rm Df}k}}$ 的通解为

$ {{\bar{\psi }}_{{\rm Df}k}}={{A}_{k}}{{I}_{0}}\left( {{r}_{{\rm De}k}}\sqrt{{{S}_{k}}\left( z \right)} \right)+{{B}_{k}}{{K}_{0}}\left( {{r}_{{\rm De}k}}\sqrt{{{S}_{k}}\left( z \right)} \right), \quad k=0, 1, \cdots, M $

(25)

式中：$A_k$ ，$B_k$ —系数，($k = 0, 1, \cdots, M$ )。

对式(25)求导，有

$ \dfrac{{\rm d}{{{\bar{\psi }}}_{{\rm Df}k}}}{{\rm d}{{r}_{{\rm De}}}}=\sqrt{{{S}_{k}}\left( z \right)}{{A}_{k}}{{I}_{1}}\left( {{r}_{{\rm De}k}}\sqrt{{{S}_{k}}\left( z \right)} \right)-\sqrt{{{S}_{k}}\left( z \right)}{{B}_{k}}{{K}_{1}}\left( {{r}_{{\rm De}k}}\sqrt{{{S}_{k}}\left( z \right)} \right),\\ \quad k=0, 1, \cdots, M $

(26)

将式(26)、(25)代入式(24)，求取$A_{k}$ 和$B_{k }$ ($k=0, 1, \cdots, M$ )

$ \left\{ {\begin{array}{*{20}{l}} {z{{\bar \psi }_{{\rm{wfD}}}} - \sqrt {{S_0}\left( z \right)} {A_0}{I_1}\left( {\sqrt {{S_0}\left( z \right)} } \right) + \sqrt {{S_0}\left( z \right)} {B_0}{K_1}\left( {\sqrt {{S_0}\left( z \right)} } \right) = \dfrac{1}{z}} \\ {{{\bar \psi }_{{\rm{wfD}}}} = {A_0}{I_0}\left( {\sqrt {{S_0}\left( z \right)} } \right) + {B_0}{K_0}\left( {\sqrt {{S_0}\left( z \right)} } \right)}&{}&{}\\ {A_k}{I_0}\left( {{r_{{\rm{De}}k}}\sqrt {{S_k}\left( z \right)} } \right) + {B_k}{K_0}\left( {{r_{{\rm{De}}k}}\sqrt {{S_k}\left( z \right)} } \right) = {A_{k + 1}}{I_0}\left( {{r_{{\rm{De}}k}}\sqrt {{S_{k + 1}}\left( z \right)} } \right) + \\{\kern 20pt}{B_{k + 1}}{K_0}\left( {{r_{{\rm{De}}k}}\sqrt {{S_{k + 1}}\left( z \right)} } \right), \\ {k = 1, 2, \cdots, M - 1}\\ {A_k}\sqrt {{S_k}\left( z \right)} {I_1}\left( {{r_{{\rm{De}}k}}\sqrt {{S_k}\left( z \right)} } \right) - {B_k}\sqrt {{S_k}\left( z \right)} {K_1}\left( {{r_{{\rm{De}}k}}\sqrt {{S_k}\left( z \right)} } \right) = \\ \dfrac{{{k_{{\rm{f}}\left( {k + 1} \right)}}}}{{{k_{{\rm{f}}k}}}} {A_{k + 1}}\cdot \\{\kern 20pt}\sqrt {{S_k}\left( z \right)} {I_1}\left( {{r_{{\rm{De}}k}}\sqrt {{S_k}\left( z \right)} } \right) - \dfrac{{{K_{{\rm{f}}\left( {k + 1} \right)}}}}{{{K_{{\rm{f}}k}}}}{B_{k + 1}}\sqrt {{S_k}\left( z \right)} {K_1}\left( {{r_{{\rm{De}}k}}\sqrt {{S_k}\left( z \right)} } \right), \\ {k = 1, 2, \cdots, M - 1}\\ {{A_M} = 0} \end{array}} \right. $

(27)

式中：$S_{{ k}}$ —拉普拉斯变量函数；

$I_0$ —零阶第一类贝塞尔函数；

$I_1$ —一阶第一类贝塞尔函数；

$K_0$ —零阶第二类贝塞尔函数；

$K_1$ —一阶第二类贝塞尔函数。

式(27)中，参数无因次有效界面半径$r_{{\rm De}k}$ 、渗透率比$K _{{ {\rm f}(k+1)}}$ /$K _{{\rm f}k}$ 和函数$S_{{ k}}$ 可由树状分形网络参数直接表示。

由式(3)和式(17)，得到无因次界面有效半径的表达式

$ {{r}_{{\rm De}k}}=\dfrac{{{l}_{0}}}{{{r}_{{\rm w}}}}\left[1+\dfrac{\alpha \left( 1-{{\alpha }^{k}} \right)\cos \theta }{1-\alpha } \right]{{{\rm e}}^{S}} $

(28)

由式(6)，渗透率比的表达式为

$ \dfrac{{{K}_{{\rm f}\left( k+1 \right)}}}{{{K}_{{\rm f}k}}}=\left\{ \begin{array}{*{35}{l}}n{{\beta }^{2}}\cos \theta, {\kern 20pt}k=0 \\n{{\beta }^{2}}, {\kern 43pt} k>0 \end{array} \right. $

(29)

$S_{{ k}}$ ($z$ )的表达式为

$ {{S}_{k}}\left( z \right)=\dfrac{\dfrac{{{K}_{{\rm f}0}}}{{{K}_{{\rm f}k}}}{{\lambda }_{k}}\left( 1-{{\omega }_{k}} \right)z}{\dfrac{{{K}_{{\rm f}0}}}{{{K}_{{\rm f}k}}}\left( 1-{{\omega }_{k}} \right){{{\rm e}}^{2S}}z+{{\lambda }_{k}}{{C}_{{\rm De}}}}+ \\ {\kern 40pt}\dfrac{{{K}_{{\rm f}0}}}{{{K}_{{\rm f}k}}}\dfrac{{{\omega }_{k}}}{{{C}_{{\rm De}}}}z $

(30)

其中

$ \dfrac{{{K}_{{\rm f}0}}}{{{K}_{{\rm f}k}}}={{\left( n{{\beta }^{2}} \right)}^{-k}} $

(31)

$ {{\lambda }_{k}}=\left\{ \begin{array}{*{35}{l}}\dfrac{32a{{K}_{{\rm m}}}r_{{\rm w}}^{2}}{Nd_{0}^{2}},&{\quad}k=0 \\ \dfrac{32a{{K}_{{\rm m}}}r_{{\rm w}}^{2}}{N{{n}^{k}}{{\beta }^{2k}}d_{0}^{2}\cos \theta },& {\quad}k>0 \\ \end{array} \right. $

(32)

$\omega_{{k}}$ 由式(10)，式(11)，式(21)联合求取。

4 典型曲线及拟压力动态特征

对式(27)进行求解，可得到拉氏空间内无因次井底拟压力${{\bar{\psi }}_{{\rm wfD}}}$ 。利用Stehfest数值反演方法对拉氏空间内无因次井底拟压力${{\bar{\psi }}_{{\rm wfD}}}$ 进行拉普拉斯逆变换，得到实空间内的裂缝性气藏气井无因次井底拟压力的数值解，从而绘制树状分形网络裂缝性气藏的拟压力动态特征典型曲线，如图 2所示。

树状分形网络裂缝性气藏中，各区域裂缝系统渗透率、裂缝体积相等的树状分形裂缝性气藏典型曲线，即满足式(13)。裂缝系统分形网络的参数取值为$\alpha$ =1, $\beta$ =0.707，$\theta$ = 1，$N=4$ ，$M$ =10，$l_{{\rm 0}}$ =10 m，$d_{{\rm 0}}=0.05$ m，其他参数取值为$r_{{\rm w}}=0.1$ m，$S= -1$ ，$C$ =0.0001，$C_{{\rm mt}}$ =2.2×10$^{{\rm -2}}$ MPa^-1，$k_{{\rm m}}$ =5 mD，$h$ =10 m，$\phi_{{\rm mp}}$ =0.1，$C_{{\rm ft}}$ =1×10$^{{\rm -2}}$ MPa^-1。从图 2可以看出，流动阶段可被划分为4个阶段：第1个典型曲线阶段为井储阶段，无因次拟压力和拟压力导数随无因次时间的增大而迅速增大，其曲线表现为倾斜为$\pi$ /4的直线。第2个典型曲线阶段为裂缝系统径向流阶段，无因次拟压力随无因次时间的增大继续增大，增加速度不断减小至零；无因次拟压力导数随无因次时间的增大继续增大，到达顶点后不断减小，直至稳定至0.5，该阶段往往被井储流动阶段所掩盖，在典型曲线上一般不容易观察到；第3个典型曲线阶段为基质系统向裂缝系统的窜流阶段，无因次拟压力导数随无因次时间的增大迅速减小，到达底端后迅速增大，拟压力导数曲线表现为先下降，然后又上升，形成“V”型曲线；第4个典型曲线阶段为整个系统的径向流阶段，无因次拟压力随无因次时间的增大而缓慢增大，无因次拟压力导数恒定为0.5。

图 3为长度比$\alpha$ 影响下树状分形网络裂缝性气藏的拟压力动态特征典型曲线。

由图 3可见，长度比$\alpha$ 决定拟压力导数曲线中第3个典型曲线阶段“V”型曲线的宽度和深度，$\alpha$ 越大，拟压力导数所形成的“V”型曲线越宽、越深，表明第3个典型曲线阶段持续的时间越长，第2个典型曲线阶段持续的时间就越短，甚至可能被井筒储集效应所掩盖。随着长度比$\alpha$ 的减小，拟压力导数所形成的“V”型曲线变窄、变浅，第3个典型曲线阶段持续的时间缩短，第2个典型曲线阶段流动阶段持续的时间增长，且逐渐表现出裂缝系统的径向流动。说明长度比$\alpha$ 主要影响第3个典型曲线阶段。

图 4为直径比$\beta$ 影响下树状分形网络裂缝性气藏的拟压力动态特征典型曲线。

由图 4可见，$\beta$ 越大，拟压力导数曲线的整体位置向下方移动。随着直径比$\beta$ 的增大，拟压力导数曲线的位置越来越靠下。此外，当$\beta>0.707$ 时，拟压力导数曲线中第4个典型曲线阶段表现为小于0.5的水平线，反之，第4个典型曲线阶段表现为大于0.5的水平线。说明直径比$\beta$ 主要影响除井筒储集阶段外的所有阶段。

图 5为分叉角度$\theta$ 影响下树状分形网络裂缝性气藏的拟压力动态特征典型曲线。

由图 5可见，分叉角度$\theta$ 决定拟压力导数曲线中第3个典型曲线阶段“V”型曲线的宽度和深度，$\theta$ 越小，拟压力导数形成的“V”型曲线越宽、越深；从2、3阶段持续时间上看，$\theta$ 越小，第2个典型曲线阶段越短，第3个典型曲线阶段越长。随着分叉角度$\theta$ 的增大，拟压力导数形成的“V”型曲线变窄变浅，第3个典型曲线阶段持续的时间缩短，裂缝系统流动阶段持续的时间增长，且逐渐表现出裂缝系统的径向流动。说明分叉角度$\theta$ 主要影响第3个典型曲线阶段。

图 6、图 7分别为$\beta$ =0.65和0.75时，总分叉级数$M$ 影响下树状分形网络裂缝性气藏的拟压力动态特征典型曲线。

由图 6、图 7可见，总分叉级数$M$ 决定拟压力导数曲线中第4个典型曲线阶段。当$\beta>0.707$ 时，$M$ 越大，拟压力导数曲线在第4个典型曲线阶段就越靠上，甚至可能出现圆形封闭边界反映。随着总分叉级数$M$ 的减小，拟压力导数曲线在第4个典型曲线阶段就越来越靠下。当$\beta$ < 0.707时，$M$ 越大，拟压力导数曲线在第4个典型曲线阶段就越靠下，随着总分叉级数$M$ 的减小，拟压力导数曲线在第4个典型曲线阶段就越来越靠上。说明总分叉级数$M$ 主要影响第4个典型曲线阶段。

图 8为初始裂缝数$N$ 影响下的树状分形网络裂缝性气藏的拟压力动态特征典型曲线。

由图 8可见，初始裂缝数$N$ 决定拟压力导数曲线中第3个典型曲线阶段“V”型曲线的位置、宽度和深度，$N$ 越小，拟压力导数形成的“V”型曲线就越靠左，且越宽越深，第3个典型曲线阶段发生的时间就越早，第2个典型曲线阶段持续的时间就越短。随着初始裂缝数$N$ 的增加，拟压力导数形成的“V”型曲线变窄变浅，向右上方移动，且逐渐表现出裂缝系统的径向流动。说明初始裂缝数$N$ 主要影响第3个典型曲线阶段。

5 结语

利用树状分形网络模拟裂缝系统，并将树状分形网络嵌入到基质中，提出了基于树状分形网络的裂缝性气藏试井模型，作出了树状分形网络的裂缝性气藏拟压力动态特征典型曲线。分析了长度比$\alpha$ 、直径比$\beta$ 、分叉角度$\theta$ 、总分叉级数$M$ 和初始裂缝数$N$ 对拟压力动态特征的影响，长度比$\alpha$ 、分叉角度$\theta$ 和初始裂缝数$N$ 主要影响动态特征典型曲线中第3个典型曲线阶段，总分叉级数$M$ 主要影响第4个典型曲线阶段，直径比$\beta$ 影响除井筒储集阶段外的所有阶段。分析基于树状分形网络的裂缝性气藏试井模型典型曲线的影响因素表明，树状分形网络能够很好地模拟裂缝性气藏的裂缝系统。