考虑表面扩散的页岩气藏多级压裂水平井数值模拟

吴明录，丁明才，姚 军，徐思南，李 轩

(中国石油大学(华东) 石油工程学院,山东 青岛 266580)

引 言

随着世界常规油气资源开采趋于饱和，以页岩气为代表的非常规资源越来越越受到人们重视。页岩气储量丰富、分布广、稳产周期长，具有巨大的开采价值。

页岩气主要以吸附气和自由气的形式赋存于纳米级的页岩孔隙中[1]。页岩基质可分为有机质和无机质，其中有机质比表面积巨大[2]，吸附着大量的天然气。国内外学者对自由气在孔隙中渗流机理做了大量研究[2-6]，考虑页岩气在孔隙中的达西流动、滑脱效应、Knudsen扩散等运移机理，推导出与页岩基质固有渗透率有关的视渗透率模型。然而这些模型忽略了吸附气在孔隙中的运移，吸附气表面扩散也是页岩气的重要运移方式[7]。一些学者从微纳尺度出发，提出了多种考虑吸附层和表面扩散的页岩气视渗透率模型[8-15]，研究表明吸附气表面扩散会显著增大页岩气在纳米级多孔介质中的视渗透率。然而目前的页岩气藏压裂水平井数值模拟研究中都没有考虑吸附气表面扩散，这使得模拟结果与实际生产数据不符。

本文考虑自由气黏性流动、Knudsen扩散、吸附气解吸和表面扩散等渗流机理，建立基质系统渗流数学模型，在此基础上建立由基质系统、水力裂缝和诱导裂缝组成的完全耦合的页岩气藏多级压裂水平井渗流数学模型。采用嵌入式离散裂缝模型[16-21]和拟牛顿迭代法[22]求解，得到页岩气藏降压开采过程中压力分布和累积产气量。根据模拟结果分析吸附气表面扩散对压力分布和累积产气量的影响，为页岩气藏多级压裂水平井降压开采的建模和求解提供理论基础。

1 页岩气藏压裂水平井渗流数学模型建立与求解

1.1 页岩基质视渗透率模型

在纳米级的基质孔隙中，页岩气的Knudsen数[2]的表达式为

(1)

式中：r为页岩基质孔隙半径，m；λ为气体平均分子自由程[23]，表达式为

(2)

式中：pm为基质系统压力，Pa；μg为气体动力黏度，Pa·s；Zg为气体压缩因子，无因次；Mg为气体摩尔质量，g/mol；R为热力学常数，8.314J·mol-1·K-1；T为温度，K。

吴克柳等[12]提出页岩孔隙中吸附气表面扩散视渗透率模型，表达式为

(3)

式中：ζms为页岩气藏纳米孔隙中吸附气修正系数，表达式为

其中，rad为孔隙中心到吸附层距离，m；rt为平均孔隙半径，m；τ为纳米毛细管的迂曲度，无因次；φ为孔隙度，无因次。

Ds为吸附气表面扩散系数[12-14]，表达式为

其中：κ为吸附气拥塞率与迁移率的比值，无因次；H(1-κ)为分段函数，当κ>1时，值为0，0<κ<1时，值为1；Ds0为吸附气覆盖率为0时的表面扩散系数，m2/s。

Csc为吸附气浓度，kg/m3，表达式为

其中：dm为孔隙直径，m；NA为阿伏伽德罗常数，6.022×1023mol-1；θ为吸附气覆盖率，无因次，可以由储层压力pm和Langmuir压力pL表示，即

(4)

页岩孔隙中的自由气黏性流动和Knudsen扩散引起的视渗透率表达式[12]分别为

(5)

(6)

式中：kvs、kk分别为自由气黏性流动、Knudsen扩散引起的视渗透率，m2；Vstd为标准状况下天然气摩尔体积，m3(标)/mol；b为滑移系数，无因次；δ孔隙表面平均颗粒直径与平均孔隙直径的比值，无因次；Df为孔隙表面的分形维数，表征孔隙的粗糙程度，无因次；α为稀薄系数，与Knudsen数有关，表达式为

ζmb为页岩气藏纳米孔隙中自由气修正系数，表达式为

其中，α0，α1，β都为拟合参数值。

考虑黏性流动、Knudsen扩散和吸附气表面扩散时的页岩基质视渗透率[12]为

kpm=kv+kk+ks。

(7)

由达西定律可得页岩气藏基质单元渗流速度为

(8)

1.2 页岩基质系统渗流数学模型

页岩气主要以吸附态和游离态的形式赋存于页岩基质系统中，考虑自由气黏性流动、Knudsen扩散、吸附气表面扩散和解吸，建立基质系统连续性方程

(9)

式中：q为源汇项，kg/(m3·s)；qad为单位体积储层中页岩气吸附量，kg/m3，表达式为

(10)

其中：VL为单位质量页岩的Langmuir体积，m3/kg；pL为Langmuir压力，Pa；ρs为岩石密度，kg/m3；Vstd为标准状况下页岩气摩尔体积，m3/mol。

式(9)中，通过脉冲函数δh、δf判断水力裂缝、诱导裂缝与基质网格的相交情况，若相交δh、δf值为1，否则为0；qm-h，qm-f分别表示基质向水力裂缝、诱导裂缝的窜流量[20,21]，其表达式为

qm-h=ρgTm-h(pm-ph)，qm-f=ρgTm-f(pm-pf)。

(11)

其中，Tm-h，Tm-f分别为基质向水力裂缝、诱导裂缝的流动系数[16-20,24]。

页岩气密度ρg表达式为

(12)

其中，Zg为气体压缩因子，由Mahmoud[25]提出的关系式计算得到，即

式中，Tpr为天然气视临界温度，ppr为天然气视临界压力。

天然气压缩系数Cg也可以由Mahmoud[25]提出的关系式求得，即

(13)

将式(8)、(10)、(12)、(13)代入式(9)，即可得到页岩气在基质系统的渗流控制方程，即

(14)

1.3 裂缝系统渗流数学模型

诱导裂缝系统的连续性方程[24]为

(15)

由达西定律得到页岩气运动方程为

(16)

式中，pf为诱导裂缝系统压力，Pa。

由诱导裂缝系统的连续性方程(15)、运动方程(16)和状态方程(12)得到裂缝系统渗流控制方程为

(17)

同理可得到水力裂缝系统中页岩气渗流控制方程为

(18)

其中，ph为水力裂缝系统压力，Pa。

1.4 页岩气藏多级压裂水平井模型

如图1所示，针对压裂改造区(SRV)和未改造区(USRV)分别建立页岩气渗流数学模型。SRV区域由基质、诱导裂缝和水力裂缝组成，其渗流控制方程由式(14)、(17)和(18)组成；未改造区域中只含基质，渗流控制方程为式(14)。

图1 页岩气藏压裂水平井示意图Fig.1 Fracturing horizontal well in shale gas reservoir

考虑初始条件：

pm|t=0=pf|t=0=ph|t=0=pi；

(19)

内边界条件(定井底流压)：

pm|Γ1=pf|Γ1=ph|Γ1=pw；

(20)

外边界条件(封闭情形)：

(21)

其中，Г1，Г2为内、外边界。

1.5 模型求解

首先，采用结构化网格将储层剖分，裂缝系统采用二维的长方形网格，基质系统采用三维长方体网格。其次，将基质系统、诱导裂缝和水力裂缝渗流控制方程差分离散得到非线性方程组[22]，采用嵌入式离散裂缝模型[20-21]将诱导裂缝和水力裂缝嵌入基岩系统，并考虑诱导裂缝与水力裂缝间的窜流，通过文献[24，26]的井处理方法处理基质系统和裂缝系统的内边界条件。最后，由拟牛顿迭代方法[26-27]求解得水平井累积产量和储层压力分布。

模拟封闭页岩气藏中心一口水平井，垂直于水平井方向压裂出12条对称的水力裂缝，其他参数借鉴文献[12]、[17]，具体见表1。

表1 页岩气藏压裂水平井的参数值Tab.1 Parameters of fracturing horizontal well in shale gas reservoir

2 模拟结果

2.1 正确性验证

为验证本文方法的正确性，将水力裂缝简化至3条并忽略诱导裂缝，分别由本文使用的嵌入式离散裂缝模型(EDFM)与商业软件Comsol求解水平井井底压力，并将2种方法的计算结果进行对比。如图2所示，EDFM与Comsol的求解结果基本一致，证明了本文方法的正确性。

图2 EDFM与Comsol求解结果对比Fig.2 Comparison between the solving results of EDFM and Comsol

2.2 结果分析

图3为考虑和不考虑吸附气表面扩散时，页岩气藏在不同时间的压力分布。由图3可以看出：考虑表面扩散时近井压力下降更加明显，这是由于在考虑吸附气表面扩散时气藏的视渗透率增大，页岩气藏流动能力增强，在井底流压恒定的情况下产量增加的缘故。

图4为考虑和不考虑表面扩散时的页岩气藏累积产量曲线。由图可以发现，考虑表面扩散时，页岩气藏累积产量与不考虑时具有较大差距，说明吸附气的表面扩散对页岩气藏产量具有重要影响，在纳米孔页岩气藏数值模拟中不可忽略。

3 参数敏感性分析

为研究Langmuir体积VL、Langmuir压力pL、表面扩散系数Ds0和各种渗流方式对页岩气累积产量的影响，在以上参数基础上，取不同的VL、pL、r0和Ds0值进行了参数敏感性分析。

3.1 各种渗流方式的影响

图5—图7为不同孔隙半径下，各种渗流方式对累积产量的影响。图5为黏性流动对页岩气累积产量的影响规律。从图5中可以看出，当只考虑表面扩散和Knudsen扩散、不考虑黏性流动时，不同孔隙半径下的计算结果基本一致，说明孔隙半径对表面扩散和Knudsen扩散视渗透率影响较小。而且孔隙半径越大，黏性流动对结果的影响越大，这主要因为孔隙半径越大，黏性流动引起的视渗透率越大，而表面扩散和Knudsen扩散引起的视渗透率基本不变的缘故。

图3 降压开采时页岩气藏在不同时间的压力分布Fig.3 Pressure distributions at different time in depressurization development of shale gas reservoir

图4 页岩气累积产量曲线Fig.4 Cumulative production curves of shale gas

图5 黏性流动对页岩气累积产量的影响Fig.5 Effect of viscous flow on cumulative production of shale gas

图6为吸附气表面扩散对页岩气累积产量的影响规律。从图中可以看出，在孔隙半径较小(5 nm)时，考虑和不考虑表面扩散的计算结果差别很大，说明孔隙半径较小时，吸附气表面扩散是页岩气的主要渗流方式。随着孔隙半径增大，考虑和不考虑吸附气表面扩散时的计算结果差异变小，在孔隙半径为100 nm时，由吸附气表面扩散引起的计算差异基本可以忽略。这是由于随着孔隙半径增大，黏性流动引起的视渗透率增大(图5)，而吸附气表面扩散引起的视渗透率基本不变的缘故。

图6 表面扩散对页岩气累积产量的影响Fig.6 Effect of surface diffusion on cumulative production of shale gas

图7为Knudsen扩散对页岩气累积产量的影响。由图中可以看出，在相同的孔隙半径条件下，Knudsen扩散对计算结果无显著影响。这主要因为：与黏性流动、表面扩散相比，Knusen扩散引起的视渗透率很小(差一数量级)；而且孔隙半径的变化对Knudsen扩散引起的视渗透率改变也很小，故Knudsen扩散对页岩气累积产量的计算结果基本没有影响，这与文献[11]中得出的结论一致。

综合以上分析，孔隙半径较小时(r0=5 nm)，页岩气的主要运移方式是表面扩散；当孔隙半径较大时(r0=100 nm)，页岩气的主要运移方式是黏性流动；两种情况下Knudsen扩散的影响都较小。而现有文献基本都忽略了表面扩散的影响，夸大了Knudsen扩散的影响[2-6]，这将使纳米孔页岩气储层的流动模拟产生很大偏差。

图7 Knudsen扩散对页岩气累积产量的影响Fig.7 Effect of Knudsen diffusion on cumulative production of shale gas

3.2 表面扩散系数的影响

在页岩基质孔隙半径分别为5 nm和50 nm时，分别取不同的表面扩散系数Ds0计算累计产量，结果如图8所示。从图中可以看出，在相同的孔隙半径条件下，表面扩散系数越大，累积产量越高；在不同的孔隙半径条件下，孔隙半径越小，表面扩散系数的影响越大。由3.1中的分析可知，这是由于随着孔隙半径变小，气体的主要运移方式由黏性流动逐渐转变为表面扩散的缘故。

图8 表面扩散系数对页岩气累积产量的影响Fig.8 Effect of surface diffusion coefficient on cumulative production of shale gas

3.3 Langmuir体积的影响

图9为Langmuir体积VL对累积产量的影响规律。由图可见，随着VL增大，页岩气的累积产量越高。这是因为VL越大，单位体积页岩基质孔隙表面吸附的气体越多，从而开采过程中由压降引起的吸附气解吸和表面扩散的气体越多，故页岩气累积产量越高。

图9 Langmuir体积VL对页岩气累积产量的影响Fig.9 Effect of Langmuir volume VL on cumulative production of shale gas

3.4 Langmuir压力的影响

图10为Langmuir压力pL对累积产量的影响。由图可见，pL越大，页岩气累积产量越低且上升越缓慢。这可从两个方面进行分析：一方面，pL越大，表面扩散引起的视渗透率越小，因此页岩气在基质中的运移速度越慢，累积产量上升也越慢；另一方面，pL越大，页岩孔隙表面气体吸附量qad越小(式(10))，在相同时间内由压降引起的吸附气解吸和表面扩散的气体量也越少，故相同时间内页岩气累积产量越低。

图10 Langmuir压力pL对页岩气累积产量的影响Fig.10 Effect of Langmuir pressure pL on cumulative production of shale gas

4 结 论

(1)当页岩储层孔隙半径较小时，页岩气主要以表面扩散方式运移；当孔隙半径较大时，页岩气主要以黏性流动的方式运移；相对于表面扩散和黏性流动，Knudsen扩散的影响较小。因此，在对纳米级页岩气储层进行流动模拟时，应考虑表面扩散的影响，可忽略Knudsen扩散的影响。

(2)储层的孔隙半径很小时，气体的运移以表面扩散为主，表面扩散系数对计算结果的影响越大，且其值越大，累积产气量越高；随着孔隙半径逐渐增大，气体的运移方式逐渐由表面扩散转变为黏性流动，表面扩散系数的影响减弱。

(3)通过对Langmuir体积、Langmuir压力的参数敏感性分析发现，二者对页岩气累积产量都有较大影响，Langmuir体积越大或Langmuir压力越小，同条件下吸附气解吸和表面扩散的量越大，累积产量也越高。

(4)建立的页岩气藏多段压裂水平井渗流数学模型具有强非线性，模型中含有大量的水力裂缝和诱导裂缝，采用嵌入式离散裂缝模型和拟牛顿迭代法进行求解具有较高的求解速度和精度，因此具有较强的实用性。

参 考 文 献：

[1] 姚军,孙海,黄朝琴,等.页岩气藏开发中的关键力学问题[J].中国科学(物理学、力学、天文学),2013,43(12):1527-1547.

YAO Jun,SUN Hai,HUANG Zhaoqin,et al.Key mechanical problems in the development of shale gas reservoirs[J].Scientia Sinica (Physica,Mechanica & Astronomica),2013,43(12):1527-1547.

[2] JAVADPOUR F,FISHER D,UNSWORTH M.Nanoscale gas flow in shale gas sediments[J].Journal of Canadian Petroleum Technology,2007,46(10):55-61.

[3] KLINKENBERG L J.The permeability of porous media to liquids and gases[J].Socar Proceedings,1941,2(2):200-213.

[4] FLORENCE F A,RUSGING J,NEWSHAM K E,et al.Improved permeability prediction relations for low permeability sands[C].Rocky Mountain Oil & Gas Technology Symposium,2007.

[5] JAVADPOUR F.Nanopores and apparent permeability of gas flow in mudrocks (shales and siltstone)[J].Journal of Canadian Petroleum Technology,2009,48(8):16-21.

[6] WANG H Y,MARONGIU-PORCU M.A unified model of matrix permeability in shale gas formations[C].SPE Reservoir Simulation Symposium.Society of Petroleum Engineers,2015.

[7] FATHI E,AKKUTLU I Y.Nonlinear sorption kinetics and surface diffusion effects on gas transport in low-permeability formations[C].SPE Annual Technical Conference and Exhibition.Society of Petroleum Engineers,2009.

[8] XIONG X,DEVEGOWDA D,VILLAZON M,et al.A fully-coupled free and adsorptive phase transport model for shale gas reservoirs including non-Darcy flow effects[C].SPE annual technical conference and exhibition.Society of Petroleum Engineers,2012.

[9] Medve'I,ERN.Surface diffusion in porous media:a critical review[J].Microporous and Mesoporous Materials,2011,142(2):405-422.

[10] 盛茂,李根生,黄中伟,等.考虑表面扩散作用的页岩气瞬态流动模型[J].石油学报,2014,35(2):347-352.

SHENG Mao,LI Gensheng,HUANG Zhongwei,et al.Shale gas transient flow model with effects of diffusion[J].Acta Petrolei Sinica,2014,35(2):347-352.

[11] 吴克柳,李相方,陈掌星.页岩纳米孔吸附气表面扩散机理和数学模型[J].中国科学(技术科学),2015,45(5):525-540.

WU Keliu,LI Xiangfang,CHEN Zhangxing.The mechanism and mathematical model for the adsorbed gas surface diffusion in nanopores of shale gas reservoirs[J].Science China Press(Technology and Science),2015,45(5):525-540.

[12] WU K.Apparent permeability for gas flow in shale reservoirs coupling effects of gas diffusion and desorption[C].Unconventional Resources Technology Conference (URTEC),2014.

[13] WU K,LI X,WANG C,et al.Model for surface diffusion of adsorbed gas in nanopores of shale gas reservoirs[J].Industrial & Engineering Chemistry Research,2015,54(12):3225-3236.

[14] WU K,LI X,GUO C,et al.Adsorbed gas surface diffusion and bulk gas transport in nanopores of shale reservoirs with real gas effect-adsorption-mechanical coupling[C].SPE Reservoir Simulation Symposium.Society of Petroleum Engineers,2015.

[15] 岳陈军,张烈辉,赵玉龙,等.考虑表面扩散的页岩气渗透率两区复合解析模型[J].水动力学研究与进展(A辑),2016,31(3):362-371.

YUE Chenjun,ZHANG Liehui,ZHAO Yulong,et al.A dual-zone composite analytic model for shale gas permeability considering surface diffusion[J].Chinese Joural of Hydrodynamics(Series A),2016,31(3):362-371.

[16] LEE S H,JENSEN C L,LOUGH M F.An efficient finite difference model for flow in a reservoir with multiple length-scale fractures[C].SPE Annual Technical Conference and Exhibition.Society of Petroleum Engineers,1999.

[17] LEE S H,LOUGH M F,JENSEN C L.Hierarchical modeling of flow in naturally fractured formations with multiple length scales[J].Water Resources Research,2001,37(3):443-455.

[18] LI L,LEE S H.Efficient field-scale simulation of black oil in a naturally fractured reservoir through discrete fracture networks and homogenized media[J].SPE Reservoir Evaluation & Engineering,2008,11(4):750-758.

[19] MOINFAR A,VARAVEI A,SEPEHRNOORI K,et al.Development of an efficient embedded discrete fracture model for 3D compositional reservoir simulation in fractured reservoirs[J].SPE Journal,2014,19(2):289-303.

[20] MOINFAR A.Development of An Efficient Embedded Discrete Fracture Model for 3D Compositional Reservoir Simulation in Fractured Reservoirs[D].The University of Texas at Austin,2013.

[21] SHAKIBA M.Modeling and Simulation of Fluid Flow in Naturally and Hydraulically Fractured Reservoirs Using Embedded Discrete Fracture Model[D].The University of Texas at Austin,2014.

[22] 关治,陆金甫.数值分析基础[M].北京:高等教育出版社,1998.

GUAN Zhi,LU Jinfu.Basis of Numerical Analysis[M].Beijing:Higher Education Press,1998.

[23] CIVAN F.Effective correlation of apparent gas permeability in tight Porous Media[J].Transport in Porous Media,2010,82(2):375-384.

[24] ZHU G,YAO J,SUN H,et al.The numerical simulation of thermal recovery based on hydraulic fracture heating technology in shale gas reservoir[J].Journal of Natural Gas Science and Engineering,2016,28(1):305-316.

[25] MAHMOUD M.Development of a new correlation of gas compressibility factor (Z-factor) for high pressure gas reservoirs[J].Journal of Energy Resources Technology,2014,136(1):2280-2285.

[26] 张烈辉.油气藏数值模拟基本原理[M].北京:石油工业出版社,2005.

ZHANG Liehui.Basic Principles of Numerical Simulation of Oil And Gas Reservoirs[M].Beijing:Petroleum Industry Press,2005.

[27] MOINFAR A,VARAVEI A,SEPEHRNOORI K,et al.Development of a coupled dual continuum and discrete fracture model for the simulation of unconventional reservoirs[C].SPE Reservoir Simulation Symposium.Society of Petroleum Engineers,2013.