技术领域
[0001] 本
发明涉及一种计算页岩气藏复杂裂缝网络非稳态产量的方法，属于矿物资源技术领域。
背景技术
[0002] 页岩气资源潜
力巨大，近年来受到了人们的广泛关注。随着钻井以及压裂工艺的迅速提高，页岩气已经成为油气勘探开发的重要组成部分。
水平井压裂形成复杂缝网是开发页岩气的关键技术。页岩储层经过压裂施工后会在压裂改造区域会形成复杂裂缝网络，而在未改造的区域仍然为基质孔隙区域。页岩气压裂水平井在生产过程中，根据页岩气的流动特征，可以把压裂后的页岩气井渗流场分为内外两个流场，内部流场为体积压裂形成复杂缝网区域的达西流动，而在外部流场为页岩基质区域受到黏性流、扩散、
解吸和滑脱等多重运移机制组成；并且页岩气压裂水平井的生产是一个不稳定渗流过程，而且初期产量递减十分明显，只有到了生产后期压力波传到边界时，产量才趋于稳定。并且在生产过程中，储
层压力下降还会显著影响页岩气体在基质和缝网中的流动。
[0003] 目前常用的计算页岩气压裂水平井产量的源函数方法一开始就假设由于生产导致的压力传播一开始就传递到储层边界，这与实际情况不吻合；且计算的页岩基质渗透率建立的是单毛细管模型，与
地层实际情况不符；复杂裂缝网络渗透率通过实验测得，不能体现实际复杂裂缝网络实际情况，从而使得计算结果与实际生产误差较大。
发明内容
[0004] 针对
现有技术中的存在的不足，本发明提供了更加贴近地层实际情况、计算结果更加精确的一种计算页岩气藏复杂裂缝网络非稳态产量的方法。
[0005] 本发明的技术方案如下：
[0006] 一种计算页岩气藏复杂裂缝网络非稳态产量的方法，其特征在于，包括以下步骤：
[0007] S1、收集页岩气藏地层、页岩气的基本参数，建立页岩基质多孔介质毛细管渗透率模型；
[0008] S2、收集页岩气藏微
地震裂缝成像数据，根据微地震裂缝成像数据，识别出裂缝条数、
位置及形态特征；运用盒子法获取表征
岩石复杂裂缝网络的分形维数；根据求取的分形维数，利用蒙特卡洛随机建模理论生成二维离
散裂缝网络模型；通过立方定律计算复杂裂缝网络渗透率；
[0009] S3、假设储层等厚、圆形边界、半径为re，体积压裂区域改造充分、圆形边界、半径为rf、与储层厚度相同；考虑页岩气藏压裂水平井外区为基质区域渗流区域、内区为复杂缝网区域流动区域，建立页岩气藏压裂水平井稳态产能模型；建立页岩气压力传播半径方程、页岩气藏物质平衡方程，运用连续拟稳定法建立页岩气藏压裂水平井非稳态模型，求解页岩气藏压裂水平井非稳态产量。
[0010] 进一步的，所述页岩气藏地层、页岩气的基本参数包括储层边界压力、储层
温度、供给半径、基质孔隙度、裂缝孔隙度、体积压裂区半径、井眼半径、储层厚度、页岩气
粘度、气体类型、分子碰撞直径、气体分子直径、气体常数、气体摩尔
质量、切向动量调节系数、气体分子
密度、平均压力、表面最大浓度、朗格缪尔压力和表面扩散系数。
[0011] 进一步的，所述步骤S1中页岩基质多孔介质毛细管渗透率模型是根据毛细管大小的分布不同分别计算不同尺度下毛细管的渗透率，然后对不同尺度下的毛细管渗透率进行
叠加，从而得到整个
岩心的渗透率大小，其表达式如下：
[0012]
[0013] 式中：kapp—整个岩心的表观渗透率，m2； —基质孔隙度，无量纲；kapp，i—不同尺度毛细管的表观渗透率，m2；λi—不同尺度毛细管的分布
频率，无量纲；i，n—计数单位，无量纲；τ—迂曲度，无量纲。其中迂曲度计算公式如下：
[0014]
[0015] 式中：m为迂曲度拟合参数，取值0.77。
[0016] 进一步的，所述步骤S1中页岩基质多孔介质毛细管渗透率模型，根据达西定律求得页岩气藏整个岩心的表观渗透率的步骤是利用页岩气储层参数计算得到克努森系数，根据克努森系数判断毛细管中的流态，利用毛细管中的流态和运移机制建立气体质量运移方程，接下来建立统一的质量运移方程，利用统一的质量运移方程得到毛细管的表观渗透率，最后进行叠加求得整个岩心的表观渗透率；
[0017] 其中克努森系数的计算公式如下：
[0018]
[0019] 式中：Kn—克努森系数，无因次；kB—玻尔兹曼常数，1.3805×10-23J/K；T—页岩储层温度，K；π—常数，3.14；δ—气体分子碰撞直径，m；p—储层压力，Pa；d—毛细管直径，m。
[0020] 进一步的，所述步骤S2中复杂裂缝网络渗透率通过微地震数据地震成像获取分形维数，通过蒙特卡洛随机建模理论建立二维离散裂缝网络模型，并进一步通过立方定律来计算压裂施工后区域内复杂缝网的渗透率，其表达式如下：
[0021]
[0022]
[0023] 式中：kx为裂缝在x方向上的渗透率，m2；fx为在x方向上正
应力产生的裂缝的等效频率，1/m；br为裂缝的初始平均开度，m；bmax为受压下裂缝最大开度，m；αx为x方向上的压力系数，无量纲；σx为x方向上正应力大小，MPa；αy为y方向上的压力系数，无量纲；σy为y方向上正应力大小，MPa；fdx为在x方向上剪应力产生的裂缝的等效频率，1/m；dmax为裂缝发生剪切破坏后的最大开度，m；γx为x方向剪切压力系数，无量纲；γy为y方向剪切压力系数，无量纲；k为x方向于y方向上的应力比值，无量纲；kc为裂缝
变形的
临界压力比，无量纲；ky为裂缝在y方向上的渗透率，m2；fy为在y方向上正应力产生的裂缝的等效频率，1/m；βx为在xy下的函数斜率，无量纲；βy为在yx下的函数斜率，无量纲；fdy为在y方向上剪应力产生的裂缝的等效频率，1/m。
[0024] 进一步的，所述步骤S2中分形维数通过盒子法计算，盒子法是用正方形的格子去
覆盖复杂缝网区域，盒子的大小是变化的；给定盒子的尺码δ，可以数出覆盖住缝网区域所需的总盒子数目N；假设第i步覆盖使用δi×δi的格子，所需盒子数目为Ni(δi)，在第i+1步需要的格子(δi+1)×(δi+1)，则需要盒子数目Ni+1(δi+1)。得到任意两个尺度下所需盒子数之比与尺码之比的关系如下：
[0025]
[0026] 推广到：
[0027] N＝aδ-D
[0028] 得到分形维数表达式为：
[0029]
[0030] 式中：D为分形维数，无量纲；N为盒子数目，无量纲；δ为盒子尺码，无量纲；i为计数符号；
[0031] 所述步骤S2的二维离散裂缝网络模型的建模方法为随机建模方法，所述二维离散裂缝网络模型的生成由以下步骤组成：
[0032] A、模拟区域生成：模拟区域根据页岩储层压裂过程中裂缝系统分布范围来确定；裂缝网络模拟时，先生成裂缝的中心点；考虑到有裂缝只有部分在模拟区域内，中心点的位置可能在模拟区域外，为了减少边界效应的影响，可以适当扩大模拟区域；为了便于区别，扩大后的模拟区域称为生成域，而原来的模拟区域称为分析区；为保证模拟结果可靠，生成域的每一个侧面距分析域的距离大于裂缝出现的最大半径；
[0033] B、裂缝数目确定：模拟区域大小确定后，根据裂缝的面密度来确定该区域内每组裂缝的数目；
[0034] C、裂缝空间位置表征：对于二维离散裂缝网络模型，裂缝位置可以由中心点及迹长，或者裂缝的两个端点完全确定，即：
[0035]
[0036] D、裂缝几何参数：a按均匀分布生成裂缝中心点坐标随机变量x0、y0；b生成裂缝的倾向及倾
角；c生成裂缝的迹长；d生成裂缝宽度，因裂缝宽度较难确定，通常把每组裂缝的宽度设为一常量；
[0037] E、裂缝网络模型生成及简化：当各组裂缝的数目、中心点坐标、产状、半径及宽度确定后，裂生成域内也就完全确定了；生成域的裂缝生成确定后，根据分析的几何边界来生成裂缝网络模型。
[0038] 进一步的，所述步骤S2中复杂裂缝网络渗透率包括正常张开型裂缝渗透率、剪切型裂缝渗透率，计算渗透率过程中，正应力和剪应力分别进行计算，最后将正应力和剪应力诱导产生裂缝所具有渗透率进行叠加再复合得到，表达式如下：
[0039] kx＝knx+kdx
[0040] ky＝kny+kdy
[0041]
[0042] 式中：kx——在x方向上复杂缝网的渗透率，10-3μm2；knx——在x方向上由于裂缝正常闭合后的渗透率，10-3μm2；kdx——在x方向上由于裂缝发生剪切破坏产生的渗透率，10-3μm2；ky——在y方向上复杂缝网的渗透率，10-3μm2；kny——在y方向上由于裂缝正常闭合后的渗透率，10-3μm2；kdy——在y方向上由于裂缝发生剪切破坏产生的渗透率，10-3μm2；Kf—体积压裂区复杂裂缝网络渗透率，m2。
[0043] 进一步的，所述步骤S3中页岩气藏压裂水平井稳态产能模型将页岩气的渗流场分为内外两个流场，内部流场为体积压裂区的达西流动，外部流场为页岩基质中的多尺度毛细管流动；
[0044] 其中外部流场的气体流动方程表达式如下：
[0045]
[0046] 其中
[0047] km＝kapp
[0048]
[0049] 内部流场的流动方程如下：
[0050]
[0051] 式中：ν—
流体流速，m/s；Kf—体积压裂区复杂裂缝网络渗透率，m2；μ—气体黏度，mPa·s； —沿程压力梯度，MPa/m；
[0052] 页岩气藏压裂水平井的稳态产量方程表达式如下：
[0053]
[0054] 其中
[0055]
[0056]
[0057] 式中：qsc—页岩气标况下的产量，m3/d；pf—体积压裂区外边界压力，MPa；pwf—井底流压，MPa； —基质孔隙度，无量纲；pe—储层边界压力，MPa；kapp，i—不同尺度毛细管的表观渗透率，m2；μ—页岩气粘度，mPa·s；z—压缩因子，无量纲；psc—标准状况下的压力，MPa；rw—井眼半径，m；rf—体积压裂区的半径，m；re—供给半径， m；Kf—体积压裂区复杂裂缝网络渗透率，m2；h—储层厚度，m；Tsc—标准状况下的温度，K；kapp—整个岩心的表观渗透率，m2。
[0058] 进一步的，所述步骤S3中页岩气藏压裂水平井非稳态模型把页岩气的流动分为体积压裂区和压裂区到储层边界页岩基质内外两个渗流场，页岩气的流动为径向流，同时考虑页岩气的解吸，得到页岩气非稳态下的产量；
[0059] 单位质量的页岩基质在压力发生改变时解吸出的页岩气体积可由Langmuir等温
吸附解吸表达式得到，如下：
[0060]
[0061] 式中：Q—单位质量页岩基质解吸页岩气体积，m3/kg；pi—原始地层压力，MPa；
[0062] 在地面标准状况下，页岩气藏基质中的游离气量为：
[0063]
[0064] 式中：Gm—地下储集页岩气体积折算到地面标准状况下的体积，m3；Vm—页岩基质3
储集空间的体积，m；Bgi—页岩气在原始地层压力和温度下的体积系数，无量纲；
[0065] 页岩气藏基质的体积为：
[0066]
[0067] 式中：Vt—页岩储层视体积，m3；Sw—束缚水
饱和度，无量纲；
[0068] 根据朗格缪尔等温吸附理论，页岩气藏的吸附气体积为：
[0069]
[0070] Ga—页岩气基质中吸附气在地面标准状况下的体积，m3；ρs—页岩基质的密度，kg/m3；
[0071] 根据上式，进一步求得：
[0072]
[0073] 随着页岩气的不断被采出，地层压力下降，同时束缚水膨胀，基质孔隙收缩，裂缝中气体膨胀，令裂缝中的游离气为Gf。则储层中剩余游离气量为：
[0074] 基质中：
[0075] 改造区内：
[0076] 式中：G'm—基质剩余游离气在地面标准状况下的体积，m3；G'f—裂缝剩余游离气在地面标准状况下的体积，m3；Gf—裂缝中的游离气在标准状况下的体积，m3；Cm—页岩基质的压缩系数，MPa-1；Bg—当前压力下的页岩气体积系数，无量纲；
[0077] 根据朗格缪尔等温吸附理论，可得到页岩气藏剩余的吸附气量：
[0078]
[0079] 根据物质平衡原理，原始地层压力下的基质中的游离气量、原始地层压力下裂缝中的游离气量、原始地层压力下的吸附气量之和等于采出气量、当前地层压力下基质中的游离气量、当前地层压力下裂缝中的游离气量、当前地层压力下的吸附气量之和，即：
[0080] Gm+Gf+Ga＝Gp+G'm+G'f+G′a
[0081] 得到页岩气藏物质平衡方程表达式如下：
[0082]
[0083] 式中：Gm—地下储集页岩气体积折算到地面标准状况下的体积，m3；Gf—裂缝中的3 3
游离气在标准状况下的体积，m ；Gp—页岩气累计产出量，m ； —基质孔隙度，无量纲； Sw—束缚水饱和度，无量纲；Bgi—页岩气在原始地层压力和温度下的体积系数，无量纲；
ρs—页岩基质的密度，kg/m3；pi—原始地层压力，MPa；VL为朗格缪尔压力，Pa；pL为朗格缪尔体积，m3/kg；
[0084] 其中：
[0085]
[0086] 页岩气压裂水平井生产过程中的压力传播半径方程表达式如下：
[0087]
[0088] 式中：R—t时刻的压力传播半径，m；K—储层渗透率，m2；t—生产时间，d； —基质孔隙度，无量纲；μ—气体黏度，mPa·s；Ct—储层的综合压缩系数，MPa-1。
[0089] 所述连续拟稳定法步骤如下：求得初始产量，假定一个时间步长，认为在该时间步长内的渗流为稳定渗流，计算出累计产量，根据物质平衡方程，算出在该时间段内压力波所波及到的区域的平均地层压力，以该平均地层压力作为下一时间步长的边界压力，计算下一时间步长下的产量，以此类推，便可得到产量与时间的关系，即非稳态下的产量。
[0090] 本发明的有益效果为：
[0091] 本发明方法充分考虑地层岩芯由多尺度毛细管组成，首先求得计算不同尺度下毛细管的渗透率，然后对不同尺度下的毛细管渗透率进行叠加，得到整个岩心的渗透率更加贴近实际情况；通过页岩气藏压裂水平井微地震裂缝成像数据，运用盒子法获取分形维数，建立二维离散裂缝网络模型，利用改进的立方定律计算复杂裂缝网络渗透率；克服源函数等方法的
缺陷，利用连续拟稳定法，考虑压力波在生产过程中逐渐往外传播，压力波的传播速度与储层物性、流体性质和生产时间紧密相关，分析了在页岩气压裂水平井生产过程中在过渡流和边界流条件下页岩储层压力的变化规律；并且使用连续拟稳定法计算非稳态产量形式简单，容易实现程序化，同时避免了源函数等数学方法再求解产量时在生产初期压力波就传递到储层边界的不合理假设。
附图说明
[0092] 为了更清楚地说明本发明实施方式的技术方案，下面将对实施方式中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些
实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它相关的附图。
[0093] 图1本发明提供的页岩基质表观渗透率；
[0094] 图2本发明提供的页岩压裂水平井微地震反演裂缝分布图；
[0095] 图3本发明提供的裂缝分布图；
[0096] 图4本发明提供的分形维数计算图；
[0097] 图5本发明提供的裂缝网络渗透率与裂缝区域平均应力的关系图；
[0098] 图6本发明提供的裂缝网络渗透率与x与y方向应力的关系图；
[0099] 图7本发明提供的产量计算结果与Eclipse模拟结果对比图；
[0100] 图8本发明提供的不同分形维数下页岩气压裂水平井累产的曲线图；
[0101] 图9本发明提供的不同压裂半径下页岩气压裂水平井累产的曲线图；
[0102] 图10本发明提供的不同裂缝开度下，页岩气压裂水平井累产的曲线图。
具体实施方式
[0103] 下面结合附图和实施例对本发明进一步说明。
[0104] 一种计算页岩气藏复杂裂缝网络非稳态产量的方法，包括以下步骤：
[0105] S1、收集页岩气藏地层、页岩气的基本参数，建立页岩基质多孔介质毛细管渗透率模型；
[0106] S2、收集页岩气藏微地震裂缝成像数据，根据微地震裂缝成像数据，识别出裂缝条数、位置及形态特征；运用盒子法获取表征岩石复杂裂缝网络的分形维数；根据求取的分形维数，利用蒙特卡洛随机建模理论生成二维离散裂缝网络模型；通过立方定律计算复杂裂缝网络渗透率；
[0107] S3、假设储层等厚、圆形边界、半径为re，体积压裂区域改造充分、圆形边界、半径为rf、与储层厚度相同；考虑页岩气藏压裂水平井外区为基质区域渗流区域、内区为复杂缝网区域流动区域，建立页岩气藏压裂水平井稳态产能模型；根据页岩气压力传播半径方程、页岩气藏物质平衡方程，运用连续拟稳定法建立页岩气藏压裂水平井非稳态模型。
[0108] 在本发明提供的实施例中，所述页岩气藏地层、页岩气的基本参数包括储层边界压力、储层温度、供给半径、基质孔隙度、裂缝孔隙度、体积压裂区半径、井眼半径、储层厚度、页岩气粘度、气体类型、分子碰撞直径、气体分子直径、气体常数、气体摩尔质量、切向动量调节系数、气体分子密度、平均压力、表面最大浓度、朗格缪尔压力和表面扩散系数。
[0109] 在本发明提供的实施例中，所述步骤S1中页岩基质多孔介质毛细管渗透率模型是根据毛细管大小的分布不同分别计算不同尺度下毛细管的渗透率，然后对不同尺度下的毛细管渗透率进行叠加，从而得到整个岩心的渗透率大小，其表达式如下：
[0110]
[0111] 式中：kapp—整个岩心的表观渗透率，m2； —基质孔隙度，无量纲；kapp，i—不同尺度毛细管的表观渗透率，m2；λi—不同尺度毛细管的分布频率，无量纲；i，n—计数单位，无量纲；τ—迂曲度，无量纲。其中迂曲度计算公式如下：
[0112]
[0113] 式中：m为迂曲度拟合参数，取值0.77。
[0114] 在本发明提供的实施例中，所述步骤S1中页岩基质多孔介质毛细管渗透率模型，根据达西定律求得页岩气藏整个岩心的表观渗透率的步骤是利用页岩气储层参数计算得到克努森系数，根据克努森系数判断毛细管中的流态，利用毛细管中的流态和运移机制建立气体质量运移方程，接下来建立统一的质量运移方程，利用统一的质量运移方程得到毛细管的表观渗透率，最后进行叠加求得整个岩心的表观渗透率；
[0115] 其中克努森系数的计算公式如下：
[0116]
[0117] 式中：Kn—克努森系数，无因次；kB—玻尔兹曼常数，1.3805×10-23J/K；p—储层压力， Pa；T—页岩储层温度，K；π—常数，3.14；δ—气体分子碰撞直径，m；d—毛细管直径， m；
[0118] 其中运移机制包括游离气运移、吸附气运移，所述游离气运移包括黏性流动、滑脱流动、 Knudsen扩散，所述吸附气运移包括吸附气解吸、吸附气表面扩散；
[0119] (1)游离态页岩气质量运移方程
[0120] 游离态页岩气(游离气)赋存于基质孔隙以及裂缝中，主要包括黏性流动、滑脱流动、 Knudsen扩散作用。
[0121] ①粘性流质量运移方程
[0122] 页岩发育有大量的
纳米级孔隙，可将纳米孔视为毛管模型，页岩则可视为毛管和基质组成。当页岩气体Knudsen数远小于1时，气体分子的运动主要受分子间碰撞支配，此时分子与壁面的碰撞较少，气体分子间的相互作用要比气体分子与孔隙表面(孔隙壁)的碰撞频繁得多，气体以连续流动为主，可采用黏性流质量运移方程描述。当不考虑吸附气存在对毛管半径的影响时，对于吼道半径为r的单根毛管，其固有渗透率计算公式如下：
[0123]
[0124] 式中：kD—固有渗透率，m2；r—页岩毛细管半径，m，r＝d/2。
[0125] 在单组分气体之间存在压力梯度所引起的黏性流动，可以用达西定律来表示描述黏性流的质量运移方程，表达式如下：
[0126]
[0127] 式中：Jvicious—黏性流质量流量，kg/(m2·s)；ρ—气体密度，kg/m3；μ—气体黏度，Pa·s； p—储层压力，Pa；—压力梯度算子符号，无量纲。
[0128] 对于气体在
纳米管中的运移，当考虑吸附态气体存在对纳米孔半径的影响时，纳米孔吼有效半径减小，因此考虑吸附气影响时纳米孔吼的有效半径可表达为：
[0129]
[0130] 式中：re—考虑吸附气影响的孔隙有效半径，m；dm—气体分子直径，m；pL—Langmuir 压力，Pa。
[0131] 将式(6)代入式(5)，可以得到：
[0132]
[0133] ②滑脱效应质量运移方程
[0134] 当页岩孔隙尺度减小，或者气体压力降低、气体分子自由程增加，气体分子自由程与孔隙直径的尺度具有可比性，气体分子与孔隙壁面的碰撞不可忽略。在0.001＜Kn＜0.1时，由于壁面页岩气分子速度不再为零，此时存在滑脱现象，且考虑滑脱效应时的页岩储层渗透率计算公式如下：
[0135]
[0136]
[0137] 式中：kslip—考虑滑脱效应渗透率，m2；paver—通过实验测试岩心渗透率时的进出口平均压力，Pa；bk—滑脱因子，Pa。
[0138] 为了能将滑脱效应在渗流方程中体现，引入滑脱因子来修正纳米孔隙滑脱效应，将式(9) 代入式(8)：
[0139]
[0140] 式中：R—气体常数，J/(mol·K)；M—气体摩尔质量，kg/mol；pavg—平均压力(在圆形单管中为进口、出口平均压力)，Pa；α—切向动量调节系数，无因次，取值为0～1。
[0141] 因此考虑滑脱效应时的质量运移方程可以写为：
[0142]
[0143] 式中：Jslip—滑脱流质量流量，kg/(m2·s)。
[0144] ③克努森扩散质量运移方程
[0145] 当孔道直径减少或者分子
平均自由程增加(在低压下)，Kn＞10时，气体分子更容易与孔隙壁面发生碰撞而不是与其他气体分子发生碰撞，这意味着气体分子达到了几乎能独立于彼此的点，称为Knudsen扩散。
[0146] Jknudsen＝α0νρ (12)
[0147] 式中：α0—无因次概率系数，无因次；ν—平均分子速度，m/s；ρ—气体分子密度，3
kg/m。
[0148] 当圆管两端都有气体，圆管传输的净流量与圆管两端的气体密度成正比，式(12)可写为：
[0149] Jknudsen＝αν(ρin-ρout) (13)
[0150] 式中：ρin——圆管进口处气体密度，kg/m3；ρout——圆管出口处气体密度，kg/m3。
[0151] 根据
气体动力学理论，气体的平均分子运动速度为：
[0152]
[0153] 对于直径为d长度为L的圆形长直管(L>>d)，α＝d/3L，将式(14)代入到式(13)中，可以得到：
[0154]
[0155] 将式(15)写为偏微分形式：
[0156]
[0157] 式(16)也可写为气体浓度C的形式，即为：
[0158]
[0159] 而纳米孔隙中的Knudsen扩散系数Dk，表达式如下所示：
[0160]
[0161] 式中：Dk—克努森扩散系数，m2/s。
[0162] 而气体密度表达形式可写为：
[0163]
[0164] 将式(18)、(17)代入式(15)，可得：
[0165]
[0166] 因此，Knudsen扩散质量运移方程可表述为：
[0167]
[0168] 式中：Jknudsen—克努森扩散质量流量，kg/(m2·s)。
[0169] (2)吸附态页岩气运移质量方程
[0170] 吸附态页岩气(吸附气)赋存于孔隙壁面和页岩固体颗粒表面，主要发生解吸附作用和表面扩散作用。假设在一定温度和压力条件下，壁面吸附气和自由气处于瞬间动态平衡，采用 Langmuir等温吸附模型，吸附质量的表达形式为：
[0171]
[0172] 式中：qads—页岩单位体积的吸附量，kg/m3；Vstd—页岩气标况下摩尔体积，m3/mol；
[0173] 页岩气井在生产过程中，地层压力逐渐下降，考虑t1时刻地层压力为p1，t2时刻地层压力为p2，则可计算出地层压力由p1下降为p2时吸附态页岩气的解吸量：
[0174]
[0175] 式中：Δqads—由于压力下降产生的吸附态页岩气解析量，kg/m3；VL—朗格缪尔体积， m3/kg。
[0176] 页岩气在微纳米孔隙表面不仅存在解吸附效应，还存在沿吸附壁面的运移，即表面扩散作用。不同于压力梯度或浓度梯度作用的其他运移方式，页岩气表面扩散在吸附势场的作用下发生运移，影响页岩气表面扩散的因素很多，包括压力、温度、纳米孔壁面属性、页岩气体分子属性、页岩气体分子与纳米孔壁面相互作用等。
[0177] 当表面扩散气体运移方程表达为浓度梯度的形式时，等于表面扩散系数与浓度梯度的乘积形式，其表面扩散质量流量计算公式如下所示：
[0178]
[0179] 式中：Jsurface—表面扩散质量流量，kg/(m2·s)；Cs—孔壁面吸附气浓度，mol/m3；Ds—表面扩散系数，m2/s；l—孔壁长度，m。
[0180] 吸附气体覆盖率θ可表示为：
[0181]
[0182] 式中：θ—吸附气体覆盖率，无因次；Csmax—吸附气最大吸附浓度，mol/m3；V—单位质量页岩实际吸附气体积，m3/kg；Csmax—吸附气最大吸附浓度，mol/m。
[0183] 式(25)可进一步改写为：
[0184]
[0185] 将式(26)代入式(24)，可得到满足Langmuir等温吸附方程的页岩气表面扩散质量运移方程：
[0186]
[0187] 在本发明提供的实施例中，考虑游离态页岩气粘性流、滑脱流、Knudsen扩散和吸附态页岩气的解吸、表面扩散作用，其总的传输质量为这几种运移模式引起的传输质量的叠加之和，并引入贡献系数ε，建立统一的的质量运移方程，式(28)：
[0188]
[0189] 其中，F为滑脱系数，其计算公式如下：
[0190]
[0191] 式中：Jtol—总的质量流量，kg/(m2·s)；Jvicious—黏性流质量流量，kg/(m2·s)；Jslip—滑脱效应质量流量，kg/(m2·s)；Jknudsen—克努森扩散质量流量，kg/(m2·s)；
Jsurface—表面扩散质量流量，kg/(m2·s)；ρ—气体密度，kg/m3；μ—气体黏度，Pa·s；kD—页岩固有渗透率，m2；dm—气体分子直径，m；r—孔隙吼道半径，m；p—储层压力，MPa；pL—Langmuir 压力，MPa；F—滑脱系数，无量纲；Dk—克努森扩散系数，m2/s；M—气体摩尔质量，g/mol； Ds—表面扩散系数，m2/s；Csmax—吸附气最大吸附浓度，mol/m3；ε—贡献系数，无因次；R—气体常数，J/(mol·K)；pavg—平均压力，所述平均压力为在圆形单管中为进口、出口平均压力，Pa；α—切向动量调节系数，无因次，取值为0～1；T—储层温度，K；π—常数，
3.14。
[0192] 页岩储层中普遍具有一定的束缚水饱和度，不同页岩尺寸页岩毛细管半径对应的有效流动半径的计算公式如下：
[0193]
[0194] 式中：rie—不同尺寸毛细管对应的有效流动半径，mm；ri—不同尺寸毛细管对应的孔隙半径，mm；Sw—页岩孔隙束缚水饱和度，无量纲；i—计数符号，无量纲。
[0195] 不同孔径尺度下的毛细管表观渗透率计算公式如下：
[0196]
[0197] 其中，
[0198]
[0199]
[0200]
[0201]
[0202] 式中：kapp，i—不同尺寸毛细管的表观渗透率，m2；i—计数符号，无量纲；dm—气体分子直径，m；rie—不同尺寸毛细管对应的有效流动半径，m；p—储层压力，Pa；pL—Langmuir 压力，Pa；Fi—不同尺寸毛细管有效流动半径对应的滑脱系数，无量纲；εi—不同尺寸毛细管有效流动半径对应的贡献系数，无量纲；Kni—不同尺寸毛细管有效流动半径对应的克努森数，无量纲；Dki—不同尺寸毛细管有效流动半径对应的克努森扩散系数，m2/s；ρ—气体密度，kg/m3；μ—气体黏度，Pa·s；Dk—克努森扩散系数，m2/s；M—气体摩尔质量，g/mol； Ds—表面扩散系数，m2/s；Csmax—吸附气最大吸附浓度，mol/m3。
[0203] 根据单根毛细管的表观渗透率模型，以及不同尺寸毛细管所占的频率，最终得到页岩基质表观渗透率模型：
[0204]
[0205] 式中：kapp—整个岩心的表观渗透率，m2； —基质孔隙度，无量纲；kapp，i—不同尺度毛细管的表观渗透率，m2；λi—不同尺度毛细管的分布频率，无量纲；i，n—计数单位，无量纲；τ—迂曲度，无量纲。其中迂曲度的计算公式如下：
[0206]
[0207] 式中：m为迂曲度拟合参数，取值0.77。
[0208] 在本发明提供的实施例中，所述步骤S2中复杂裂缝网络渗透率通过微地震数据地震成像获取分形维数，通过蒙特卡洛随机建模理论建立二维离散裂缝网络模型，并进一步通过改进的立方定律来计算压裂施工后区域内复杂缝网的渗透率。
[0209] 在本发明提供的实施例中，所述步骤S2中分形维数通过盒子法计算，盒子法是用正方形的格子去覆盖复杂缝网区域，盒子的大小是变化的；给定盒子的尺码δ，可以数出覆盖住缝网区域所需的总盒子数目N；假设第i步覆盖使用δi×δi的格子，所需盒子数目为Ni(δi)，在第i+1步需要的格子(δi+1)×(δi+1)，则需要盒子数目Ni+1(δi+1)。得到任意两个尺度下所需盒子数之比与尺码之比的关系如下：
[0210]
[0211] 推广到：
[0212] N＝aδ-D (39)
[0213] 得到分形维数表达式为：
[0214]
[0215] 式中：D为分形维数，无量纲；N为盒子数目，无量纲；δ为盒子尺码，无量纲；i为计数符号。
[0216] 在本发明提供的实施例中，所述步骤S2的二维离散裂缝网络模型的建模方法为随机建模方法，所述二维离散裂缝网络模型的生成由以下步骤组成：
[0217] A、模拟区域生成：模拟区域根据页岩储层压裂过程中裂缝系统分布范围来确定；裂缝网络模拟时，先生成裂缝的中心点；考虑到有的裂缝只有部分在模拟区域内，中心点的位置可能在模拟区域外；为了减少边界效应的影响，可以适当应扩大模拟区域；为了便于区别，扩大后的模拟区域称为生成域，而原来的模拟区域称为分析区；为保证模拟结果可靠，生成域的侧面距分析域的距离大于裂缝可能出现的最大半径。
[0218] B、裂缝数目确定：模拟区域大小确定后，根据裂缝的面密度来确定该区域内每组裂缝的数目。
[0219] C、裂缝空间位置表征：对于二维离散裂缝网络模型，裂缝位置可以由中心点及迹长，或者裂缝的两个端点完全确定，即：
[0220]
[0221] D、裂缝几何参数：a按均匀分布生成裂缝中心点坐标随机变量x0、y0，b生成裂缝的倾向及倾角，c生成裂缝的迹长；d生成裂缝宽度，因裂缝宽度较难确定，通常把每组裂缝的宽度设为一常量。
[0222] E、裂缝网络模型生成及简化。当各组裂缝的数目、中心点坐标、产状、半径及宽度确定后，裂生成域内也就完全确定了，生成域的裂缝生成确定后，根据分析的几何边界来生成裂缝网络模型。
[0223] 在本发明提供的实施例中，所述步骤S5中复杂裂缝网络渗透率包括正常张开型裂缝渗透率、剪切型裂缝渗透率，计算渗透率过程中，正应力和剪应力分别进行计算，最后将正应力和剪应力诱导产生裂缝所具有渗透率进行叠加再进行复合，表达式如下：
[0224] kx＝knx+kdx (42)
[0225] ky＝kny+kdy (43)
[0226]
[0227] 式中：kx——在x方向上复杂缝网的渗透率，10-3μm2；knx——在x方向上由于裂缝正常闭合后的渗透率，10-3μm2；kdx——在x方向上由于裂缝发生剪切破坏产生的渗透率，10-3μm2； ky——在y方向上复杂缝网的渗透率，10-3μm2；kny——在y方向上由于裂缝正常闭合后的渗透率，10-3μm2；kdy——在y方向上由于裂缝发生剪切破坏产生的渗透率，10-3μm2；Kf—体积压裂区复杂裂缝网络渗透率，m2。
[0228] 在离散裂缝网络中，由于单条裂缝并不能反映出不同长度不同开度的裂缝，因此，在发明中使用等效裂缝频率(fx,fy)和等效开度(bx,by)的概念，使用这种形式来相对接近的表述离散裂缝网络中的复杂裂缝变形。通过等效频率和等效开度，正应力产生的张开或者闭合裂缝在x和y方向上(knx,kny)的渗透率可以通过下式计算所得：
[0229]
[0230]
[0231] 等效开度是指用一个值来近似代表整个离散裂缝网络所有裂缝的宽度，然后根据式(45) 的立方定律来表达离散裂缝网络的水力特性。其中裂缝的等效频率可以通过数值实验和初始开度反推得到：
[0232]
[0233]
[0234] 式中：bi为裂缝的初始平均开度，μm；
[0235] 裂缝的开度与地应力之间的关系，可以由如下公式表示：
[0236] b＝br+bm＝br+bmaxexp(-α'σ) (49)
[0237] 式中：br为裂缝的初始平均开度，μm；bmax为受压下裂缝最大开度，μm；α'为正常闭合下正常剪切下的函数斜率，无量纲；σ为正应力大小，MPa。
[0238] 根据式(49)可以从一维的公式推导到二维条件，在x方向上和在y方向上裂缝开度的大小：
[0239] bx＝br+bmax＝br+bmaxexp{-(αxσx+αyσy)} (50)
[0240] 式中：αx为x方向上应力系数，无量纲；αy为y方向上应力系数，无量纲；σx为x方向上正应力大小，MPa；σy为y方向上正应力大小，MPa。
[0241] 同理：
[0242] by＝br+bmax＝br+bmaxexp{-(βxσx+βyσy)} (51)
[0243] 式中：βx为正常闭合下正常剪切下的函数斜率，无量纲；βy为正常闭合下正常剪切下的函数斜率，无量纲；σx为x方向上正应力大小，MPa；σy为y方向上正应力大小，MPa。
[0244] 接下来计算
剪切应力条件下裂缝渗透率，计算公式如下：
[0245]
[0246]
[0247] 式中：fdx为在x方向上剪切裂缝的等效频率，1/m；fdy为在y方向上剪切裂缝的等效频率，1/m；dx为在x方向上剪切裂缝的等效开度，μm；dy为在y方向上剪切裂缝的等效开度，μm。
[0248] 由于只有在临界压力附近的裂缝才会剪切破坏，因此只有部分裂缝会出现剪切变形，其他大部分裂缝都不会出现剪切变形。所以发生剪切破坏的裂缝产生的渗透率会小于之前计算所得的正常闭合裂缝产生的渗透率。由裂缝剪切变形产生的开度大小可以通过下式计算：对于k
[0249] dx＝0 (54)
[0250] dy＝0 (55)
[0251] 对于k≥kc：
[0252] dx＝dmax[1-exp{-γx(k-kc)}] (56)
[0253] dy＝dmax[1-exp{-γy(k-kc)}] (57)
[0254] 式中：kc为裂缝变形的临界压力比，无量纲；dmax为裂缝发生剪切破坏后的最大开度，μm。
[0255] 导致裂缝破坏的临界压力比和临界方向值可以通过库伦破坏准则来计算：
[0256]
[0257]
[0258] 式中：φ为岩石内摩擦角，°；φf为岩石临界破坏角，°。
[0259] 最后，x和y方向上的渗透率(kx,ky)大小即可通过叠加正常闭合和剪切变形的裂缝渗透率得到，如下：
[0260]
[0261]
[0262] 整个裂缝网络区域的综合渗透率通过x与y方向渗透率均方根求得：
[0263]
[0264] 式中：kf—体积压裂区复杂裂缝网络渗透率，m2，kx为裂缝在x方向上的渗透率， -3 210 μm ；fx为在x方向上正应力产生的裂缝的等效频率，1/m；br为裂缝的初始平均开度，μm；
bmax为受压下裂缝最大开度，μm；αx为x方向上的压力系数，无量纲；σx为x方向上正应力大小，MPa；αy为y方向上的压力系数，无量纲；σy为y方向上正应力大小，MPa；fdx为在x方向上剪应力产生的裂缝的等效频率，1/m；dmax为裂缝发生剪切破坏后的最大开度，μm；γx为x方向剪切压力系数，无量纲；γy为y方向剪切压力系数，无量纲；k为x方向于 y方向上的应力比值，无量纲；kc为裂缝变形的临界压力比，无量纲；ky为裂缝在y方向上的渗透率，10-3μm2；fy为在y方向上正应力产生的裂缝的等效频率，1/m；βx为在xy下的函数斜率，无量纲；βy为在yx下的函数斜率，无量纲；fdy为在y方向上剪应力产生的裂缝等效频率，1/m。
[0265] 在本发明提供的实施例中，所述步骤S3中页岩气藏压裂水平井稳态产能模型将页岩气的渗流场分为内外两个流场，内部流场为体积压裂区的达西流动，外部流场为页岩基质中的克努森流动；
[0266] 其中外部流场的理想气体流动方程表达式如下：
[0267]
[0268] 其中
[0269] km＝kapp (64)
[0270]
[0271] 内部流场的流动方程如下：
[0272]
[0273] 式中：ν—流体流速，m/s；kf—体积压裂区复杂裂缝网络渗透率，m2；μ—气体黏度，mPa·s； —沿程压力梯度，MPa/m；
[0274] 页岩气藏压裂水平井的稳态产量方程表达式如下：
[0275]
[0276] 其中
[0277]
[0278]
[0279] 式中：qsc—页岩气标况下的产量，m3/d；pf—体积压裂区外边界压力，MPa；pwf—井底流压，MPa； —基质孔隙度，无量纲；pe—储层边界压力，MPa；kapp，i—不同尺度毛细管的表观渗透率，m2；μ—页岩气粘度，mPa·s；z—页岩气压缩因子，无量纲；psc—标准状况下的压力，MPa；rw—井眼半径，m；rf—体积压裂区的半径，m；re—供给半径，m；kf—体积压裂区复杂裂缝网络渗透率，m2；h—储层厚度，m；Tsc—标准状况下的温度，K；kapp—整个岩心的表观2
渗透率，m。
[0280] 在本发明提供的实施例中，所述步骤S3中页岩气藏压裂水平井非稳态模型把页岩气的流动分为体积压裂区和压裂区到储层边界页岩基质内外两个渗流场，页岩气的流动为径向流，同时考虑页岩气的解吸，得到页岩气非稳态下的产量；
[0281] 单位质量的页岩基质在压力发生改变时解吸出的页岩气体积可由Langmuir等温吸附解吸表达式得到，如下：
[0282]
[0283] 式中：Q—单位质量页岩基质解吸页岩气体积，m3/kg；pi—原始地层压力，MPa；
[0284] 在地面标准状况下，页岩气藏基质中的游离气量为：
[0285]
[0286] 式中：Gm—地下储集页岩气体积折算到地面标准状况下的体积，m3；Vm—页岩基质储集空间的体积，m3；Bgi—页岩气在原始地层压力和温度下的体积系数，无量纲；
[0287] 页岩气藏基质的体积为：
[0288]
[0289] 式中：Vt—页岩储层视体积，m3；Sw—束缚水饱和度，无量纲；
[0290] 根据朗格缪尔等温吸附理论，页岩气藏的吸附气体积为：
[0291]
[0292] Ga—页岩气基质中吸附气在地面标准状况下的体积，m3；ρs—页岩基质的密度，kg/m3；
[0293] 根据式(71)、(72)和(73)求得：
[0294]
[0295] 随着页岩气的不断被采出，地层压力下降，同时束缚水膨胀，基质孔隙收缩，裂缝中气体膨胀，令裂缝中的游离气为Gf。则储层中剩余游离气量为：
[0296] 基质中：
[0297] 改造区内：
[0298] 式中：G'm—基质剩余游离气在地面标准状况下的体积，m3；G'f—裂缝剩余游离气在地面标准状况下的体积，m3；Gf—裂缝中的游离气在标准状况下的体积，m3；Cm—页岩基质-1的压缩系数，MPa ；Bg—当前压力下的页岩气体积系数，无量纲；
[0299] 根据朗格缪尔等温吸附理论，可得到页岩气藏剩余的吸附气量：
[0300]
[0301] 根据物质平衡原理，原始地层压力下的基质中的游离气量、原始地层压力下裂缝中的游离气量、原始地层压力下的吸附气量之和等于采出气量、当前地层压力下基质中的游离气量、当前地层压力下裂缝中的游离气量、当前地层压力下的吸附气量之和，即：
[0302] Gm+Gf+Ga＝Gp+G'm+G'f+G′a (78)
[0303] 得到页岩气藏物质平衡方程表达式如下：
[0304]
[0305] 式中：Gm—地下储集页岩气体积折算到地面标准状况下的体积，m3；Gf—裂缝中的游离气在标准状况下的体积，m3；Gp—页岩气累计产出量，m3； —基质孔隙度，无量纲； Sw—束缚水饱和度，无量纲；Bgi—页岩气在原始地层压力和温度下的体积系数，无量纲；ρs—页岩基质的密度，kg/m3；pi—原始地层压力，MPa；VL为朗格缪尔压力，Pa；pL为朗格缪尔体积，m3/kg；
[0306] 其中：
[0307]
[0308] 页岩气压裂水平井生产过程中的压力传播半径方程表达式如下：
[0309]
[0310] 式中：R—t时刻的压力传播半径，m；K—储层渗透率，m2；t—生产时间，d； —基质孔隙度，无量纲；μ—气体黏度，mPa·s；Ct—储层的综合压缩系数，MPa-1；
[0311] 所述连续拟稳定法步骤如下：求得初始产量，假定一个时间步长，认为在该时间步长内的渗流为稳定渗流，计算出累计产量，根据物质平衡方程，算出在该时间段内压力波所波及到的区域的平均地层压力，以该平均地层压力作为下一时间步长的边界压力，计算下一时间步长下的产量，以此类推，便可得到产量与时间的关系，即非稳态下的产量。
[0312] 由于页岩气藏是致密储层，为获得商业气流往往需要大规模的体积压裂。所以，在页岩储层中，压力波传播的第一阶段又可以分为两个时段，即压力波先从井底传到体积压裂区的外边界，然后再从体积压裂区的外边界传到储层边界。压力波传播的距离即压力波及半径仅与时间、储层和流体的物性有关。因此，先设定一个时间步长Δt，可以得到压力波及半径R1，再根据稳态下产量公式式得到初始产量q1。由于页岩气藏生产初期所产气体是体积压裂区的游离气，所以时间步长Δt要取小，以免压力波及超出体积压裂区的外边界，造成较大误差。此后，再取相同时间步长Δtf，得到Δt+Δtf时刻压力波及R2，根据容积法可算出在当前压力波及半径内的地质储量，假定在Δtf时间内页岩气稳定产出，得到在Δtf时间内的累计产量，根据物质平衡方程计算当前压力波及半径内的地层平均压力p，然后根据稳态下页岩气产量公式计算Δt+Δtf时刻的产量q2，重复以上工作，得到页岩气非稳态的产量。
[0313] 实施例：
[0314] 某井场的页岩储层单根毛细管表观渗透率的
基础数据如表1所示：
[0315] 表1页岩气储层相关数据表
[0316]参数名 符号 单位 数值
气体类型 CH4 — —
-9
分子碰撞直径 δ m 0.42×10
气体分子直径 dm m 3.8×10-10
气体常数 R J/(mol·K) 8.314
温度 T K 423
气体摩尔质量 M g/mol 16
气体黏度 μ Pa·s 1.84×10-5
压缩因子 z 无量纲 0.85
切向动量调节系数 α 无量纲 0.8
气体分子密度 ρ kg/m3 0.655
原始地层压力,储层平均压力 P,pavg Pa 30×106
表面最大浓度 Csmax mol/m3 25040
Langmuir压力 pL Pa 2.46×10-6
表面扩散系数 Ds m2/s 2.89×10-10
[0317] 现场取回岩心分析得到基质孔隙度、束缚水饱和度及毛细管管径分布结果见表2。
[0318] 表2页岩岩心基质孔隙度及毛细管
孔径分布[0319]
[0320] 根据页岩气井压裂现场微地震数据反演以及实验测试获得的分形维数及缝网基础参数结果见表3。
[0321] 表3缝网区域基础参数
[0322]参数 单位 数值
分形维数 - 1.47
最大裂缝长度rmax m 150
最小裂缝长度rmin m 3
裂缝区域宽度 m 100
裂缝区域长度 m 100
裂缝面密度 条/m2 3
bx上的应力系数(αx，αy)* - 0.034，0.146
*
by上的应力系数(βx，βy) - 0.101，0.071
剩余开度(br) μm 50
最大开度(bmax) μm 300
最大剪切开度(dmax) μm 500
裂缝初始开度(bi) μm 300
摩擦角* ° 24.9
临界压力比 无因次 2.45
有效水平主应力(σx) MPa 37.8～65
有效水平主应力(σy) MPa 32.5～67
[0323] 计算页岩气压裂水平井非稳态产量需要的基础数据见表4。
[0324] 表4地层基础参数
[0325]
[0326] 根据以上参数，通过式(36)计算得到整个岩心的表观渗透率，如图1所示，图1表示页岩基质表观渗透率计算结果。
[0327] 根据以上表3数据，得到页岩压裂水平井微地震反演裂缝分布图像，如图2所示；根据微地震反演出的裂缝分布图像进行灰度处理，识别出该区域内的裂缝分布，如图3所示；通过盒子法计算出裂缝区域的分形维数图，如图4所示。
[0328] 根据以上数据，通过式(62)计算得到复杂裂缝网络的渗透率，分析得到x方向与y方向水平主应力相同时渗透率的变化情况，如图5所示；x方向与y方向水平主应力不同时渗透率的变化情况，如图6所示。
[0329] 根据以上数据，通过连续拟稳定法得到页岩气藏压裂水平井复杂裂缝网络的非稳态产量，且本发明预测结果与Eclipse模型对比结果如图7所示，从图7可以看出本发明计算结果与该井的实际生产情况吻合，且误差较小；通过对计算结果分析，得到不同分形维数下页岩气压裂水平井累产的变化情况，如图8所示；得到不同压裂半径下，页岩气压裂水平井累产的变化情况，如图9所示；得到不同裂缝开度下，页岩气压裂水平井累产的变化情况，如图10所示。
[0330] 以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单
修改、等同变化与修饰，均仍属于本发明技术方案的范围内。
