引言

致密气作为一种非常规油气资源，表现出巨大的勘探开发前景，其探明储量占到中国天然气储量的33%。2019年，中国天然气产量为1 777$\times$10$^8$ m$^3$，其中，致密气产量为450$\times$10$^8$ m$^3$，占比达25.3%，致密气成为中国天然气增储上产的重要领域。目前，现场广泛应用多级压裂水平井技术开发致密气藏，通过人造裂缝改善近井渗流环境，提高流动能力，从而提高开发效果^[1-4]。

压裂井压后生产动态的预测模型大致分为3类。第一类是以线性流模型为主的解析模型^[5-8]，其优势是物理模型意义明确、求解简单、计算速度快。但该类模型大多未考虑复杂渗流机理，而且对于复杂裂缝形态刻画多以均匀对称为主，裂缝间干扰采用设置不渗透边界进行分区渗流研究，因此，该方法存在较大的局限性。第二类是在源函数基础上建立的半解析模型^[9-16]，该类方法通过划分微元体离散裂缝，建立储层及裂缝系数矩阵，迭代计算稀疏矩阵进行求解，该类方法考虑了裂缝的形态对生产动态的影响，但在裂缝三维非平面研究方面存在局限。第三类是数值模型^[17-23]，该类方法利用网格离散技术和数值计算方法，可以准确地描述压裂后复杂缝网形态及储层渗流特征，常用方法有基于正交网格的有限差分法、基于非结构网格的有限元法和控制体积-有限元法。非结构网格降低了正交网格取向性，避免局部裂缝网格加密导致的不收敛，适合复杂形态描述。控制体积-有限元法则结合有限元和有限体积法的优点，在局部守恒和饱和度计算时能保证较高的精度并能提高计算速度。前人有关致密气藏多级压裂水平井生产动态数值模拟主要考虑垂直裂缝条数、间距、长度等因素对压后生产动态的影响，可实现传统垂直裂缝压裂水平井稳态和非稳态产量计算功能，但缺少对非对称空间分布压裂裂缝的定量描述，且未采用基于三维非结构网格的控制体积有限元方法对复杂裂缝形态和流动特征进行准确表征。本文综合考虑致密气储层应力敏感、气水两相非线性流动特征和压裂裂缝气体高速非达西流动规律，建立双重连续介质-离散裂缝耦合渗流数学模型，采用非结构三维四面体网格和控制体积-有限元法建立全隐式数值模型，分析含水饱和度、应力敏感系数、压裂裂缝压开程度和空间非对称分布等对致密气藏多级压裂水平井生产动态的影响。

1 物理模型与假设条件

致密气藏具有低孔、低渗特征，采用传统的直井、水平井等开采方式难以获得工业产能。通过水力压裂改造，在水平井周围形成具有一定长度、高度和导流能力的高渗透压裂裂缝，改善近井地带流动特征，增大渗流体积。模型及空间网格离散如图 1所示。

渗流物理模型基本假设如下：(1)渗流过程为气水两相等温渗流；(2)采用双重介质模型描述天然裂缝广泛发育的致密储层中的流动，流体由基质向天然裂缝系统的流动为拟稳态窜流；(3)压裂缝具有有限导流和空间非对称分布特征，各条裂缝在长度、开度、方位角和裂缝高度等属性上存在差异；(4)流体由基质系统流向天然裂缝系统，并通过压裂缝流入水平井筒进而产出。

2 数学模型 2.1 天然裂缝应力敏感效应

在生产过程中储层天然裂缝会随着储层压力的改变，发生闭合或渗透率降低等现象，这种现象在致密气藏中表现尤为明显。本文采用达西渗流公式对天然裂缝系统中的流动特征进行描述

$ v = - \dfrac{{{K_{{\rm{fe}}}}}}{{\mu _{\rm g} B_{\rm g}}}\nabla p $

(1)

考虑应力敏感效应的天然裂缝修正渗透率$K_{\rm{fe}}$计算式为^[14]

$ {K_{{\rm{fe}}}} = {K_{{\rm{fi}}}}{{\mathop{\rm e}\nolimits} ^{ - {d_{\rm{f}}}({p_{{\rm{fi}}}} - {p_{\rm{f}}})}} $

(2)

2.2 气相表观渗透率

气体在通过致密储层微纳尺度孔隙时表现出滑脱、分子扩散等非达西流动特征。采用直径、长度、大小不等的毛细管束对微纳尺度孔隙进行表征，基于文献[25]，推导单毛细管中耦合黏性流和分子扩散的摩尔质量通量

$ J\left( \lambda \right) =\\{\kern 30pt} {\dfrac{{{ {\rm{ \mathsf{ π} }}}{\lambda ^4}{\rho _{\rm{g}}}}}{{128{\mu _{\rm{g}}}M}}\dfrac{{\Delta p}}{{{L_{\rm{t}}}}}} \dfrac{{{ {\rm{ \mathsf{ π} }}}{\lambda ^3}}}{{12}}{\rm{ + }} {\dfrac{{{K{ n}}}}{{{K{ n}} + 1}}\sqrt {\dfrac{{8Z}}{{{ {\rm{ \mathsf{ π} }}}M{\rm R}T}}} \dfrac{p}{Z}{c_{\rm{g}}}\dfrac{{\Delta p}}{{{L_{\rm{t}}}}}} $

(3)

基于分形渗流理论，将式(3)中的单根毛细管质量通量$J(\lambda)$沿着最小和最大毛细管直径($\lambda_{\rm{min}}\sim \lambda_{\rm{max}}$)积分，获取岩样3D空间中总数量为$N$的弯曲毛细管的总摩尔质量流量

$ {Q_{\rm{T}}} = - \int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {J\left( \lambda \right){\rm{d}}N} $

(4)

对比达西渗流摩尔质量流量公式，可得到致密气藏考虑黏性流动、分子扩散的真实气体分形表观渗透率

$ {K_{{\rm{am}}}} = {\dfrac{{{ {\rm{ \mathsf{ π} }}}\lambda _{\max }^{3 + {D_{\rm{T}}}}}}{{128L_0^{{D_{\rm{T}}} + 1}}}\dfrac{{{D_{\rm{f}}}}}{{3 - {D_{\rm{f}}} + {D_{\rm{T}}}}}} + {\dfrac{{{{10}^{ - 9}} { {\rm{ \mathsf{ π} }}}{\mu _{\rm{g}}}}}{{12}}\sqrt {\dfrac{{8Z{\rm{R}}T}}{{{ {\rm{ \mathsf{ π} }}}M}}} \dfrac{{{D_{\rm{f}}}\lambda _{\max }^{{D_{\rm{f}}}}}}{{L_0^{{D_{\rm{T}}} + 1}}}{c_{\rm{g}}}\int_{{\lambda _{\min }}}^{{\lambda _{\max }}} {\dfrac{{{K{ n}}}}{{{K{ n}} + 1}}{\lambda ^{1 - {D_{\rm{f}}} + {D_{\rm{T}}}}}{\rm{d}}\lambda } } $

(5)

2.3 致密储层渗流模型

基质系统与天然裂缝系统间的气相拟稳态窜流公式^[17]

$ q_{\rm{g, fm}} = \dfrac{{\alpha {K_{{\rm{am}}}}{K_{{\rm{rgm}}}}}}{{{\mu _{{\rm{gm}}}}{B_{{\rm{gm}}}}}}\left( {{p_{{\rm{gm}}}} - {p_{{\rm{gf}}}}} \right) $

(6)

基质系统与天然裂缝系统间的水相拟稳态窜流公式^[17]

$ {q_{{\rm{w}}, {\rm{fm}}}} = \dfrac{{\alpha {K_{\rm{m}}}{K_{{\rm{rwm}}}}}}{{{\mu _{{\rm{wm}}}}{B_{{\rm{wm}}}}}}\left( {{p_{{\rm{wm}}}} - {p_{{\rm{wf}}}}} \right) $

(7)

(1) 基质系统

气相渗流控制方程

$ \nabla \cdot \left( {\dfrac{{{K_{{\rm{am}}}}{K_{{\rm{rgm}}}}}}{{{\mu _{{\rm{gm}}}}{B_{{\rm{gm}}}}}}\nabla {p _{{\rm{gm}}}}} \right) - {q_{{\rm{g, fm}}}} = \dfrac{\partial }{{\partial t}}\left( {\dfrac{{{\phi _{_{\rm{m}}}}{S_{{\rm{gm}}}}}}{{{B_{{\rm{gm}}}}}}} \right) $

(8)

水相渗流控制方程

$ \nabla \cdot \left( {\dfrac{{{K_{\rm{m}}}{K_{{\rm{rwm}}}}}}{{{\mu _{{\rm{wm}}}}{B_{{\rm{wm}}}}}}\nabla {p _{{\rm{wm}}}}} \right) - {q_{{\rm{w, fm}}}} = \dfrac{\partial }{{\partial t}}\left( {\dfrac{{{\phi _{\rm{m}}}{S_{{\rm{wm}}}}}}{{{B_{{\rm{wm}}}}}}} \right) $

(9)

(2) 裂缝系统

气相渗流控制方程

$ \nabla \cdot \left( {\dfrac{{{K_{{\rm{fe}}}}{K_{{\rm{rgf}}}}}}{{{\mu _{{\rm{gf}}}}{B_{{\rm{gf}}}}}}\nabla {p _{{\rm{gf}}}}} \right) + {q_{{\rm{g, fm}}}} - {q_{{\rm{g, Ff}}}} = \\[5pt]{\kern 40pt} \dfrac{\partial }{{\partial t}}\left( {\dfrac{{{\phi _{\rm{f}}}{S_{{\rm{gf}}}}}}{{{B_{{\rm{gf}}}}}}} \right) $

(10)

水相渗流控制方程

$ \nabla \cdot \left( {\dfrac{{{K_{{\rm{fe}}}}{K_{{\rm{rwf}}}}}}{{{\mu _{{\rm{wf}}}}{B_{{\rm{wf}}}}}}\nabla {p _{{\rm{wf}}}}} \right) + {q_{{\rm{w, fm}}}} - {q_{{\rm{w, Ff}}}}= \\[5pt]{\kern 40pt} \dfrac{\partial }{{\partial t}}\left( {\dfrac{{{\phi _{\rm{f}}}{S_{{\rm{wf}}}}}}{{{B_{{\rm{wf}}}}}}} \right) $

(11)

2.4 压裂缝渗流模型

气体在通过高渗透压裂裂缝时表现出高速非达西流动特征。考虑含水饱和度影响的扩展Forchheimer高速非达西公式为^[21]

$ \nabla {p_{{\rm{gF}}}} = \dfrac{{{\mu _{{\rm{gF}}}}}}{{{K_{\rm{F}}}{K_{{\rm{rgF}}}}}}v + \beta {\rho _{{\rm{gF}}}}{v^2} $

(12)

同时，定义一个气相高速非达西修正系数$\varepsilon_n$，则满足^[21]

$ \left\{ \begin{array}{l} {v_n} = {\varepsilon _n}{v_{{\rm{d}}n}}, n=x, y, z \\[5pt] {\varepsilon _n} = \dfrac{2}{{1 + \sqrt {1 + 4{\rho _{\rm{g}}}{\beta }{K_{\rm F}}{K_{{\rm{rgF}}}}{v_{{\rm{d}}n}}/{\mu _{\rm{gF}}}} }} \end{array}\right. $

(13)

由此，压裂缝中气相渗流控制方程为

$ \nabla \cdot \left( {{\varepsilon _n}\dfrac{{{K_{\rm{F}}}{K_{{\rm{rgF}}}}}}{{{\mu _{{\rm{gF}}}}{B_{{\rm{gF}}}}}}\nabla {p _{{\rm{gF}}}}} \right) + {q_{{\rm{g, Ff}}}} + {q_{{\rm{gsct}}}} =\\[5pt]{\kern 40pt} \dfrac{\partial }{{\partial t}}\left( {\dfrac{{{\phi _{\rm{F}}}{S_{{\rm{gF}}}}}}{{{B_{{\rm{gF}}}}}}} \right) $

(14)

压裂缝中水相渗流控制方程为

$ \nabla \cdot \left( {\dfrac{{{K_{\rm{F}}}{K_{{\rm{rwF}}}}}}{{{\mu _{\rm{w}}}{B_{\rm{w}}}}}\nabla {p _{{\rm{wF}}}}} \right) + {q_{{\rm{w, Ff}}}} + {q_{{\rm{wsct}}}} = \\[5pt]{\kern 40pt} \dfrac{\partial }{{\partial t}}\left( {\dfrac{{{\phi _{\rm{F}}}{S_{{\rm{wF}}}}}}{{{B_{\rm{w}}}}}} \right) $

(15)

2.5 双重介质-离散裂缝耦合模型

式(1)$\sim$式(15)构成了致密气藏多级压裂水平井的综合渗流方程组，基于离散裂缝模型原理^[20]，通过式(16)将水力压裂缝在“降维”空间的渗流方程统一到双重连续介质渗流空间中，构成了全储层区域渗流方程体系

$ \int_{{\varOmega _{\rm{d}}}} F{\rm{ d}}{{\varOmega _{\rm{d}}}} {\rm{ = }}\int_{{\varOmega _{\rm{f}}}} F {\rm{ d}}{\varOmega _{\rm{f}}} + {w_{\rm{F}}} \int_{{\varOmega _{\rm{F}}}} F {\rm{ d}}{\varOmega _{\rm{F}}} $

(16)

初始条件满足

$ \left\{ \begin{array}{l} {p_{\rm{f}}} = {p_{{\rm{fi}}}}\\ {p_{\rm{m}}} = {p_{{\rm{mi}}}} \end{array} \right. $

(17)

封闭外边界条件为

$ \left\{ \begin{array}{l} \dfrac{{\partial {p_{\rm{f}}}}}{{\partial n}} = 0\\[5pt] \dfrac{{\partial {p_{\rm{m}}}}}{{\partial n}} = 0 \end{array} \right. $

(18)

3 模型求解

基于储层三维四面体网格和压裂裂缝三角形网格，采用控制体积-有限元法进行数值求解^[21]。首先采用线性插值函数，用四面体顶点值来插值逼近四面体单元平均值，得到

$ \begin{array}{l} \overline p \left( {x, y} \right) = \sum\limits_{v = 1}^4 {{N_v}\left( {x, y} \right){p_v}} \\ {\overline S _{\rm{w}}}\left( {x, y} \right) = \sum\limits_{v = 1}^4 {{N_v}\left( {x, y} \right){S_{{\rm{w}}v}}} \end{array} $

(19)

插值形函数$N_v$可表示为

$ {N_v} = \dfrac{1}{{6V}}\left( {{a_v} + {b_v}x + {c_v}y + {d_v}z} \right), v=1, 2, 3, 4 $

(20)

同时，采用牛顿-拉夫逊迭代格式对$m+1$时间步下的变量$p^{m+1}$和$S_{\rm w}^{m+1}$采用第$k$迭代步的值来近似逼近，可以得到

$ p^{m+1} \approx {p^{k + 1}} = {p^k} + \delta {p^k} $

(21)

$ S_{\rm{w}}^{m+1} \approx S_{\rm{w}}^{k+1} = S_{\rm{w}}^k + \delta S_{\rm{w}}^k $

(22)

选取储层空间任意一个网格单元，构建包含致密储层渗流和压裂裂缝流动的全隐式迭代单元矩阵

$ \left[ \begin{array}{l} {\pmb{k}_{11}} \ \ {\pmb{k}_{12}} \ \ {\pmb{k}_{13}} \ \ {\pmb{k}_{14}}\\ {\pmb{k}_{21}} \ \ {\pmb{k}_{22}} \ \ {\pmb{k}_{23}} \ \ {\pmb{k}_{24}}\\ {\pmb{k}_{31}} \ \ {\pmb{k}_{32}} \ \ {\pmb{k}_{33}} \ \ {\pmb{k}_{34}}\\ {\pmb{k}_{41}} \ \ {\pmb{k}_{42}} \ \ {\pmb{k}_{43}} \ \ {\pmb{k}_{44}} \end{array} \right]\left[ \begin{array}{l} {δ}\pmb{P}_{{\rm{gf}}}^k\\ {δ}\pmb{P}_{{\rm{gm}}}^k\\ {δ}\pmb{S}_{{\rm{wf}}}^k\\ {δ}\pmb{S}_{{\rm{wm}}}^k \end{array} \right] = \left[ \begin{array}{l} \pmb{R}_{{\rm{gf}}}^k\\ \pmb{R}_{{\rm{gm}}}^k\\ \pmb{R}_{{\rm{wf}}}^k\\ \pmb{R}_{{\rm{wm}}}^k \end{array} \right] $

(23)

其中

${\pmb{k}_{11}}{\rm{ = }}\pmb{T}_{{\rm{gf}}}^k + \pmb{T}_{{\rm{gf}}}^k {\rm d}{p_{{\rm{gf}}}}\pmb{P}_{{\rm{gf}}}^k - \pmb{W}_{\rm{g}}^k -\\ {\hspace{4em}} \dfrac{{\pmb{N}_{{\rm{gf}}}^k {\rm d}{p_{{\rm{gf}}}}}}{{\Delta t}} + \dfrac{{\pmb{N}_{{\rm{gf}}{{\rm{S}}_{{\rm{wf}}}}}^k {\rm d}{p_{{\rm{gf}}}}}}{{\Delta t}}\\ {\pmb{k}_{12}}{\rm{ = }}\pmb{W}_{\rm{g}}^k {\rm d}{p_{{\rm{gm}}}} + \pmb{W}_{\rm{g}}^k\\ {\pmb{k}_{13}} = \pmb{T}_{{\rm{gf}}}^k {\rm d}{S_{{\rm{wf}}}}\pmb{P}_{{\rm{gf}}}^k + \dfrac{{\pmb{N}_{{\rm{gf}}}^k}}{{\Delta t}}\\ {\pmb{k}_{14}}{\rm{ = }}\pmb{W}_{\rm{g}}^k {\rm d}{\pmb{S}_{{\rm{wm}}}}\\ {\pmb{k}_{21}}{\rm{ = }}\pmb{W}_{\rm{g}}^k\\ {\pmb{k}_{22}}{\rm{ = }}\pmb{T}_{{\rm{gm}}}^k + \pmb{T}_{{\rm{gm}}}^k {\rm d}{p_{{\rm{gm}}}}\pmb{P}_{{\rm{gm}}}^k - \pmb{W}_{\rm{g}}^k {\rm d}{p_{{\rm{gm}}}} - \pmb{W}_{\rm{g}}^k -\\ {\hspace{4em}} \dfrac{{\pmb{N}_{{\rm{gm}}}^k {\rm d}{p_{{\rm{gm}}}}}}{{\Delta t}} + \dfrac{{\pmb{N}_{{\rm{gm}}{{\rm{S}}_{{\rm{wm}}}}}^k {\rm d}{p_{{\rm{gm}}}}}}{{\Delta t}}\\ {\pmb{k}_{23}} = 0\\ {\pmb{k}_{24}}{\rm{ = }}\pmb{T}_{{\rm{gm}}}^k {\rm d}{S_{{\rm{wm}}}}\pmb{P}_{{\rm{gm}}}^k - \pmb{W}_{\rm{g}}^k {\rm d}{S_{{\rm{wm}}}} + \dfrac{{\pmb{N}_{{\rm{gm}}}^k}}{{\Delta t}}\\ {\pmb{k}_{31}}{\rm{ = }}\pmb{T}_{{\rm{wf}}}^k - \pmb{W}_{\rm{w}}^k\\ {\pmb{k}_{32}}{\rm{ = }}\pmb{W}_{\rm{w}}^k\\ {\pmb{k}_{33}} = \pmb{T}_{{\rm{wf}}}^k {\rm d}{S_{{\rm{wf}}}}\pmb{P}_{{\rm{gf}}}^k - \dfrac{{\pmb{N}_{{\rm{wf}}}^k}}{{\Delta t}}\\ {\pmb{k}_{34}}{\rm{ = }}\pmb{W}_{\rm{w}}^k {\rm d}{S_{{\rm{wm}}}}\\ {\pmb{k}_{41}}{\rm{ = }}\pmb{W}_{\rm{w}}^k\\ {\pmb{k}_{42}}{\rm{ = }}\pmb{T}_{{\rm{wm}}}^k - \pmb{W}_{\rm{w}}^k\\ {\pmb{k}_{43}} = 0\\ {\pmb{k}_{44}}{\rm{ = }}\pmb{T}_{{\rm{wm}}}^k {\rm d}{S_{{\rm{wm}}}}\pmb{P}_{{\rm{gm}}}^k - \dfrac{{\pmb{N}_{{\rm{wm}}}^k}}{{\Delta t}}\\ $

假设网格剖分形成$N$个节点，即$N$个控制体积网格，通过叠加每一个单元矩阵，最终形成$4N\times 4N$的全区域方程组来求解4$N$个包含$\delta p_{\rm{gf1}}$, …, $\delta p_{{\rm{gf}}N}$, $\delta p_{\rm{gm1}}$, …, $\delta p_{{\rm{gm}}N}$, $\delta S_{\rm{wf1}}$, …, $\delta S_{{\rm{wf}}N}$, $\delta S_{\rm{wm1}}$, …, $\delta S_{{\rm{wm}}N}$的变量，整体求解矩阵表达式为

$ {\left[ \boldsymbol{{k}} \right]_{4N \times 4N}}{\left[ {\delta \boldsymbol{{X}}} \right]_{4N \times 1}} = {\left[ \boldsymbol{{R}} \right]_{4N \times 1}} $

(24)

同时，由离散裂缝模型将压裂裂缝“降维”处理成宽度为$w_{{\rm F}i}$的裂缝面，结合拟稳态流动方程，将储层中流体进入压裂裂缝的流动看作径向流入高度为$w_{{\rm F}i}$，底部半径为$r_{0i}$的等效垂直井的流动过程。其等效半径$r_{0i}$为

$ {r_{0i}} = \sqrt {\dfrac{{{A_i}}}{{ {\rm{ \mathsf{ π} }}}}} $

(25)

气-水两相产量在$m+1$时间步下的拟稳态井模型满足

$ \left\{\begin{array}{l} q_{{\rm{gsct}}}^{m + 1} = \sum\limits_{i = 1}^{{N_{\rm F}}} {J_{{\rm{g}}i}^{m + 1}{{\left( {{p_{{\rm{bh}}}} - {p_{{\rm{ave}}i}}} \right)}^{m + 1}}} \\ q_{{\rm{wsct}}}^{m + 1} = \sum\limits_{i = 1}^{{N_{\rm F}}} {J_{{\rm{w}}i}^{m + 1}{{\left( {{p_{{\rm{bh}}}} - {p_{{\rm{ave}}i}}} \right)}^{m + 1}}} \end{array}\right. $

(26)

$ J_{{\rm{g}}i} = \dfrac{{1.086{K_{{\rm{F}}i}}{K_{{\rm{rg}}i}}{w_{{\rm{F}}i}}}}{{{\mu _{\rm{g}}}{B_{\rm{g}}}\left( {\ln {\dfrac{{{r_{0i}}}}{{{r_{\rm w}}}}} + {S_{{\rm{c}}i}}} \right)\left( {1 + \sqrt {1 + 4b{q_{d0}}} } \right)}} $

(27)

$ J_{{\rm{w}}i} = \dfrac{{0.543{K_{{\rm{F}}i}}{K_{{\rm{rw}}i}}{w_{{\rm{F}}i}}}}{{{\mu _{\rm{w}}}{B_{\rm{w}}}\left( {\ln {\dfrac{{{r_{0i}}}}{{{r_{\rm{w}}}}}} + {S_{{\rm{c}}i}}} \right)}} $

(28)

同理，由第$k+1$迭代步的值来逼近$m+1$时间步下的值，针对定井底流压生产和定产生产等方式，展开产量项并代入式(24)中，构建完整的求解矩阵。对于每一时间步，通过循环迭代，直到$\delta \boldsymbol{{X}}$满足精度。并由式(21)和式(22)获得$m+1$时间步下的稳定压力值和饱和度值，并开始下一时间步下的循环。

4 致密气藏多级压裂水平井生产动态分析

针对国内X致密气藏某区块中的一口多级压裂水平井，开展生产动态数值模拟分析。计算区域体积为1 600 m$\times$600 m$\times$50 m，水平井段长度1 200 m，通过加砂水力压裂改造形成9条垂直于水平井筒的主裂缝。

采用Comsol多场耦合方法对计算区域进行网格剖分，共形成了1 417个非结构四面体单元(图 2)。基质系统和裂缝系统的相对渗透率曲线和毛管力曲线如图 3所示。多级压裂水平井以15 MPa的井底流压定压生产。模型基本参数见表 1。

4.1 井型对比

利用数值模拟模型，对比了在致密气藏相同计算区域内，分别采用传统的非压裂水平井和多级压裂水平井进行开采时的生产动态变化情况，结果如图 4所示，可以看出，经过相同的10 a生产时间，多级压裂水平井的压力降落范围更广，日产气量和累积产气量都明显高于传统非压裂水平井生产情况(图 5)。上述现象源于经过水力压裂改造，近井筒区域渗流特征得到显著改善，压裂裂缝形成的高渗流通道，有效沟通了远端原始致密储层。因此，对于致密气藏，采用多级压裂水平井进行开发的生产效益更高。

4.2 含水饱和度的影响

致密气藏受周围水体分布和原生水的作用而出现的气水两相流动对压裂水平井生产动态影响显著。本文假设了3组含水饱和度值，利用数值模拟模型研究饱和度对产量曲线的影响(图 6，图 7)，含水饱和度越小，产气量越大，产水量越小。

具体来说，含水饱和度越小，多级压裂水平井的高产、稳产周期越长，以定井底流压开采的日产量越高。这是因为较小的含水饱和度意味着较大的气体饱和度，也代表了较大的流体流动性和储层条件下较大的气体体积。

4.3 应力敏感效应的影响

随着多级压裂水平井开井生产后压力下降，原始致密储层中无支撑剂作用的天然裂缝系统的渗透率呈现下降趋势。数值模拟计算表明，随着应力敏感系数的增加，日产气量和累产气量均下降(图 8)。根据式(2)，当压力下降时，较高的应力敏感性导致裂缝渗透率降低越多。当应力敏感性系数$d_{\rm f}$=0.1 MPa$^{-1}$时，开井生产10 a后的累产气量为3.07$\times$10$^8$ m$^3$，比不考虑应力敏感时减少约0.2$\times$10$^8$ m$^3$。

4.4 压裂裂缝空间分布的影响

对致密气藏进行水力压裂改造时，在储层物性非均质性，地应力分布特征等因素的作用下，压裂裂缝在三维空间的延伸扩展中可能出现不沿着垂直井筒方向扩展，从而形成倾斜裂缝的现象(图 9和图 10)。本文假设了4种压裂裂缝空间非对称分布情况并运用本文的模拟方法进行生产动态预测。通过与对称、垂直压裂裂缝模型的模拟计算结果对比分析可知，针对致密气藏，通过水力压裂改造尽可能形成三维空间垂直于水平井筒的压裂主缝，所获得的产气量最高。这是因为垂直压裂裂缝能够获得最大的纵向有效缝高，并保证气体经压裂裂缝流入井筒的流量最大。

同时，致密气藏广泛发育的天然裂缝也将极大地影响压裂裂缝的延伸规律。在压裂施工中，水力压裂裂缝与天然裂缝相互作用主要存在压裂裂缝沿天然裂缝生长，穿透天然裂缝沿原路径扩展，在相交处分叉、转向和发生剪切破坏等。可进一步构建如图 11所示分支缝网压裂水平井模型，并采用三维四面体网格对分支缝网进行空间离散。通过与不存在分支缝网的传统单一主缝压裂水平井模型日产气量和累积产气量进行对比分析(图 12)，可以看出，分支缝网模型的日产气量和累产气量均高于单一主缝模型，表明对于致密储层，在原始天然裂缝分布影响下，通过压裂技术手段形成的分支缝网系统可以有效地改善近井地带渗流特征，扩大压裂改造范围，从而提高致密气藏多级压裂水平井产量。

4.5 压开程度的影响

对致密气藏中水平井进行水力压裂施工时，存在部分拟压裂段(簇)难以被压开或压开后裂缝扩展范围低于压前预期设计要求，即存在部分段(簇)压裂改造不充分的情形。利用本文所构建的“降维”离散裂缝模型以及非结构三维四面体网格，可以简便快捷地模拟存在垂向上部分压开裂缝的多级压裂水平井生产动态特征，结果如图 13所示。可以看出，部分压开时形成裂缝对产量曲线影响显著，同完全压开时形成裂缝的井相比，部分压开裂缝井的日产量和累产气量均较低。在开井生产10 a之后，完全压开裂缝的水平井累产气量为2.50$\times$10$^8$ m$^3$，比部分压开裂缝井多出约2.81$\times$10$^7$ m$^3$的产量。因此，在致密气藏压裂改造中，需要结合地质参数、地应力和天然裂缝分布特征，尽可能实现压裂完全压开目标储层，以获取最大产能。

5 结论

（1）基于双重连续介质-离散裂缝耦合模型的多级压裂水平井数值模拟能较好地刻画致密储层和水力压裂裂缝系统流动特征，有效预测致密气藏多级压裂水平井生产动态。

（2）致密储层含水饱和度越低，气井前期产量越高，稳产期越长；应力敏感系数越大，裂缝渗透率降低越多，气井累产气量降低明显。

（3）在压裂改造过程中，垂直裂缝与倾斜裂缝相比，纵向上有效缝高越大，压后日产气量越高；分支缝网系统可以有效改善近井地带渗流特征，扩大压裂改造范围，提高气井产量；提高储层的压开程度能显著提高单井产量。

符号说明

$v$—天然裂缝中的渗流速度，m/s；

$\mu _{\rm g}$—气体黏度，mPa$\cdot$s；

$B_{\rm g}$—气体体积系数，无因次；

$\nabla p$—压力梯度，MPa/m；

$K_{\rm{fi}}$—页岩储层天然裂缝原始渗透率，mD；

$d_{\rm f}$—应力敏感系数，MPa$^{-1}$；

$p_{\rm{fi}}$—天然裂缝系统原始压力，MPa；

$p_{\rm f}$—天然裂缝系统压力，MPa；

$K_{\rm{fe}}$—考虑应力敏感效应的天然裂缝修正渗透率，mD；

$J$—单毛细管中耦合黏性流和分子扩散的摩尔质量通量，mol/s；

$\lambda$—孔隙毛细管截面直径，m；

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

$M$—摩尔质量，kg/mol；

$\Delta p$—生产压差，MPa；

$L_{\rm t}$—真实毛细管长度，m；

$K{ n}$—Knudsen数，无因次；

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

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

$T$—温度，K；

$p$—气体压力，MPa；

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

$Q_{\rm{T}}$—弯曲毛细管总摩尔质量流量，mol/s；

$\lambda_{\min}$—最小毛细管直径，m；

$\lambda_{\max}$—最大毛细管直径，m；

$N$—弯曲毛细管总数；

$K_{{\rm{am}}}$—基质系统表观渗透率，mD；

$D_{\rm T}$—迂曲分形维度，无因次；

$L_0$—岩样直线边长，m；

$D_{\rm f}$—毛管弯曲形状分形维度，无因次；

$q_{\rm{w, fm}}$—基质与天然裂缝系统水相窜流量，1/s；

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

$K_{\rm{rwm}}$—基质系统水相相对渗透率，无因次；

$\mu_{\rm{wm}}$—基质系统水相黏度，mPa$\cdot$s；

$B_{\rm{wm}}$—基质系统水相体积系数，无因次；

$p_{\rm{wm}}$—基质系统水相压力，MPa；

$p_{\rm{wf}}$—天然裂缝系统水相压力，MPa；

$q_{\rm{g, fm}}$—基质与天然裂缝系统气相窜流量，1/s；

$\alpha$—窜流形状因子，m$^{-2}$；

$K_{\rm{rgm}}$—基质系统气相相对渗透率，无因次；

$\mu_{\rm{gm}}$—基质系统气相黏度，mPa$\cdot$s；

$B_{\rm{gm}}$—基质系统气相体积系数，无因次；

$p_{\rm{gm}}$—基质系统气相压力，MPa；

$p_{\rm{gf}}$—天然裂缝系统气相压力，MPa；

$\nabla {p_{{\rm{wm}}}}$—基质系统水相压力梯度，MPa/m；

$S_{\rm{wm}}$—基质系统含水饱和度，无因次；

$\nabla {p_{{\rm{gm}}}}$—基质系统气相压力梯度，MPa/m；

$\phi_{\rm{m}}$—基质系统平均孔隙度，无因次；

$S_{\rm{gm}}$—基质系统含气饱和度，无因次；

$ K _{\rm{rgf}}$—天然裂缝系统气相相对渗透率，无因次；

$\mu _{\rm{gf}}$—天然裂缝系统气相黏度，mPa$\cdot$s；

$B _{\rm{gf}}$—天然裂缝系统气相体积系数，无因次；

$\nabla {p _{{\rm{gf}}}}$—天然裂缝系统气相压力梯度，MPa/m；

${q_{{\rm{g, Ff}}}}$—天然裂缝与压裂裂缝气相窜流量，1/s；

$\phi_{{\rm{f}}}$—天然裂缝系统平均孔隙度，无因次；

$S_{\rm{gf}}$—天然裂缝系统含气饱和度，无因次；

$\nabla {p_{{\rm{gF}}}}$—压裂裂缝气相压力梯度，MPa/m；

$\mu _{\rm{gF}}$—压裂裂缝气相黏度，mPa$\cdot$s；

$K _{\rm{F}}$—压裂裂缝平均渗透率，mD；

$K _{\rm{rgF}}$—压裂裂缝气相相对渗透率，无因次；

$\beta$—气水两相紊流系数，1/m；

$\rho _{\rm{gF}}$—压裂裂缝气体密度，kg/m$^3$；

$v_n$—高速非达西渗流速度，m/s；

$v_{{\rm d}n}$—达西渗流速度，m/s；

$K _{\rm{rwf}}$—天然裂缝系统水相相对渗透，无因次；

$\mu _{\rm{wf}}$—天然裂缝系统水相黏度，mPa$\cdot$s；

$B _{\rm{wf}}$—天然裂缝系统水相体积系数，无因次；

$\nabla {p _{{\rm{wf}}}}$—天然裂缝系统水相压力梯度，MPa/m；

$S_{\rm{wf}}$—天然裂缝系统含水饱和度，无因次；

$B _{\rm{gF}}$—压裂裂缝气体体积系数，无因次；

$\nabla {p _{{\rm{gF}}}}$—压裂裂缝气相压力梯度，MPa/m；

$q_{\rm{gsct}}$—地面条件下的产气量，m$^3$/d；

$\phi_{{\rm{F}}}$—压裂裂缝平均孔隙度，无因次；

$S_{\rm{gF}}$—压裂裂缝含气饱和度，无因次；

$ K _{\rm{rwF}}$—压裂裂缝水相相对渗透率，无因次；

$\mu _{\rm{wF}}$—压裂裂缝水相黏度，mPa$\cdot$s；

$B _{\rm{wF}}$—压裂裂缝水相体积系数，无因次；

$\nabla {p _{{\rm{wF}}}}$—压裂裂缝水相压力梯度，MPa/m；

${q_{{\rm{w, Ff}}}}$—天然裂缝与压裂裂缝水相窜流量，1/s；

$q_{\rm{wsct}}$—地面条件下的产水量，m$^3$/d；

$S_{\rm{wF}}$—压裂裂缝含水饱和度，无因次；

$\varOmega _{\rm{d}}$—全储层渗流空间，m$^3$；

$\varOmega _{\rm{f}}$—天然裂缝系统渗流空间，m$^3$；

$\varOmega _{\rm{F}}$—压裂裂缝渗流空间，m$^3$；

$F$—综合渗流偏微分方程组；

$w_{\rm F}$—压裂裂缝宽度，m；

$p_{\rm m}$—基质系统压力，MPa；

$p_{\rm{mi}}$—基质系统初始压力，MPa；

$v=1, 2, 3, 4$—四面体顶点；

$p_v$—任意顶点压力，MPa；

$p_{{\rm w}v}$—任意顶点含水饱和度；

$V$—四面体体积，m$^3$；

$m$—时间步；

$k$—迭代步；

$r_{0i}$—第$i$井源汇点的等效井半径，m；

$\boldsymbol{{k}}$—系数矩阵；

$δ\boldsymbol{{X}}$—求解变量矩阵；

$\boldsymbol{{R}}$—余项矩阵；

$A_i$—水平井筒与裂缝面的交点所构成的控制体积的底面积，m$^2$；

$J_{{\rm g}i}$—考虑高速非达西效应的采气指数，m$^3$/(d$\cdot$MPa)；

$J_{{\rm{w}}i}$—产水指数，m$^3$/(d$\cdot$MPa)；

$p_{\rm{bh}}$—井底流压，MPa；

$p_{{\rm{ave}}i}$—井源汇点所在网格块的平均压力，MPa；

$N_{\rm F}$—压裂裂缝与水平井筒的总交点个数；

$S_{{\rm{c}}i}$—第$i$井源汇点的表皮系数，无因次。