考虑应力敏感和复杂运移的页岩气藏压力动态分析

姜瑞忠，张福蕾，崔永正，潘 红，张 旭，张春光，沈泽阳

（1.中国石油大学（华东）石油工程学院，山东青岛266580；2.中国石油大港油田分公司采油工艺研究院，天津300280）

0 引言

页岩气藏作为重要的非常规油气资源之一，近年来已成为国内外学者研究的热点［1-2］。页岩气存在吸附解吸、基质中扩散、基质向裂缝间窜流和裂缝中渗流等多种气体运移机制［3-4］。页岩储层渗透率极低，需对其进行水力压裂改造［5-7］。压裂改造可形成高导流区，把包括裂缝在内的这一区域称为改造储层体积（SRV）［8-9］，其形态可简化为一椭圆形区域［10］。一些学者［11-12］利用复合气藏模型分析了考虑椭圆SRV的压裂井压力变化特征，但都未考虑基质向裂缝间的窜流。Zhang等［11］建立了考虑椭圆SRV的双重介质页岩气藏压裂井复合模型，依照Fick第一扩散定律在基质中拟稳态扩散。姜瑞忠等［12］提出依照Fick第二扩散定律的压裂改造复合页岩气藏压力分析方法，运用椭圆坐标系模型，综合分析了非稳态扩散、解吸、渗流等传质机理。页岩气藏渗透率低，储层表现出很强的应力敏感性，其对压力的影响不容忽视［13］。尹洪军等［14］研究了均质页岩气藏垂直裂缝井椭圆流压力动态问题，模型考虑了应力敏感效应的影响，但也未涉及SRV的分析。

综上所述，目前还未有考虑应力敏感性、基质向裂缝间非稳态窜流、椭圆SRV等的数学模型。本次研究根据页岩气吸附解吸、Knudsen扩散、非稳态窜流和渗流等多种气体运移机制对压力的影响，来建立综合应力敏感性和椭圆SRV的复合页岩气藏压裂井试井模型，并利用Mathieu函数、Pedrosa变量代换、正则摄动理论等方法对数学模型进行求解，以期提出的压力动态分析方法可为水力压裂开采页岩气藏提供一定的理论支持。

1 物理模型

考虑SRV的压裂改造页岩气藏模型如图1所示，内外两区均为椭圆形，外区为原始储层，内区为压裂改造区（SRV），其孔隙度和渗透率较外区高。模型假设如下：①外区边界无限大，内区半径为ξ1；②气藏储层厚度为h，初始地层压力为Pi；③基质假设为球状，半径为Rm，基质中存在吸附气和自由气，裂缝中为自由气；④压裂井位于气藏中心，以定产量qsc生产；⑤考虑基质中气体的吸附解吸，气体扩散满足Knudsen定律，并且考虑从基质到裂缝间的非稳态窜流；⑥只考虑裂缝受应力敏感的影响，由于应力敏感主要发生在近井地带，只在内区考虑其对渗流的影响［15］；⑦忽略重力和毛管压力的影响。

图1 考虑椭圆SRV的复合页岩气藏压裂井模型Fig.1 Fractured well model in composite shale gas reservoir considering elliptical SRV

2 数学模型建立

页岩气按照Knudsen扩散模型在纳米级孔隙中扩散［16］，由于浓度差作用的扩散所产生的质量通量为

式中：Jk为由于浓度差作用的扩散所产生的质量通量，kg/（m2·s）；D为扩散系数，m2/s；Mg为气体摩尔质量，kg/mol；Cm为基质中气体浓度，mol/m3；rm为基质中径向距离，m。

气体状态方程为

式中：pm为基质压力，Pa；V 为气体体积，m3；Z 为气体压缩因子；n为气体物质的量，mol；R为通用气体常数，Pa·m3/（mol·K）；T为气藏储层温度，K。

由式（2）可求得气体浓度为

式中：cgm为基质中气体压缩系数，Pa-1；ρ为气体密度，kg/m3。

将式（3）代入式（1）并联立式（4）整理可得

式中：Jp为压力差作用的渗流所产生的质量通量，kg/（m2·s）；km为基质渗透率，mD；μ为气体黏度，mPa∙s。

式（2）变形可得

则浓度差和压力差共同作用所产生的质量通量为

式中：ka为基质的表观渗透率，mD。

球状基质表面处的质量流速［17］为

式中：v为气体质量流速，kg/（m2·s）；Rm为基质块半径，m。

单位体积基质在单位时间内向裂缝的非稳态窜流量为qm，基质球面处的流速等于单位时间基质的流出量与基质表面积的比值，即

式中：qm为气体窜流质量流量，kg/（m3·s）。

将式（11）代入式（10）整理可得

式中：qd为气体解吸质量流速，kg/（m3·s）；φf为裂缝孔隙度，%；φm为基质孔隙度，%；GL为兰格缪尔体积，m3/m3；pL为兰格缪尔压力，Pa；t为时间，s。

引入气体拟压力表达式

式中：ψ为拟压力，Pa/（mPa·s）；

内区裂缝考虑应力敏感性，渗透率模量定义如下

式中：γ 为渗透率模量，mPa·s/Pa2；kf1为内区裂缝渗透率，mD；ψf1为内区裂缝拟压力，Pa2/（mPa·s）。式（15）可变形为

求解式（16）可得

式中：式中：φf1为内区裂缝孔隙度，%；ctf1为内区裂缝总压缩系数，Pa-1；Rm1为内区基质块半径，m；ka1为内区基质的表观渗透率，mD；ψm1为内区基质拟压力，Pa2/（mPa·s）；rm1为内区基质中径向距离，m。

内区基质的渗流方程为

式中：φm1为内区基质孔隙度，%；ctm1为内区基质总压缩系数，Pa-1；ka1为内区基质的表观渗透率，mD；cd1为内区考虑解吸的附加压缩系数，Pa-1，表达式如下

式中：cd为考虑解吸的附加压缩系数，Pa-1；psc为标准条件下的压力，Pa；Tsc为标准条件下的温度，K；ψL为兰格缪尔拟压力，Pa2/（mPa·s）。

式中：L为参考长度，m；qsc为恒定地面气体产量，m3/s；C为井筒存储系数，m3/Pa；h为气藏厚度，m；β为内外区储容比；θ为渗透率比；φm2为外区基质孔隙度，%；ctm2为外区基质总压缩系数，Pa-1；φf2为外区裂缝孔隙度，%；ctf2为外区裂缝总压缩系数，Pa-1；ωf为裂缝储容比；ωm为裂缝储容比；ωd为考虑解吸的储容比；λ为窜流系数；M为流度比；kf2为外区裂缝渗透率，mD。

将无因次参数代入式（18）—（19）中，得无因次条件下裂缝和基质的渗流方程分别为

式中：ωf1为内区裂缝储容比；λ1为内区窜流系数；θ1为内区渗透率比；ωm1为内区基质储容比；ωd1为内区考虑解吸的储容比。

利用Pedrosa变量代换处理裂缝渗流方程的强非线性［19］，表达式如下

式中：χ1D为摄动变换函数。

将式（23）代入式（21）可得

将式（24）进行整理，等式两边约去相同项可得

无因次渗透率模量通常很小（γD＜＜1），0阶近似解满足工程精度要求，式（25）变为

式中：s为拉普拉斯变量。

球形基质的初始及内外边界条件经无因次化、Pedrosa变量代换及拉普拉斯变换后为

将式（31）—（33）与式（30）联立可得

将式（34）代入式（29）可得

笛卡尔坐标转换为椭圆坐标的方法为

式中：ξ，η为椭圆坐标。

将式（36）代入式（35）可得椭圆坐标系的渗流方程为

外区无因次化及拉普拉斯变换后的裂缝和基质渗流方程分别为

式中：ωf2为外区裂缝储容比；λ2为外区窜流系数；θ2为外区渗透率比；ωm2为外区基质储容比；ωd2为外区考虑解吸的储容比。

利用与内区相同的求解方法对式（38）、（39）进行求解，并利用式（36）将坐标转换为椭圆坐标为

式中：ξw为椭圆坐标下井筒半径。

无限大外边界条件为

ξ=ξ1时，内外区压力和流量相等，交界面条件为

3 数学模型求解

根据式（42）和 Mathieu 函数性质［20］得式（37），（40）的解为

式中：ce2n为2 n阶第一类Mathieu函数角函数；Ge2n为2 n阶第一类Mathieu函数径函数；Ue2n为2 n阶第二类Mathieu函数径函数。

由式（41），（43）和（44）得

其中A2n0为Mathieu函数展开式的系数。

联立式（47）—（49）可得 D2n，I2n和 H2n，Pedrosa变量代换后的产量为

由角度Mathieu函数的性质［20］，有下式

考虑井筒存储效应和表皮系数时，式（53）变为

式中：S为表皮系数。

利用Stehfest数值反演法［22］可得真实空间下的压力χwD，则考虑应力敏感的井底拟压力为

4 典型曲线分析

图2为考虑椭圆SRV的复合页岩气藏拟压力及导数典型曲线，基本参数取值为：CD=10，S=0.1，ξ=11，kf1=0.02 mD，kf2=0.005 mD，D=10-6（m2/s），cd=0.013 MPa-1。依据曲线形态可将流动阶段划分为9段：①井筒存储阶段，拟压力及导数曲线重合，斜率为1；②表皮因子影响阶段，导数曲线形状为一驼峰；③内区裂缝径向流阶段，表现为值为0.5的水平线；④内区窜流阶段，导数曲线为一凹子；⑤内区整体径向流阶段，导数曲线水平；⑥过渡流阶段，压力传播到内外区交界面处；⑦外区裂缝径向流阶段；⑧外区窜流阶段；⑨储层整体径向流阶段。

5 敏感性分析

对渗透率模量、SRV半径、外区裂缝渗透率、扩散系数、解吸压缩系数等参数进行敏感性分析，基本参数取值为：CD=10，S=0.1，γD=0.01，ξ=11，kf1=0.02 mD，kf2=0.005 mD，D=10-6（m2/s），cd=0.013 MPa-1。

图2 拟压力及导数典型曲线Fig.2 Typical curves of pseudo-pressure and derivative

5.1 渗透率模量影响

图2反映了渗透率模量γD对复合页岩气藏压力的影响，渗透率模量的增加使得拟压力及导数曲线上升，并且外区曲线所受到的影响更加强烈。在整体径向流阶段，导数曲线开始上翘，这是应力敏感效应产生的典型曲线特征。当气藏考虑应力敏感性时，储层渗透率随压力的减小而降低，页岩气的流动阻力增大，导致流动所需压差增大。

5.2 SRV半径影响

图3反映了SRV半径ξ1对复合页岩气藏压力的影响，SRV半径对内区整体径向流阶段和过渡流阶段产生的影响较大。随SRV半径的增大，拟压力曲线向下移动，导数曲线向右移动，当SRV半径足够小时，内区整体径向流阶段消失。因为SRV半径越大，储层改造的体积越大，越有利于气井生产的进行，气井生产所需的压差减小。

5.3 外区裂缝渗透率影响

图4反映了外区裂缝渗透率kf2对复合页岩气藏压力的影响，外区裂缝渗透率影响了内区整体径向流阶段之后的曲线形态。外区裂缝渗透率的增加与流度比的减小相对应，使得外区拟压力及导数曲线下移。流度比减小，则外区未改造的原始储层条件相对内区变好，外区流动所需的压降减小，则压裂改造提高裂缝渗透率可以很好地提高开发效果。

图3 SRV半径对复合页岩气藏压力的影响Fig.3 Influence of SRV radius on pressure of composite shale gas reservoir

图4 外区裂缝渗透率对复合页岩气藏压力的影响Fig.4 Influence of fracture permeability in outer zone on pressure of composite shale gas reservoir

5.4 扩散系数影响

图5反映了扩散系数D对复合页岩气藏压力的影响，扩散系数主要影响内外区窜流阶段。扩散系数的增加导致了页岩基质表观渗透率的增加，使得窜流系数增大，窜流阶段提前，在导数曲线上表现为内外区凹子向左移动，拟压力曲线向下移动，气井生产所需的压差减小。

5.5 解吸压缩系数影响

图6反映了解吸压缩系数cd对复合页岩气藏压力的影响，解吸压缩系数影响内区裂缝径向流之后的阶段。解吸压缩系数增大使得考虑解吸的储容比ωd增大，窜流阶段的凹子变深且变宽，影响了导数曲线窜流阶段前后径向流的持续时间，解吸压缩系数增大导致过渡流阶段右移，驼峰升高，气井生产所需压差降低。

图5 扩散系数对复合页岩气藏压力的影响Fig.5 Influence of diffusion coefficient on pressure of composite shale gas reservoir

图6 解吸压缩系数对复合页岩气藏压力的影响Fig.6 Influence of desorption compressibility on pressure of composite shale gas reservoir

6 结论

（1）在综合考虑应力敏感效应和椭圆SRV的基础上，建立了压裂改造复合页岩气藏压力动态分析模型，并考虑了页岩气吸附解吸、Knudsen扩散、非稳态窜流和渗流等多种气体运移机制对压力的影响。

（2）采用Mathieu函数、Pedrosa变量代换、正则摄动理论、拉普拉斯变换和Stehfest数值反演法等多种方法相结合对数学模型进行求解，并绘制压力动态典型曲线，将其划分为9个流动阶段。

（3）渗透率模量的增加使得无因次曲线上移，气井生产所需压差增大；SRV半径越大，储层改造的体积越大，气井生产所需压差减小；外区裂缝渗透率的增加与流度比的减小相对应，流度比越小，外区流动所需的压降减小，显示出压裂改造对页岩气藏开发产生的良好效果；扩散系数的增加导致了页岩基质表观渗透率的增加，使得窜流系数增大；解吸压缩系数增大使得考虑解吸的储容比增大，窜流段的凹子变深且变宽。