朱维耀，亓倩，马千，邓佳，岳明，刘玉章

（1.北京科技大学土木与环境工程学院；2.中国石油勘探开发研究院）

页岩气不稳定渗流压力传播规律和数学模型

朱维耀1，亓倩1，马千1，邓佳1，岳明1，刘玉章2

（1.北京科技大学土木与环境工程学院；2.中国石油勘探开发研究院）

利用稳定状态依次替换法，研究了页岩基质储集层内压力扰动的传播规律，得到动边界随时间变化的关系，考虑解吸、扩散、滑移作用及动边界的影响，建立了页岩气不稳定渗流数学模型。采用拉普拉斯变换，求解了内边界定产、外边界为动边界条件下的不稳定渗流压力特征方程。结合中国南方某海相页岩气藏储集层参数，应用MATLAB编程，计算分析了页岩气不稳定渗流压力特征及其影响因素。研究表明：页岩气开采过程中，压力传播具有动边界效应，动边界随时间延续向外传播，且传播速度逐渐减慢；动边界使压力传播速度变慢，储集层压力下降减缓；页岩气解吸使压力传播速度减慢，地层压力下降减缓；扩散系数越大，地层压力下降越慢，且扩散系数影响逐渐减小。在气藏开采过程中，扩散、滑移对产气量贡献逐渐增加，占主要地位；渗流及解吸对产气量贡献逐渐减小后趋于平稳。图9参18

页岩气；不稳定渗流；压力传播规律；数学模型；动边界

0 引言

页岩气具有多尺度流动的特征，并以吸附、游离或溶解状态赋存于纳米—微米级页岩孔隙及裂缝中，气体产出机理主要为解吸、扩散、滑移等，为低速强非线性渗流[1-5]；PASCAL[6]和刘慈群等[7-8]认为：低渗透油藏渗流具有启动压力梯度，压力扰动的传播并非瞬时到达无穷远，其渗流规律就是一个动边界问题；页岩储集层具有与低渗油藏类似的动边界压力传播特性；压力方程的建立和求解较难，考虑动边界问题则更难，至今尚未见这方面的研究报道。虽然目前国内外对页岩气开发的研究已进入快速发展阶段，然而由于页岩气储集层条件复杂，对页岩气渗流和产能递减规律的研究，大多数研究成果给出的仅是页岩气渗流规律及其影响因素，尽管也有部分文献提出了具体的渗流方程，但考虑的因素尚少[9-12]，方程过于简化，难于更好地反映页岩气的低速强非线性渗流规律。为此，有必要揭示其页岩气的非线性流动规律，研究页岩气不稳定渗流特征，以便选取合理有效的开发方式和增产手段，为页岩气开发提供理论依据。

本文基于压力传播的稳定状态依次替换法对页岩气渗流动边界问题进行研究，进而建立考虑解吸、扩散、滑移及动边界影响的页岩气不稳定渗流数学模型，并推导和求解。结合中国南方某海相页岩储集层参数，分析页岩气不稳定渗流压力特征及其影响因素。

1 页岩气压力传播规律

在解决不稳定渗流压力动态的问题时，可以把不稳定渗流过程的每一瞬间状态看作是稳定的，这种方法称为稳定状态依次替换法[13]。

当页岩气投入开发、页岩储集层被打开后，形成的压力降将逐渐向外传播，设某时刻t，压力降传到R(t)处，在R(t)范围内形成压降漏斗，R(t)即为渗流过程中的扰动边界。R(t)随时间逐渐增大，压力降波及到的边缘为条件影响边缘，在该边缘上压力等于原始地层压力。

1.1 非线性渗流对动边界的影响

页岩储集层非常致密，主要为纳米—微米级孔隙，管壁和流体之间的微观作用力使气体在纳米—微米孔隙中流动时出现类似油藏启动压力梯度现象。页岩气的流动不仅有渗流过程，还存在扩散、滑移、解吸流动，气体流动总体表现为非线性流动，流动阻力比常规天然气大。滑脱效应附加了一种滑脱动力，但在驱动力小于气固间吸附作用所产生的阻力后，气体同样不能流动，即压力传播具有一定的动用范围。因此在不稳定渗流过程中压力扰动随时间延续逐渐向外传播，其边界条件也是一个动边界问题。

对于纳米—微米孔隙页岩储集层，气体在其中流动时，由于储集层渗透率极低，流动已偏离达西定律，扩散、滑移作用对储集层内气体流动影响增加。朱维耀等[14]建立了考虑扩散、滑移的纳米—微米孔隙气体流动方程：

则任一瞬间，地层压力为：

如图1所示，在地层中半径为r处取出厚度为h，宽度为dr的微小圆环体，其体积为2πrhdr，此单元体中游离态气体的原始质量为2πrhφiρidr。在给定时刻t，该单元体中残留气体质量为2πrhφρdr，因此，从单元体孔隙中采出的游离态气体质量为2πrh(φiρi-φρ)dr。

图1 页岩气平面径向流示意图

考虑吸附态气体的解吸，采出的总气体量为：

代入气体状态方程，联立（2）式、（3）式进一步求得：

气井流量按稳定渗流公式可写为：

假设页岩气井采用定产量生产，即qsc1为常数，则，联立（4）式、（5）式可得：

由此得出页岩储集层压力扰动传播动边界随时间变化的关系为：

中国南方下志留统龙马溪组海相页岩C气藏孔隙度为0.07，绝对渗透率0.000 5×10-3μm2；地层温度396.15 K，标准状态下气体压缩因子为1，真实气体压缩因子0.89；气体黏度0.027 mPa·s；外边界压力24 MPa；井筒半径0.1 m；气藏厚度30.5 m；平均解吸量3.370 7×10-10kg/(m3·s)；扩散系数8.406 7×10-7cm2/s，定产量200 m3/d。

利用C气藏基本参数，采用（7）式绘制不同渗透率条件下页岩气储集层动边界随时间变化关系图。由图2可见，同一时刻，渗透率越大，压力扰动传播动边界越远。在生产初期，压力扰动传播边界扩展较快，随着生产时间的延续，传播速度逐渐减慢。当渗透率大于0.5×10-3μm2时，压力快速传播到动边界；当渗透率小于等于0.5×10-3μm2时，压力随着时间的推进逐渐向外传播，因此页岩储集层中压力传播需考虑动边界的影响。

图2 不同渗透率条件下动边界随时间变化关系

1.2 页岩储集层压力传播动边界问题

部分学者认为在低渗透气藏渗流问题中，由于启动压力梯度的影响，压力并不是瞬间传播到无穷远，而是随着时间的推进逐渐向外传播，压力的传播边界称之为动边界。分析启动压力梯度和非线性渗流对压力传播边界的影响，认为只有在外加压力梯度大于启动压力梯度时，气体才发生流动。在页岩气开采过程中，由于启动压力梯度的存在，压力传播存在一定的范围。建立考虑启动压力梯度的稳态径向流常微分方程及其内边界定产、外边界定压的边界条件：

求解得到启动压力梯度下页岩储集层中直井稳态渗流表达式：

地层中采出的总气体量：

将（10）式代入考虑启动压力梯度的页岩储集层直井稳态渗流表达式（（9）式），得页岩气储集层压力扰动传播动边界随时间变化的关系为：

当Gt→0时，即得到不考虑启动压力梯度的页岩储集层压力扰动传播动边界与时间的变化关系：

图3为不同启动压力梯度下动边界随时间变化关系。可见，在同一时刻，启动压力梯度越大，压力扰动传播动边界越小。

图3 不同启动压力梯度下动边界随时间变化关系

2 页岩气不稳定渗流模型及求解

页岩气在储集层中的流动包括3个过程：甲烷从页岩基质表面解吸；甲烷通过页岩基质和微孔隙流动；甲烷在裂缝系统流动。本文考虑页岩储集层气体解吸、扩散、滑移等非线性渗流特征在不稳定渗流过程中对压力扰动传播动边界的影响，利用压力传播的稳定状态依次替换法，建立页岩气不稳定渗流模型，并进行推导和求解。

2.1 页岩气吸附-解吸模型

通过吸附-解吸实验得到不同平衡压力下页岩吸附-解吸过程的含气量变化。将得到的实验结果与吸附-解吸模型对比，结果表明，Langmuir模型拟合程度很高且形式简单，适于描述页岩气的吸附过程[15]。

Langmuir吸附模型[15]的数学表达式如下：

考虑瞬时平衡条件，解吸量可以表示为：

2.2 页岩气不稳定渗流模型

在考虑解吸、扩散、滑移作用的页岩气非线性渗流方程基础上[16-17]，基于天然气渗流的连续性方程、运动方程和状态方程，考虑不稳定渗流过程中压力扰动传播动边界的影响，引入动边界的模型，建立页岩气储集层不稳定渗流控制方程：

引入拟压力函数，并定义如下：

由（15）式、（17）式得到用拟压力表示的页岩气不稳定渗流的基本微分方程：

则（18）式可转换为：

2.3 页岩气不稳定渗流模型求解

当页岩气直井以某一恒定产量生产时，内边界定产，外边界定压，受动边界影响，泄压半径为当前时刻压力扰动传播到的距离，则其定解条件如下：

为方便求解，将（23）式用无因次量表示，无因次距离：

无因次拟时间：

无因次拟压力：

无因次流量：

则（23）式可化为：

引入拉普拉斯变换函数：

则（28）式变换为：

代入对应的定解条件：

经过逆变换，求得地层任意一点压力变化规律为：

其中βn为下式第n个正根：

当r=rw时，得井底压力变化规律：

3 页岩气不稳定渗流储集层压力分布特征

根据前面推导出的考虑解吸、扩散的页岩气不稳定渗流压力分布规律，结合C气藏参数，应用MATLAB编程计算，对页岩气不稳定渗流压力分布及其影响因素进行分析。

图4为不同时间地层压力分布曲线。由图可见，随着生产时间的延长，地层压力逐渐向外传播。在动边界影响范围内，页岩气储集层及气体释放弹性能，形成一个压降漏斗。动边界影响范围以外的地区，由于没有压力扰动，气体并不流动，且动边界的传播速度逐渐减慢。

图5为不同产量条件下井底压力随时间变化曲线。由图可见，井底压力随时间增加而降低，且减小趋势逐渐减缓，在前50 d，井底流压下降较快。同一时刻，产量越大，地层压力下降越多。

图4 不同生产时间地层压力分布曲线

图5 不同产气量条件下井底压力随生产时间变化曲线

图6为动边界对地层压力分布的影响。由图可见，考虑动边界影响时，压力扰动范围减小，地层压力分布下降减缓。对于超致密的纳微米孔隙页岩储集层，压力扰动随时间逐渐向外传播，且速度较慢，因此，考虑动边界影响的压力分布更贴近实际，更能准确地指导页岩气的生产。

图6 动边界对地层压力分布的影响

图7为解吸量对地层压力分布的影响。由图可见，考虑解吸量时地层压力较高，吸附气体的解吸使地层压力传播速度减慢，地层压力下降减缓。

图7 解吸气体对地层压力分布的影响

图8为扩散系数对地层压力分布的影响。由图可见，定产条件下，在动边界影响范围内，扩散系数越大，地层压力下降越慢，且扩散系数的影响逐渐减小。

图9为解吸、渗流、扩散、滑移及对产气量贡献率。由图可见，扩散、滑移在页岩气藏开采过程中占主要地位。在气藏开采过程中，随着开采的进行，渗流及解吸对产气量的贡献逐渐减小后趋于平稳，而扩散、滑移的贡献逐渐增加。生产初期，页岩气藏压降较小，基质中的游离气滑移、扩散，对产气量的贡献逐渐增加；随着生产的进行，解吸气不断释放，贡献逐渐趋于平稳。

图8 扩散系数对地层压力分布的影响

图9 渗流、解吸、扩散、滑移对产气量的贡献率

4 结论

基于页岩气吸附-解吸模型及纳米—微米孔隙气体流动模型，建立了考虑解吸、扩散、滑移的页岩气不稳定渗流数学模型。引入拟压力、拟时间函数及天然气扩散、解吸压缩系数，得到了页岩气不稳定渗流控制方程。

利用稳定状态依次替换法推导了考虑解吸、扩散、滑移的页岩气储集层压力扰动传播动边界随时间变化的关系式。研究表明，动边界随时间增加逐渐增大，渗透率越大，动边界传播越远。

通过拉普拉斯变换，推导求解了内边界定产、外边界为动边界的页岩气储集层不稳定渗流压力特征方程，得出了井底压力变化规律。结合中国南方某海相页岩气藏参数，计算分析了页岩气不稳定渗流储集层压力分布及其影响因素。研究结果表明，地层压力分布在动边界影响范围内形成压降漏斗，在动边界影响范围以外，没有压力扰动，气体并不流动；动边界随时间向外传播，且传播速度逐渐减慢。井底压力随时间增加而降低，且降低趋势逐渐减缓。与不考虑动边界的压降曲线对比，考虑动边界影响的地层压力传播慢，压力下降减缓。气体解吸扩散影响压力的传播速度，解吸使压力传播速度减慢，地层压力下降减缓；扩散系数越大，地层压力下降越慢，且扩散系数的影响逐渐减小。在气藏开采过程中，扩散、滑移对产气量的贡献逐渐增加，且占主要地位；渗流及解吸对产气量的贡献逐渐减小后趋于平稳。用本模型方法描述页岩气储集层的压力分布特征更符合实际。

符号注释：

a——与努森数Kn有关的修正系数（当0≤Kn＜0.001，a=0；0.001≤Kn＜0.1，a=1.2；0.1≤Kn＜10，a=1.34）；C——气体的等温压缩系数，Pa-1；Cd——解吸压缩系数，Pa-1；Ck——扩散压缩系数，Pa-1；Ct*——综合压缩系数，Pa-1；DK——扩散系数，m2/s；Ei——幂积分函数；G——启动压力梯度，Pa/m；h——气层厚度，m；I0——零阶第一类虚变量贝塞尔函数；I1——一阶第一类虚变量贝塞尔函数；J0——零阶第一类贝塞尔函数；J1——一阶第一类贝塞尔函数；K0——储集层绝对渗透率，m2；K0——零阶第二类虚变量贝塞尔函数；K1——一阶第二类虚变量贝塞尔函数；m——拟压力函数，Pa/s；p——储集层压力，Pa；pa——某一已知压力，Pa；pe——外边界压力，Pa；pL——Langmuir压力，表示吸附量为最大吸附量一半时的压力，Pa；psc——标准压力，Pa；pw——内边界压力，Pa；qd——单位体积页岩单位时间的解吸量，kg/（m3·s）；Qsc——地层中采出的总气体量，m3；qsc1——标准条件下气井流量，m3/s；qsc2——启动压力梯度影响标准条件下气井流量，m3/s；r——距井筒距离，m；re——气井供给半径，m；rw——井筒半径，m；R（t）——动边界，m；R1（t）——微尺度效应影响动边界，m；R2（t）——启动压力梯度影响动边界，m；s——拉普拉斯算子；t——生产时间，s；ta*——拟时间，s；T——地层温度，K；Tsc——标准状态下温度，K；v——气体渗流速度，m/s；Vd——单位体积基质累计解吸量，m3/m3；VE——气体的吸附量，m3/m3；VL——Langmuir体积，表示最大吸附量，m3/m3；x——渗流距离，m；Y0——零阶第二类贝塞尔函数；Y1——一阶第二类贝塞尔函数；Z——气体压缩因子，无因次；Zsc——标准状态下气体压缩因子，无因次；μ——气体黏度，Pa·s；ρ——气体密度，kg/m3；ρgsc——标准状态下气体密度，kg/m3；φ——储集层孔隙度，f；Ψ——非线性拟压力函数，Pa/s；——拉普拉斯变换非线性拟压力函数，Pa/s；Ψw——井底非线性拟压力，Pa/s。下标：D——无因次；i——原始地层条件；w——井底。

[1] 宁正福,王波,杨峰,等.页岩储集层微观渗流的微尺度效应[J].石油勘探与开发,2014,41(4):445-452.NING Zhengfu,WANG Bo,YANG Feng,et al.Microscale effect of microvadose in shale reservoirs[J].Petroleum Exploration and Development,2014,41(4):445-452.

[2] 邹才能,董大忠,王玉满,等.中国页岩气特征、挑战及前景（一）[J].石油勘探与开发,2015,42(6):689-701.ZOU Caineng,DONG Dazhong,WANG Yuman,et al.Shale gas in China:Characteristics,challenges and prospects (Ⅰ)[J].Petroleum Exploration and Development,2015,42(6):689-701.

[3] 杨峰,宁正福,胡昌蓬,等.页岩储层微观孔隙结构特征[J].石油学报,2013,34(2):301-311.YANG Feng,NING Zhengfu,HU Changpeng,et al.Characterization of microscopic pore structures in shale reservoirs[J].Acta Petrolei Sinica,2013,34(2):301-311.

[4] 钟太贤.中国南方海相页岩孔隙结构特征[J].天然气工业,2012,32(9):1-4.ZHONG Taixian.Characteristics of pore structure of marine shales in South China[J].Natural Gas Industry,2012,32(9):1-4.

[5] 李治平,李智锋.页岩气纳米级孔隙渗流动态特征[J].天然气工业,2012,32(4):50-53.LI Zhiping,LI Zhifeng.Dynamic characteristics of shale gas flow in nanoscale pores[J].Natural Gas Industry,2012,32(4):50-53.

[6] PASCAL H.Non-steady flow through porous media in the presence of a threshold gradient[J].Acta Mechanica,1981,39(3/4):207-224.

[7] 刘慈群.有起始比降固结问题的近似解[J].岩土工程学报,1982,4(3):107-109.LIU Ciqun.Approximate solution of a starting gradient ratio consolidation problem[J].Chinese Journal of Geotechnical Engineering,1982,4(3):107-109.

[8] 李凡华,刘慈群.含启动压力梯度的不定常渗流的压力动态分析[J].油气井测试,1997(1):1-4.LI Fanhua,LIU Ciqun.Pressure transient analysis for unsteady porous flow with start-up pressure derivative[J].Well Testing,1997(1):1-4.

[9] 尹虎,王新海,姜永,等.页岩气藏渗流数值模拟及井底压力动态分析[J].长江大学学报(自然科学版):理工,2012,9(8):68-71.YIN Hu,WANG Xinhai,JIANG Yong,et al.Numerical simulation of shale gas seepage and transient analysis of bottom-hole pressure[J].Journal of Yangtze University(Natural Science Edition):Sci & Eng,2012,9(8):68-71.

[10] 于荣泽,张晓伟,卞亚南,等.页岩气藏流动机理与产能影响因素分析[J].天然气工业,2012,32(9):10-15.YU Rongze,ZHANG Xiaowei,BIAN Ya’nan,et al.Flow mechanism of shale gas reservoirs and influential factors of their productivity[J].Natural Gas Industry,2012,32(9):10-15.

[11] 程远方,董丙响,时贤,等.页岩气藏三孔双渗模型的渗流机理[J].天然气工业,2012,32(9):44-47.CHENG Yuanfang,DONG Bingxiang,SHI Xian,et al.Seepage mechanism of a triple porosity/dual permeability model for shale gas reservoirs[J].Natural Gas Industry,2012,32(9):44-47.

[12] 李亚洲,李勇明,罗攀,等.页岩气渗流机理与产能研究[J].断块油气田,2013,20(2):186-190.LI Yazhou,LI Yongming,LUO Pan,et al.Study on seepage mechanism and productivity of shale gas[J].Fault-Block Oil & Gas Field,2013,20(2):186-190.

[13] 葛家理.油气层渗流力学[M].北京:石油工业出版社,1982.GE Jiali.The oil and gas flow through porous media[M].Beijing:Petroleum Industry Press,1982.

[14] 朱维耀,马千,邓佳,等.纳微米级孔隙气体流动数学模型及应用[J].北京科技大学学报,2014,36(6):709-715.ZHU Weiyao,MA Qian,DENG Jia,et al.Mathematical model and application of gas flow in nano-micron pores[J].Journal of University of Science and Technology Beijing,2014,36(6):709-715.

[15] 郭为,熊伟,高树生,等.页岩气等温吸附/解吸特征[J].中南大学学报(自然科学版),2013,44(7):2836-2840.GUO Wei,XIONG Wei,GAO Shusheng,et al.Isothermal adsorption/desorption characteristics of shale gas[J].Journal of Central South University(Science and Technology),2013,44(7):2836-2840.

[16] GERAMI S,POOLADI-DARVISH M,MORAD K,et al.Type curves for dry CBM reservoirs with equilibrium desorption[J].Journal of Canadian Petroleum Technology,2008,47(7):48-56.

[17] DENG Jia,ZHU Weiyao,MA Qian.A new seepage model for shale gas reservoir and productivity analysis of fractured well[J].Fuel,2014,124(15):232-240.

[18] AGARWAL R G,GARDNER D C,FUSSELL D D.Analyzing well production data using combined-type-curve and decline-curve analysis concepts[J].SPE Reservoir Evaluation & Engineering,1999,2(5):478-486.

（编辑 郭海莉）

Unstable seepage modeling and pressure propagation of shale gas reservoirs

ZHU Weiyao1,QI Qian1,MA Qian1,DENG Jia1,YUE Ming1,LIU Yuzhang2
(1.Civil and Environmental Engineering School,University of Science and Technology Beijing,Beijing 100083,China; 2.PetroChina Research Institute of Petroleum Exploration & Development,Beijing 100083,China)

Pressure disturbance propagation was investigated using the steady state replacement method,the relationship between moving boundary and time was obtained.An unstable seepage model in shale gas reservoirs was established considering the effects of desorption,diffusion,slip and moving boundary.Using Laplace transform,the pressure characteristics equation was solved for the condition of internal boundary being constant production and outer boundary being the moving boundary.Subsequently,combining the parameters of shale gas in southern China,unstable seepage pressure characteristics and its influence factors of shale gas reservoir were analyzed using MATLAB software.The results indicate that the pressure propagation is characterized by moving boundary effect during shale gas exploitation,which means that moving boundary is propagated outwards with the propagation velocity decreasing gradually.Under the effect of moving boundary or shale gas desorption,the pressure propagation velocity decreases and the reservoir pressure drop slows down.With the increasing of the diffusion coefficient,the reservoir pressure drop slows down and the effect of diffusion coefficient decreases gradually.In the process of gas reservoir exploitation,diffusion and slip contribute more and more to gas production,acting as the dominant factors,while the contribution of flow and desorption level off after decreasing.

shale gas; unstable seepage; pressure propagation law; mathematical model; moving boundary

国家重点基础研究（973）发展计划（2013CB228002）

TE37

A

1000-0747(2016)02-0261-07

10.11698/PED.2016.02.12

朱维耀（1960-），男，辽宁沈阳人，博士，北京科技大学教授，现从事渗流力学、非常规油气田开发方面的研究工作。地址：北京市海淀区学院路30号，北京科技大学土木与环境工程学院，邮政编码：100083。E-mail：weiyaook@sina.corn

2015-07-01

2016-02-20