根据国际能源机构（IEA）的统计资料显示，世界上煤层气储量约为260×10¹³ m³，其中，中国的煤层气储量超过10×10¹³ m^3[1-2]。煤层气作为一种优质的清洁能源，是国家大力提倡开发的新能源之一，其开发利用对缓解中国的能源紧张形势具有重大的战略意义。由于煤层气藏储层的低孔、低渗透特性，常规直井开采效果不好，而多分支水平井能够增加井筒与煤层气藏储层的接触面积，增加各个分支的波及面积和泄气面积，从而大大提高单井产量和产气速度^[3-5]，又由于钻井技术的飞速发展和成本的降低，多分支水平井技术在煤层气藏中的应用较为广泛。

诸多中外学者依据不同井型对煤层气藏的生产动态开展研究。张先敏等基于分形理论，考虑了低渗透煤层气藏启动压力梯度的影响，建立了分形煤层气藏直井渗流数学模型，分析了分形维数等对生产动态的影响^[6]。宋洪庆等基于渗流力学基本理论，推导了耦合煤层气藏解吸效应和低速非达西渗流特征的压力和产能特征方程，求解给出了径向流解析解^[7]。罗毅等采用裂缝波及区和未波及区等复合两区模型，推导了煤层气藏压裂水平井的产能公式，分析了水平长度、储层厚度以及缝网波及宽度等对产能的影响^[8]。ZHAO 等建立了煤层气藏无限导流能力压力直井的渗流数学模型，同样将整个煤层气藏分为复合两区，包括受压裂井影响形成的压裂改造的内区和未压裂改造的外区，求解绘制了产量递减动态曲线图，同时对比了拟稳态和非稳态的影响^[9]。ZHANG 等建立了复合两区煤层气藏水平井渗流数学模型，通过求解井底压力和产量，分别绘制了井底压力动态曲线和产量动态曲线图，进行了相关的敏感性参数分析^[10]。

虽然在煤层气藏开采中多分支水平井技术应用广泛，但是相关煤层气藏多分支水平井非稳态产能模型的研究较少^[11]。同时，大量研究发现，煤层气藏储层渗透率往往存在着较强的应力敏感特性^[12-16]，在进行煤层气藏生产动态研究中，应考虑到该特性。为此，笔者建立了煤层气藏多分支水平井非稳态产能模型，耦合了割理系统的渗透率应力敏感特性和基质系统的非稳态扩散，求解绘制出产量递减动态曲线，并分析了相关因素的影响规律。

1 模型建立

1.1 物理模型

假设在无限大顶底边界封闭的地层中心，存在着以定产量生产的一口多分支水平井（图1），该多分支水平井由 1 个水平井筒和 N 个分支井筒组成，每个分支井筒的长度为 L _a，各个分支井筒之间的距离为 d，各个分支井筒与水平井筒之间的夹角为 χ_a。考虑煤层气藏的吸附解吸和扩散特性，分别用 LANGMUIR 等温吸附理论^[17] 和 FICK 第二扩散定律进行表征，气体流动途径如图2所示。同时，忽略重力、毛细管力、黏滞力以及温度变化，假设分支井筒和水平井筒导流能力无限大。

1.2 数学模型

为了方便模型推导和求解，在此引入一系列无因次变量（表1），从而建立点源数学模型进行求解。

1.2.1 割理系统

结合质量守恒方程、状态方程、运动方程以及 PEDROSA^[18] 提出的表征渗透应力敏感特性的渗透率模量表达式，推导得到割理系统的流动方程为：

结合拟压力和按拟压力定义的渗透率模量，按照表1 中的无因次变量进行无因次化之后，得到割理系统的无因次化流动方程为：

初始条件为：

外边界条件为：

顶底封闭界面条件为：

内边界定产条件为：

1.2.2 基质系统

采用 FICK 第二定律来描述基质系统的煤层气非稳态扩散现象，假设基质块为球形，球形半径为 R，径向距离为r _m （图3）^[19-20]。

结合表1 中的无因次变量，基质系统中的非稳态扩散方程可以表示为：

初始条件为：

外边界条件为：

球形对称内边界条件为：

2 数学模型求解

对于数学模型的求解，整体思路为：先求解出基质系统的气体体积浓度表达式，然后代入割理系统的流动方程中，通过一系列数学方法进行求解。

2.1 基质系统求解

将基质系统扩散方程进行Laplace变换后，可以得到：

为了求解方便，在此使用变量替换法，令：

则（12）式变为：

（14）式为二阶常系数微分方程组，可以求解得到其通解为：

因此：

从而：

2.2 割理系统求解

由于割理系统流动方程的非线性较强，在此采用Pedrosa变换和摄动变换消除非线性，其表达式分别为：

通过 Pedrosa 变换和摄动变换后，进行 Laplace 变换，可以得到：

将（17）式代入（22）式，可得：

其中：

对（23）式进行有限余弦变换，则：

其中：

（25）式的通解为：

由外边界条件可得A =0，则（27）式可写为：

结合外边界条件可以求得：

进而点源函数解为：

结合有限余弦反演得到 Laplace 空间下的点源函数解为：

为了求解得到井底压力，在此利用数值离散和叠加原理的方法，将每个分支井筒划分为 M个微元段，同时微元段线流量密度为q _a，b。这样第a个分支的第b个微元段在任意一点形成的压力均可写为：

这样，所有微元段在第 α 个分支第 β 个微元段所形成的压力为井底压力（因为假设分支井筒和水平井筒导流能力无限大）：

另外，因各微元段流量总和为水平井筒总流量，有：

（33）式和（34）式可以构成一个矩阵方程组，求解可得井底压力，进而结合杜哈美原理可得耦合表皮系数和无因次井筒储集系数的井底压力解为：

然后，通过Stehfest数值反演^[21] 和（35）式可以得到真实空间下的井底压力解为：

依据 VAN EVERDINGEN 等的研究成果^[22]，可以得到Laplace空间下产量的表达式为：

3 结果与讨论

3.1 模型验证

将本文模型中的单相气换为单相油，将模型中的f（u）替换为文献^[23]中的f（u），并与其中的模型进行对比（图4）。基本参数取值为：S =0.01，C _D=1× 10^-6，L _aD=1，d _D=10，N =4，ω=0.1，λ=1×10^-4 ，χ_a=60°， γ_pD=0。从图4中可以看出，本文模型计算结果与文献^[23]模型结果吻合较好，从而验证了本文模型的可靠性。

3.2 典型曲线分析

煤层气藏多分支水平井拟压力曲线和产能曲线（图5）的基本参数取值为：S =0.2，C _D=5×10^-6，L _aD= 1，d _D=8，N =4，ω=0.1，ω_ad=5，λ=1×10^-9，χ_a=60°，γ_pD=0。根据曲线特征，将其划分为9个流动阶段：①井筒储集阶段。拟压力和拟压力导数曲线重合，产能和产能导数曲线重合。②表皮效应阶段。拟压力导数曲线呈驼峰状，产能和产能导数曲线继续下降。③ 分支井筒第一径向流阶段。拟压力导数呈水平线，产能导数曲线呈下降趋势。④分支井筒第一线性流阶段，拟压力导数曲线斜率为 0.5，产能导数曲线呈上升趋势。⑤分支井筒的第二径向流阶段。拟压力导数呈水平线，产能与产能导数曲线平行。⑥多分支水平井的整体第二线性流阶段。拟压力导数曲线斜率为0.5，产能导数呈水平线。⑦割理系统拟径向流阶段。拟压力导数呈水平线，产能与产能导数曲线平行。⑧基质系统向割理系统的窜流阶段。拟压力导数和产能导数曲线均呈现出一个宽的凹子。⑨整体径向流阶段，此时整个系统的流动达到平衡状态。

从考虑非稳态扩散和拟稳态扩散的产能曲线对比（图6）中可知：前期，两者的产能曲线重合；而后期，拟稳态扩散呈现出的凹子窄且深，而非稳态扩散呈现的凹子宽且浅，采用拟稳态扩散模型时，窜流强度较高，所以就计算产能而言精确度不高。

3.3 敏感性参数分析

无因次渗透率模量  从图7 中可以看出，无因次渗透率模量主要影响产能曲线最后部分。其原因为，在生产初期井底压力和地层原始压力相差较小，割理系统渗透率的应力敏感性较差。随着时间的推移，应力敏感性逐渐增加，割理系统渗透率下降，气体流动阻力增强，从而对产能造成影响。无因次渗透率模量取值越小，后期产能较高；相反，若取值越大，则产能较低。

分支井筒间距  从图8 中可以看出，分支井筒间距造成的影响在产能曲线后部分，随着分支井筒间距的减小，产能和产能导数曲线下降。其原因为，在生产前期，各个分支井筒之间的干扰并未形成；随着生产的进行，分支井筒之间的干扰开始形成，同时分支井筒间距越小，干扰发生的时间越早，干扰程度越大，从而较小的分支井筒间距会造成产能下降的时间提前且下降的速度变快。

分支井筒个数  分支井筒个数反映多分支水平井的生产规模，个数越多，流入水平井筒的气体越多，从而生产能力越强。从图9中可以看出，随着分支井个数的增加，在产能曲线中间段部分，产能越大，且产能降低速度越小，在中期整体表现为产能和产能导数曲线的上升，但在后期可以看出增加分支井筒个数的增产效果不明显。

窜流系数  从窜流系数对产能曲线的影响（图10）可以看出，窜流系数越小，说明基质系统和割理系统之间的物理性质差别越大，基质系统解吸出来的气体越难扩散窜流进入割理系统中被产出。在产能曲线上表现为窜流系数越小，窜流开始的时间越晚，凹子右移且下移。

吸附系数  从图11 中可以看出，吸附系数越大，窜流阶段的凹子就越深且越宽。这是因为吸附系数代表了基质系统以吸附形式存在的气体体积。因此，吸附系数越大，解吸出来的气体扩散窜流至割理系统的量越多，从而产能曲线的位置越高，窜流阶段的凹子越明显。

弹性储容比  从弹性储容比对产能曲线的影响（图12）可以看出，弹性储容比对产能曲线中期部分影响较大。弹性储容比越大，表明割理系统的储集能力越强，气体量越大，在生产压差不变的情况下，产出的气体越多。具体表现为窜流阶段的凹子越浅越窄，同时凹子出现时间推迟，产能曲线中部位置抬升越高。

4 结论

综合考虑割理系统的应力敏感特性及基质系统的非稳态扩散建立了煤层气藏多分支水平井非稳态产能数学模型，利用一系列数学方法求解得到产能表达式。通过模型对比验证了所提出模型的可靠性，并且绘制了典型的拟压力曲线以及产能曲线。

分别考虑拟稳态扩散和非稳态扩散的产能曲线对比，说明了非稳态扩散模型就产能计算而言精确度更高。无因次渗透率模量的增加会降低产能； 分支井筒间距和分支井筒个数主要影响产能曲线的中间位置，分支井筒间距越小产能曲线位置越低，分支井筒个数的增加会引起产能的增加；窜流系数和吸附系数主要影响窜流阶段；弹性储容比除了影响窜流阶段之外，其取值越大，产能越大。此外，基于本文模型还可以进行考虑煤层气藏非线性渗流的产能研究。

符号解释

L _a ——分支井筒的长度，m；d ——分支井筒间距，m； χ_a——分支井筒与水平井筒之间夹角，（°）；h ——储层厚度， m；x，y，z ——笛卡尔坐标，m；m ——拟压力，Pa²/ （Pa·s）；下标 D ——无因次；K _fhi ——水平方向割理系统初始渗透率，m²； T _sc——标准状态下的温度，K；p _sc——标准状态下的压力， Pa；Q ——多分支水平井产量，m³ /s；T ——温度，K；m_i ——初始条件下的拟压力，Pa² /（Pa·s）；γ_p——拟渗透率模量，mPa· s/Pa²；V _m——基质系统吸附气体体积浓度，m³ /m³；V _Ei ——标准状态下的气体平衡体积浓度，m³ /m³；C _m——基质系统下的气体体积浓度，m³ /m³；C _mi ——基质系统下的初始气体体积浓度，m³ /m³；t ——时间，s；Λ——无因次中间变量；L ——参考长度，m；q ——点源产量，m³ /s；λ——窜流系数；D ——扩散系数，m² /s；R ——基质单元半径，m；ω_ad——吸附系数；m _L—— Langmuir 拟压力，Pa²/ （Pa·s）；m _f ——割理系统拟压力，Pa² / （Pa·s）；ω——弹性储容比；ϕ_f ——割理系统下的孔隙度，f； C _tf ——割理系统压缩系数，Pa^-1；μ——黏度，Pa·s；r ——径向距离，m；κ——无穷小垂向距离，m；σ——无穷小径向距离， m；K _fvi ——垂直方向割理系统初始渗透率，m²；r _m——基质单元径向距离，m；γ——渗透率模量，Pa^-1 ；p _i ——储层初始压力，Pa；p _f ——割理系统压力，Pa；Z ——气体压缩因子，小数；C _t ——储层综合压缩系数，Pa^-1；z _wD——无因次点源纵向坐标；$\overline{C_{\mathrm{m}\mathrm{D}}}$——Laplace 空间下无因次基质系统气体体积浓度；u ——Laplace因子；$\overline{m_{\mathrm{f}\mathrm{D}}}$——Laplace空间下无因次割理系统拟压力；$\overline{V_{\mathrm{m}\mathrm{D}}}$——Laplace空间下无因次基质系统吸附气体体积浓度；$\overline{q_{\mathrm{D}}}$——Laplace空间下无因次点源产量；τ_fD——Pe⁃ drosa变换后的无因次割理系统拟压力；$\overline{{\tau }_{\mathrm{f}\mathrm{D}0}}$——Laplace空间下无因次割理系统拟压力；$\stackrel{~}{\overline{{\tau }_{\mathrm{f}\mathrm{D}0}}}$——有限余弦变换后的 La⁃ place 空间下无因次割理系统拟压力；g_n ——无因次中间变量；n ——取0，1，2，3，···；A ——中间系数；I ₀——第一类零阶虚宗量Bessel函数；B ——中间系数；K ₀——第二类零阶虚宗量Bessel函数；I ₁——第一类一阶虚宗量Bessel函数；K ₁—— 第二类一阶虚宗量 Bessel函数；M ——每个分支井筒离散微元段数；q _a，b——第 a 个分支第 b 个微元段的流量，m³ / s；$\overline{{\tau }_{{\mathrm{f}\mathrm{D}}_{a，b}}}$ ——Laplace 空间下第 a 个分支第 b 个微元段的压力响应；Γ_Da，b——第 a 个分支第 b 个微元段的无因次积分长度；$X_{{\mathrm{D}}_{\alpha ，\beta }}，Y_{{\mathrm{D}}_{\alpha ，\beta }}，Z_{{\mathrm{D}}_{\alpha ，\beta }}$ ——第 α 个分支第 β 个微元段无因次中点坐标；N——分支井筒个数；$\overline{q_{{\mathrm{D}}_{a，b}}}$——Laplace 空间下第 a个分支第 b 个微元段的线流量密度；C _D——无因次井筒储集系数； m _wD——真实空间下无因次井底压力；$\mathrm{\Delta }{l}_{{\mathrm{D}}_{a，b}}$——第 a 个分支第 b 个微元段的无因次长度；S ——表皮系数；$\overline{m_{\mathrm{w}\mathrm{D}}}$——La⁃ place空间下的无因次井底压力。

参考文献