CN108518212B - 一种计算页岩气藏复杂裂缝网络非稳态产量的方法 - Google Patents
一种计算页岩气藏复杂裂缝网络非稳态产量的方法 Download PDFInfo
- Publication number
- CN108518212B CN108518212B CN201810311636.6A CN201810311636A CN108518212B CN 108518212 B CN108518212 B CN 108518212B CN 201810311636 A CN201810311636 A CN 201810311636A CN 108518212 B CN108518212 B CN 108518212B
- Authority
- CN
- China
- Prior art keywords
- permeability
- shale
- shale gas
- gas
- reservoir
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- 230000035699 permeability Effects 0.000 claims abstract description 94
- 239000011159 matrix material Substances 0.000 claims abstract description 68
- 238000004364 calculation method Methods 0.000 claims abstract description 26
- 238000003384 imaging method Methods 0.000 claims abstract description 11
- 239000000463 material Substances 0.000 claims abstract description 11
- 238000011105 stabilization Methods 0.000 claims abstract description 8
- 230000004856 capillary permeability Effects 0.000 claims abstract description 5
- 206010017076 Fracture Diseases 0.000 claims description 137
- 208000010392 Bone Fractures Diseases 0.000 claims description 135
- 230000015572 biosynthetic process Effects 0.000 claims description 38
- 238000009792 diffusion process Methods 0.000 claims description 35
- 230000005012 migration Effects 0.000 claims description 29
- 238000013508 migration Methods 0.000 claims description 29
- 238000004519 manufacturing process Methods 0.000 claims description 28
- 239000011148 porous material Substances 0.000 claims description 28
- 238000001179 sorption measurement Methods 0.000 claims description 26
- 238000009826 distribution Methods 0.000 claims description 16
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 11
- 239000011435 rock Substances 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 7
- 238000010008 shearing Methods 0.000 claims description 6
- 230000007246 mechanism Effects 0.000 claims description 5
- 208000002565 Open Fractures Diseases 0.000 claims description 3
- 230000004308 accommodation Effects 0.000 claims description 3
- 230000000877 morphologic effect Effects 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims description 2
- 230000009286 beneficial effect Effects 0.000 abstract description 2
- 229910052500 inorganic mineral Inorganic materials 0.000 abstract description 2
- 239000011707 mineral Substances 0.000 abstract description 2
- 239000007789 gas Substances 0.000 description 256
- 238000004088 simulation Methods 0.000 description 19
- 230000000694 effects Effects 0.000 description 14
- 238000003795 desorption Methods 0.000 description 11
- 238000004458 analytical method Methods 0.000 description 6
- 230000008859 change Effects 0.000 description 6
- 230000001186 cumulative effect Effects 0.000 description 6
- 238000013461 design Methods 0.000 description 5
- 230000002829 reductive effect Effects 0.000 description 5
- 239000012530 fluid Substances 0.000 description 4
- 230000005540 biological transmission Effects 0.000 description 3
- 238000007906 compression Methods 0.000 description 3
- 230000007423 decrease Effects 0.000 description 3
- 238000002474 experimental method Methods 0.000 description 3
- 230000009471 action Effects 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 2
- 230000006835 compression Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 230000000977 initiatory effect Effects 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 230000033001 locomotion Effects 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 238000009827 uniform distribution Methods 0.000 description 2
- 208000003044 Closed Fractures Diseases 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005553 drilling Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 238000012067 mathematical method Methods 0.000 description 1
- 239000002071 nanotube Substances 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH OR ROCK DRILLING; MINING
- E21B—EARTH OR ROCK DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B43/00—Methods or apparatus for obtaining oil, gas, water, soluble or meltable materials or a slurry of minerals from wells
- E21B43/25—Methods for stimulating production
- E21B43/26—Methods for stimulating production by forming crevices or fractures
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
- G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
- G06Q50/02—Agriculture; Fishing; Forestry; Mining
Landscapes
- Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Mining & Mineral Resources (AREA)
- Geology (AREA)
- Physics & Mathematics (AREA)
- Business, Economics & Management (AREA)
- Environmental & Geological Engineering (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geochemistry & Mineralogy (AREA)
- Fluid Mechanics (AREA)
- Human Resources & Organizations (AREA)
- Marketing (AREA)
- Animal Husbandry (AREA)
- Health & Medical Sciences (AREA)
- Economics (AREA)
- General Health & Medical Sciences (AREA)
- Agronomy & Crop Science (AREA)
- Marine Sciences & Fisheries (AREA)
- Primary Health Care (AREA)
- Strategic Management (AREA)
- Tourism & Hospitality (AREA)
- General Business, Economics & Management (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geophysics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明属于矿物资源技术领域，具体而言，涉及一种计算页岩气藏复杂裂缝网络非稳态产量的方法，包括以下步骤：S1、收集页岩气藏地层、页岩气基本参数，建立页岩基质多孔介质毛细管渗透率模型；S2、收集页岩气藏微地震裂缝成像数据，运用盒子法获取分形维数，建立二维离散裂缝网络模型，利用改进的立方定律计算复杂裂缝网络渗透率；S3、考虑页岩气藏压裂水平井外区为基质区域渗流区域、内区为复杂缝网区域流动区域，建立页岩气藏压裂水平井稳态产能模型；根据页岩气压力传播半径方程、页岩气藏物质平衡方程，运用连续拟稳定法建立页岩气藏压裂水平井非稳态模型。本发明的有益效果：考虑因素贴近地层实际情况、计算结果更加精确。
Description
技术领域
本发明涉及一种计算页岩气藏复杂裂缝网络非稳态产量的方法，属于矿物资源技术领域。
背景技术
页岩气资源潜力巨大，近年来受到了人们的广泛关注。随着钻井以及压裂工艺的迅速提高，页岩气已经成为油气勘探开发的重要组成部分。水平井压裂形成复杂缝网是开发页岩气的关键技术。页岩储层经过压裂施工后会在压裂改造区域会形成复杂裂缝网络，而在未改造的区域仍然为基质孔隙区域。页岩气压裂水平井在生产过程中，根据页岩气的流动特征，可以把压裂后的页岩气井渗流场分为内外两个流场，内部流场为体积压裂形成复杂缝网区域的达西流动，而在外部流场为页岩基质区域受到黏性流、扩散、解吸和滑脱等多重运移机制组成；并且页岩气压裂水平井的生产是一个不稳定渗流过程，而且初期产量递减十分明显，只有到了生产后期压力波传到边界时，产量才趋于稳定。并且在生产过程中，储层压力下降还会显著影响页岩气体在基质和缝网中的流动。
目前常用的计算页岩气压裂水平井产量的源函数方法一开始就假设由于生产导致的压力传播一开始就传递到储层边界，这与实际情况不吻合；且计算的页岩基质渗透率建立的是单毛细管模型，与地层实际情况不符；复杂裂缝网络渗透率通过实验测得，不能体现实际复杂裂缝网络实际情况，从而使得计算结果与实际生产误差较大。
发明内容
针对现有技术中的存在的不足，本发明提供了更加贴近地层实际情况、计算结果更加精确的一种计算页岩气藏复杂裂缝网络非稳态产量的方法。
本发明的技术方案如下：
一种计算页岩气藏复杂裂缝网络非稳态产量的方法，其特征在于，包括以下步骤：
S1、收集页岩气藏地层、页岩气的基本参数，建立页岩基质多孔介质毛细管渗透率模型；
S2、收集页岩气藏微地震裂缝成像数据，根据微地震裂缝成像数据，识别出裂缝条数、位置及形态特征；运用盒子法获取表征岩石复杂裂缝网络的分形维数；根据求取的分形维数，利用蒙特卡洛随机建模理论生成二维离散裂缝网络模型；通过立方定律计算复杂裂缝网络渗透率；
S3、假设储层等厚、圆形边界、半径为re，体积压裂区域改造充分、圆形边界、半径为rf、与储层厚度相同；考虑页岩气藏压裂水平井外区为基质区域渗流区域、内区为复杂缝网区域流动区域，建立页岩气藏压裂水平井稳态产能模型；建立页岩气压力传播半径方程、页岩气藏物质平衡方程，运用连续拟稳定法建立页岩气藏压裂水平井非稳态模型，求解页岩气藏压裂水平井非稳态产量。
进一步的，所述页岩气藏地层、页岩气的基本参数包括储层边界压力、储层温度、供给半径、基质孔隙度、裂缝孔隙度、体积压裂区半径、井眼半径、储层厚度、页岩气粘度、气体类型、分子碰撞直径、气体分子直径、气体常数、气体摩尔质量、切向动量调节系数、气体分子密度、平均压力、表面最大浓度、朗格缪尔压力和表面扩散系数。
进一步的，所述步骤S1中页岩基质多孔介质毛细管渗透率模型是根据毛细管大小的分布不同分别计算不同尺度下毛细管的渗透率，然后对不同尺度下的毛细管渗透率进行叠加，从而得到整个岩心的渗透率大小，其表达式如下：
式中：kapp—整个岩心的表观渗透率，m2；
—基质孔隙度，无量纲；kapp，i—不同尺度毛细管的表观渗透率，m2；λi—不同尺度毛细管的分布频率，无量纲；i，n—计数单位，无量纲；τ—迂曲度，无量纲。其中迂曲度计算公式如下：
式中：m为迂曲度拟合参数，取值0.77。
进一步的，所述步骤S1中页岩基质多孔介质毛细管渗透率模型，根据达西定律求得页岩气藏整个岩心的表观渗透率的步骤是利用页岩气储层参数计算得到克努森系数，根据克努森系数判断毛细管中的流态，利用毛细管中的流态和运移机制建立气体质量运移方程，接下来建立统一的质量运移方程，利用统一的质量运移方程得到毛细管的表观渗透率，最后进行叠加求得整个岩心的表观渗透率；
其中克努森系数的计算公式如下：
式中：Kn—克努森系数，无因次；kB—玻尔兹曼常数，1.3805×10-23J/K；T—页岩储层温度，K；π—常数，3.14；δ—气体分子碰撞直径，m；p—储层压力，Pa；d—毛细管直径，m。
进一步的，所述步骤S2中复杂裂缝网络渗透率通过微地震数据地震成像获取分形维数，通过蒙特卡洛随机建模理论建立二维离散裂缝网络模型，并进一步通过立方定律来计算压裂施工后区域内复杂缝网的渗透率，其表达式如下：
式中：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。
进一步的，所述步骤S2中分形维数通过盒子法计算，盒子法是用正方形的格子去覆盖复杂缝网区域，盒子的大小是变化的；给定盒子的尺码δ，可以数出覆盖住缝网区域所需的总盒子数目N；假设第i步覆盖使用δi×δi的格子，所需盒子数目为Ni(δi)，在第i+1步需要的格子(δi+1)×(δi+1)，则需要盒子数目Ni+1(δi+1)。得到任意两个尺度下所需盒子数之比与尺码之比的关系如下：
推广到：
N＝aδ-D
得到分形维数表达式为：
式中：D为分形维数，无量纲；N为盒子数目，无量纲；δ为盒子尺码，无量纲；i为计数符号；
所述步骤S2的二维离散裂缝网络模型的建模方法为随机建模方法，所述二维离散裂缝网络模型的生成由以下步骤组成：
A、模拟区域生成：模拟区域根据页岩储层压裂过程中裂缝系统分布范围来确定；裂缝网络模拟时，先生成裂缝的中心点；考虑到有裂缝只有部分在模拟区域内，中心点的位置可能在模拟区域外，为了减少边界效应的影响，可以适当扩大模拟区域；为了便于区别，扩大后的模拟区域称为生成域，而原来的模拟区域称为分析区；为保证模拟结果可靠，生成域的每一个侧面距分析域的距离大于裂缝出现的最大半径；
B、裂缝数目确定：模拟区域大小确定后，根据裂缝的面密度来确定该区域内每组裂缝的数目；
C、裂缝空间位置表征：对于二维离散裂缝网络模型，裂缝位置可以由中心点及迹长，或者裂缝的两个端点完全确定，即：
D、裂缝几何参数：a按均匀分布生成裂缝中心点坐标随机变量x0、y0；b生成裂缝的倾向及倾角；c生成裂缝的迹长；d生成裂缝宽度，因裂缝宽度较难确定，通常把每组裂缝的宽度设为一常量；
E、裂缝网络模型生成及简化：当各组裂缝的数目、中心点坐标、产状、半径及宽度确定后，裂生成域内也就完全确定了；生成域的裂缝生成确定后，根据分析的几何边界来生成裂缝网络模型。
进一步的，所述步骤S2中复杂裂缝网络渗透率包括正常张开型裂缝渗透率、剪切型裂缝渗透率，计算渗透率过程中，正应力和剪应力分别进行计算，最后将正应力和剪应力诱导产生裂缝所具有渗透率进行叠加再复合得到，表达式如下：
kx＝knx+kdx
ky＝kny+kdy
式中：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。
进一步的，所述步骤S3中页岩气藏压裂水平井稳态产能模型将页岩气的渗流场分为内外两个流场，内部流场为体积压裂区的达西流动，外部流场为页岩基质中的多尺度毛细管流动；
其中外部流场的气体流动方程表达式如下：
其中
km＝kapp
内部流场的流动方程如下：
式中：ν—流体流速，m/s；Kf—体积压裂区复杂裂缝网络渗透率，m2；μ—气体黏度，mPa·s；
—沿程压力梯度，MPa/m；
页岩气藏压裂水平井的稳态产量方程表达式如下：
其中
式中：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。
进一步的，所述步骤S3中页岩气藏压裂水平井非稳态模型把页岩气的流动分为体积压裂区和压裂区到储层边界页岩基质内外两个渗流场，页岩气的流动为径向流，同时考虑页岩气的解吸，得到页岩气非稳态下的产量；
单位质量的页岩基质在压力发生改变时解吸出的页岩气体积可由Langmuir等温吸附解吸表达式得到，如下：
式中：Q—单位质量页岩基质解吸页岩气体积，m3/kg；pi—原始地层压力，MPa；
在地面标准状况下，页岩气藏基质中的游离气量为：
式中：Gm—地下储集页岩气体积折算到地面标准状况下的体积，m3；Vm—页岩基质储集空间的体积，m3；Bgi—页岩气在原始地层压力和温度下的体积系数，无量纲；
页岩气藏基质的体积为：
式中：Vt—页岩储层视体积，m3；Sw—束缚水饱和度，无量纲；
根据朗格缪尔等温吸附理论，页岩气藏的吸附气体积为：
Ga—页岩气基质中吸附气在地面标准状况下的体积，m3；ρs—页岩基质的密度，kg/m3；
根据上式，进一步求得：
随着页岩气的不断被采出，地层压力下降，同时束缚水膨胀，基质孔隙收缩，裂缝中气体膨胀，令裂缝中的游离气为Gf。则储层中剩余游离气量为：
基质中：
改造区内：
式中：G'm—基质剩余游离气在地面标准状况下的体积，m3；G'f—裂缝剩余游离气在地面标准状况下的体积，m3；Gf—裂缝中的游离气在标准状况下的体积，m3；Cm—页岩基质的压缩系数，MPa-1；Bg—当前压力下的页岩气体积系数，无量纲；
根据朗格缪尔等温吸附理论，可得到页岩气藏剩余的吸附气量：
根据物质平衡原理，原始地层压力下的基质中的游离气量、原始地层压力下裂缝中的游离气量、原始地层压力下的吸附气量之和等于采出气量、当前地层压力下基质中的游离气量、当前地层压力下裂缝中的游离气量、当前地层压力下的吸附气量之和，即：
Gm+Gf+Ga＝Gp+G'm+G'f+G′a
得到页岩气藏物质平衡方程表达式如下：
式中：Gm—地下储集页岩气体积折算到地面标准状况下的体积，m3；Gf—裂缝中的游离气在标准状况下的体积，m3；Gp—页岩气累计产出量，m3；
—基质孔隙度，无量纲；Sw—束缚水饱和度，无量纲；Bgi—页岩气在原始地层压力和温度下的体积系数，无量纲；ρs—页岩基质的密度，kg/m3；pi—原始地层压力，MPa；VL为朗格缪尔压力，Pa；pL为朗格缪尔体积，m3/kg；
其中：
页岩气压裂水平井生产过程中的压力传播半径方程表达式如下：
式中：R—t时刻的压力传播半径，m；K—储层渗透率，m2；t—生产时间，d；
—基质孔隙度，无量纲；μ—气体黏度，mPa·s；Ct—储层的综合压缩系数，MPa-1。
所述连续拟稳定法步骤如下：求得初始产量，假定一个时间步长，认为在该时间步长内的渗流为稳定渗流，计算出累计产量，根据物质平衡方程，算出在该时间段内压力波所波及到的区域的平均地层压力，以该平均地层压力作为下一时间步长的边界压力，计算下一时间步长下的产量，以此类推，便可得到产量与时间的关系，即非稳态下的产量。
本发明的有益效果为：
本发明方法充分考虑地层岩芯由多尺度毛细管组成，首先求得计算不同尺度下毛细管的渗透率，然后对不同尺度下的毛细管渗透率进行叠加，得到整个岩心的渗透率更加贴近实际情况；通过页岩气藏压裂水平井微地震裂缝成像数据，运用盒子法获取分形维数，建立二维离散裂缝网络模型，利用改进的立方定律计算复杂裂缝网络渗透率；克服源函数等方法的缺陷，利用连续拟稳定法，考虑压力波在生产过程中逐渐往外传播，压力波的传播速度与储层物性、流体性质和生产时间紧密相关，分析了在页岩气压裂水平井生产过程中在过渡流和边界流条件下页岩储层压力的变化规律；并且使用连续拟稳定法计算非稳态产量形式简单，容易实现程序化，同时避免了源函数等数学方法再求解产量时在生产初期压力波就传递到储层边界的不合理假设。
附图说明
为了更清楚地说明本发明实施方式的技术方案，下面将对实施方式中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它相关的附图。
图1本发明提供的页岩基质表观渗透率；
图2本发明提供的页岩压裂水平井微地震反演裂缝分布图；
图3本发明提供的裂缝分布图；
图4本发明提供的分形维数计算图；
图5本发明提供的裂缝网络渗透率与裂缝区域平均应力的关系图；
图6本发明提供的裂缝网络渗透率与x与y方向应力的关系图；
图7本发明提供的产量计算结果与Eclipse模拟结果对比图；
图8本发明提供的不同分形维数下页岩气压裂水平井累产的曲线图；
图9本发明提供的不同压裂半径下页岩气压裂水平井累产的曲线图；
图10本发明提供的不同裂缝开度下，页岩气压裂水平井累产的曲线图。
具体实施方式
下面结合附图和实施例对本发明进一步说明。
一种计算页岩气藏复杂裂缝网络非稳态产量的方法，包括以下步骤：
S1、收集页岩气藏地层、页岩气的基本参数，建立页岩基质多孔介质毛细管渗透率模型；
S2、收集页岩气藏微地震裂缝成像数据，根据微地震裂缝成像数据，识别出裂缝条数、位置及形态特征；运用盒子法获取表征岩石复杂裂缝网络的分形维数；根据求取的分形维数，利用蒙特卡洛随机建模理论生成二维离散裂缝网络模型；通过立方定律计算复杂裂缝网络渗透率；
S3、假设储层等厚、圆形边界、半径为re，体积压裂区域改造充分、圆形边界、半径为rf、与储层厚度相同；考虑页岩气藏压裂水平井外区为基质区域渗流区域、内区为复杂缝网区域流动区域，建立页岩气藏压裂水平井稳态产能模型；根据页岩气压力传播半径方程、页岩气藏物质平衡方程，运用连续拟稳定法建立页岩气藏压裂水平井非稳态模型。
在本发明提供的实施例中，所述页岩气藏地层、页岩气的基本参数包括储层边界压力、储层温度、供给半径、基质孔隙度、裂缝孔隙度、体积压裂区半径、井眼半径、储层厚度、页岩气粘度、气体类型、分子碰撞直径、气体分子直径、气体常数、气体摩尔质量、切向动量调节系数、气体分子密度、平均压力、表面最大浓度、朗格缪尔压力和表面扩散系数。
在本发明提供的实施例中，所述步骤S1中页岩基质多孔介质毛细管渗透率模型是根据毛细管大小的分布不同分别计算不同尺度下毛细管的渗透率，然后对不同尺度下的毛细管渗透率进行叠加，从而得到整个岩心的渗透率大小，其表达式如下：
式中：kapp—整个岩心的表观渗透率，m2；
—基质孔隙度，无量纲；kapp，i—不同尺度毛细管的表观渗透率，m2；λi—不同尺度毛细管的分布频率，无量纲；i，n—计数单位，无量纲；τ—迂曲度，无量纲。其中迂曲度计算公式如下：
式中：m为迂曲度拟合参数，取值0.77。
在本发明提供的实施例中，所述步骤S1中页岩基质多孔介质毛细管渗透率模型，根据达西定律求得页岩气藏整个岩心的表观渗透率的步骤是利用页岩气储层参数计算得到克努森系数，根据克努森系数判断毛细管中的流态，利用毛细管中的流态和运移机制建立气体质量运移方程，接下来建立统一的质量运移方程，利用统一的质量运移方程得到毛细管的表观渗透率，最后进行叠加求得整个岩心的表观渗透率；
其中克努森系数的计算公式如下：
式中：Kn—克努森系数，无因次；kB—玻尔兹曼常数，1.3805×10-23J/K；p—储层压力，Pa；T—页岩储层温度，K；π—常数，3.14；δ—气体分子碰撞直径，m；d—毛细管直径，m；
其中运移机制包括游离气运移、吸附气运移，所述游离气运移包括黏性流动、滑脱流动、Knudsen扩散，所述吸附气运移包括吸附气解吸、吸附气表面扩散；
(1)游离态页岩气质量运移方程
游离态页岩气(游离气)赋存于基质孔隙以及裂缝中，主要包括黏性流动、滑脱流动、Knudsen扩散作用。
①粘性流质量运移方程
页岩发育有大量的纳米级孔隙，可将纳米孔视为毛管模型，页岩则可视为毛管和基质组成。当页岩气体Knudsen数远小于1时，气体分子的运动主要受分子间碰撞支配，此时分子与壁面的碰撞较少，气体分子间的相互作用要比气体分子与孔隙表面(孔隙壁)的碰撞频繁得多，气体以连续流动为主，可采用黏性流质量运移方程描述。当不考虑吸附气存在对毛管半径的影响时，对于吼道半径为r的单根毛管，其固有渗透率计算公式如下：
式中：kD—固有渗透率，m2；r—页岩毛细管半径，m，r＝d/2。
在单组分气体之间存在压力梯度所引起的黏性流动，可以用达西定律来表示描述黏性流的质量运移方程，表达式如下：
式中：Jvicious—黏性流质量流量，kg/(m2·s)；ρ—气体密度，kg/m3；μ—气体黏度，Pa·s；p—储层压力，Pa；
—压力梯度算子符号，无量纲。
对于气体在纳米管中的运移，当考虑吸附态气体存在对纳米孔半径的影响时，纳米孔吼有效半径减小，因此考虑吸附气影响时纳米孔吼的有效半径可表达为：
式中：re—考虑吸附气影响的孔隙有效半径，m；dm—气体分子直径，m；pL—Langmuir压力，Pa。
将式(6)代入式(5)，可以得到：
②滑脱效应质量运移方程
当页岩孔隙尺度减小，或者气体压力降低、气体分子自由程增加，气体分子自由程与孔隙直径的尺度具有可比性，气体分子与孔隙壁面的碰撞不可忽略。在0.001＜Kn＜0.1时，由于壁面页岩气分子速度不再为零，此时存在滑脱现象，且考虑滑脱效应时的页岩储层渗透率计算公式如下：
式中：kslip—考虑滑脱效应渗透率，m2；paver—通过实验测试岩心渗透率时的进出口平均压力，Pa；bk—滑脱因子，Pa。
为了能将滑脱效应在渗流方程中体现，引入滑脱因子来修正纳米孔隙滑脱效应，将式(9)代入式(8)：
式中：R—气体常数，J/(mol·K)；M—气体摩尔质量，kg/mol；pavg—平均压力(在圆形单管中为进口、出口平均压力)，Pa；α—切向动量调节系数，无因次，取值为0～1。
因此考虑滑脱效应时的质量运移方程可以写为：
式中：Jslip—滑脱流质量流量，kg/(m2·s)。
③克努森扩散质量运移方程
当孔道直径减少或者分子平均自由程增加(在低压下)，Kn＞10时，气体分子更容易与孔隙壁面发生碰撞而不是与其他气体分子发生碰撞，这意味着气体分子达到了几乎能独立于彼此的点，称为Knudsen扩散。
Jknudsen＝α0νρ (12)
式中：α0—无因次概率系数，无因次；ν—平均分子速度，m/s；ρ—气体分子密度，kg/m3。
当圆管两端都有气体，圆管传输的净流量与圆管两端的气体密度成正比，式(12)可写为：
Jknudsen＝αν(ρin-ρout) (13)
式中：ρin——圆管进口处气体密度，kg/m3；ρout——圆管出口处气体密度，kg/m3。
根据气体动力学理论，气体的平均分子运动速度为：
对于直径为d长度为L的圆形长直管(L>>d)，α＝d/3L，将式(14)代入到式(13)中，可以得到：
将式(15)写为偏微分形式：
式(16)也可写为气体浓度C的形式，即为：
而纳米孔隙中的Knudsen扩散系数Dk，表达式如下所示：
式中：Dk—克努森扩散系数，m2/s。
而气体密度表达形式可写为：
将式(18)、(17)代入式(15)，可得：
因此，Knudsen扩散质量运移方程可表述为：
式中：Jknudsen—克努森扩散质量流量，kg/(m2·s)。
(2)吸附态页岩气运移质量方程
吸附态页岩气(吸附气)赋存于孔隙壁面和页岩固体颗粒表面，主要发生解吸附作用和表面扩散作用。假设在一定温度和压力条件下，壁面吸附气和自由气处于瞬间动态平衡，采用Langmuir等温吸附模型，吸附质量的表达形式为：
式中：qads—页岩单位体积的吸附量，kg/m3；Vstd—页岩气标况下摩尔体积，m3/mol；
页岩气井在生产过程中，地层压力逐渐下降，考虑t1时刻地层压力为p1，t2时刻地层压力为p2，则可计算出地层压力由p1下降为p2时吸附态页岩气的解吸量：
式中：Δqads—由于压力下降产生的吸附态页岩气解析量，kg/m3；VL—朗格缪尔体积，m3/kg。
页岩气在微纳米孔隙表面不仅存在解吸附效应，还存在沿吸附壁面的运移，即表面扩散作用。不同于压力梯度或浓度梯度作用的其他运移方式，页岩气表面扩散在吸附势场的作用下发生运移，影响页岩气表面扩散的因素很多，包括压力、温度、纳米孔壁面属性、页岩气体分子属性、页岩气体分子与纳米孔壁面相互作用等。
当表面扩散气体运移方程表达为浓度梯度的形式时，等于表面扩散系数与浓度梯度的乘积形式，其表面扩散质量流量计算公式如下所示：
式中：Jsurface—表面扩散质量流量，kg/(m2·s)；Cs—孔壁面吸附气浓度，mol/m3；Ds—表面扩散系数，m2/s；l—孔壁长度，m。
吸附气体覆盖率θ可表示为：
式中：θ—吸附气体覆盖率，无因次；Csmax—吸附气最大吸附浓度，mol/m3；V—单位质量页岩实际吸附气体积，m3/kg；Csmax—吸附气最大吸附浓度，mol/m。
式(25)可进一步改写为：
将式(26)代入式(24)，可得到满足Langmuir等温吸附方程的页岩气表面扩散质量运移方程：
在本发明提供的实施例中，考虑游离态页岩气粘性流、滑脱流、Knudsen扩散和吸附态页岩气的解吸、表面扩散作用，其总的传输质量为这几种运移模式引起的传输质量的叠加之和，并引入贡献系数ε，建立统一的的质量运移方程，式(28)：
其中，F为滑脱系数，其计算公式如下：
式中：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。
页岩储层中普遍具有一定的束缚水饱和度，不同页岩尺寸页岩毛细管半径对应的有效流动半径的计算公式如下：
式中：rie—不同尺寸毛细管对应的有效流动半径，mm；ri—不同尺寸毛细管对应的孔隙半径，mm；Sw—页岩孔隙束缚水饱和度，无量纲；i—计数符号，无量纲。
不同孔径尺度下的毛细管表观渗透率计算公式如下：
其中，
式中：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。
根据单根毛细管的表观渗透率模型，以及不同尺寸毛细管所占的频率，最终得到页岩基质表观渗透率模型：
式中：kapp—整个岩心的表观渗透率，m2；
—基质孔隙度，无量纲；kapp，i—不同尺度毛细管的表观渗透率，m2；λi—不同尺度毛细管的分布频率，无量纲；i，n—计数单位，无量纲；τ—迂曲度，无量纲。其中迂曲度的计算公式如下：
式中：m为迂曲度拟合参数，取值0.77。
在本发明提供的实施例中，所述步骤S2中复杂裂缝网络渗透率通过微地震数据地震成像获取分形维数，通过蒙特卡洛随机建模理论建立二维离散裂缝网络模型，并进一步通过改进的立方定律来计算压裂施工后区域内复杂缝网的渗透率。
在本发明提供的实施例中，所述步骤S2中分形维数通过盒子法计算，盒子法是用正方形的格子去覆盖复杂缝网区域，盒子的大小是变化的；给定盒子的尺码δ，可以数出覆盖住缝网区域所需的总盒子数目N；假设第i步覆盖使用δi×δi的格子，所需盒子数目为Ni(δi)，在第i+1步需要的格子(δi+1)×(δi+1)，则需要盒子数目Ni+1(δi+1)。得到任意两个尺度下所需盒子数之比与尺码之比的关系如下：
推广到：
N＝aδ-D (39)
得到分形维数表达式为：
式中：D为分形维数，无量纲；N为盒子数目，无量纲；δ为盒子尺码，无量纲；i为计数符号。
在本发明提供的实施例中，所述步骤S2的二维离散裂缝网络模型的建模方法为随机建模方法，所述二维离散裂缝网络模型的生成由以下步骤组成：
A、模拟区域生成：模拟区域根据页岩储层压裂过程中裂缝系统分布范围来确定；裂缝网络模拟时，先生成裂缝的中心点；考虑到有的裂缝只有部分在模拟区域内，中心点的位置可能在模拟区域外；为了减少边界效应的影响，可以适当应扩大模拟区域；为了便于区别，扩大后的模拟区域称为生成域，而原来的模拟区域称为分析区；为保证模拟结果可靠，生成域的侧面距分析域的距离大于裂缝可能出现的最大半径。
B、裂缝数目确定：模拟区域大小确定后，根据裂缝的面密度来确定该区域内每组裂缝的数目。
C、裂缝空间位置表征：对于二维离散裂缝网络模型，裂缝位置可以由中心点及迹长，或者裂缝的两个端点完全确定，即：
D、裂缝几何参数：a按均匀分布生成裂缝中心点坐标随机变量x0、y0，b生成裂缝的倾向及倾角，c生成裂缝的迹长；d生成裂缝宽度，因裂缝宽度较难确定，通常把每组裂缝的宽度设为一常量。
E、裂缝网络模型生成及简化。当各组裂缝的数目、中心点坐标、产状、半径及宽度确定后，裂生成域内也就完全确定了，生成域的裂缝生成确定后，根据分析的几何边界来生成裂缝网络模型。
在本发明提供的实施例中，所述步骤S5中复杂裂缝网络渗透率包括正常张开型裂缝渗透率、剪切型裂缝渗透率，计算渗透率过程中，正应力和剪应力分别进行计算，最后将正应力和剪应力诱导产生裂缝所具有渗透率进行叠加再进行复合，表达式如下：
kx＝knx+kdx (42)
ky＝kny+kdy (43)
式中：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。
在离散裂缝网络中，由于单条裂缝并不能反映出不同长度不同开度的裂缝，因此，在发明中使用等效裂缝频率(fx,fy)和等效开度(bx,by)的概念，使用这种形式来相对接近的表述离散裂缝网络中的复杂裂缝变形。通过等效频率和等效开度，正应力产生的张开或者闭合裂缝在x和y方向上(knx,kny)的渗透率可以通过下式计算所得：
等效开度是指用一个值来近似代表整个离散裂缝网络所有裂缝的宽度，然后根据式(45)的立方定律来表达离散裂缝网络的水力特性。其中裂缝的等效频率可以通过数值实验和初始开度反推得到：
式中：bi为裂缝的初始平均开度，μm；
裂缝的开度与地应力之间的关系，可以由如下公式表示：
b＝br+bm＝br+bmaxexp(-α'σ) (49)
式中：br为裂缝的初始平均开度，μm；bmax为受压下裂缝最大开度，μm；α'为正常闭合下正常剪切下的函数斜率，无量纲；σ为正应力大小，MPa。
根据式(49)可以从一维的公式推导到二维条件，在x方向上和在y方向上裂缝开度的大小：
bx＝br+bmax＝br+bmaxexp{-(αxσx+αyσy)} (50)
式中：αx为x方向上应力系数，无量纲；αy为y方向上应力系数，无量纲；σx为x方向上正应力大小，MPa；σy为y方向上正应力大小，MPa。
同理：
by＝br+bmax＝br+bmaxexp{-(βxσx+βyσy)} (51)
式中：βx为正常闭合下正常剪切下的函数斜率，无量纲；βy为正常闭合下正常剪切下的函数斜率，无量纲；σx为x方向上正应力大小，MPa；σy为y方向上正应力大小，MPa。
接下来计算剪切应力条件下裂缝渗透率，计算公式如下：
式中：fdx为在x方向上剪切裂缝的等效频率，1/m；fdy为在y方向上剪切裂缝的等效频率，1/m；dx为在x方向上剪切裂缝的等效开度，μm；dy为在y方向上剪切裂缝的等效开度，μm。
由于只有在临界压力附近的裂缝才会剪切破坏，因此只有部分裂缝会出现剪切变形，其他大部分裂缝都不会出现剪切变形。所以发生剪切破坏的裂缝产生的渗透率会小于之前计算所得的正常闭合裂缝产生的渗透率。由裂缝剪切变形产生的开度大小可以通过下式计算：对于k<kc：
dx＝0 (54)
dy＝0 (55)
对于k≥kc：
dx＝dmax[1-exp{-γx(k-kc)}] (56)
dy＝dmax[1-exp{-γy(k-kc)}] (57)
式中：kc为裂缝变形的临界压力比，无量纲；dmax为裂缝发生剪切破坏后的最大开度，μm。
导致裂缝破坏的临界压力比和临界方向值可以通过库伦破坏准则来计算：
式中：φ为岩石内摩擦角，°；φf为岩石临界破坏角，°。
最后，x和y方向上的渗透率(kx,ky)大小即可通过叠加正常闭合和剪切变形的裂缝渗透率得到，如下：
整个裂缝网络区域的综合渗透率通过x与y方向渗透率均方根求得：
式中：kf—体积压裂区复杂裂缝网络渗透率，m2，kx为裂缝在x方向上的渗透率，10-3μ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方向上的渗透率，10-3μm2；fy为在y方向上正应力产生的裂缝的等效频率，1/m；βx为在xy下的函数斜率，无量纲；βy为在yx下的函数斜率，无量纲；fdy为在y方向上剪应力产生的裂缝等效频率，1/m。
在本发明提供的实施例中，所述步骤S3中页岩气藏压裂水平井稳态产能模型将页岩气的渗流场分为内外两个流场，内部流场为体积压裂区的达西流动，外部流场为页岩基质中的克努森流动；
其中外部流场的理想气体流动方程表达式如下：
其中
km＝kapp (64)
内部流场的流动方程如下：
式中：ν—流体流速，m/s；kf—体积压裂区复杂裂缝网络渗透率，m2；μ—气体黏度，mPa·s；
—沿程压力梯度，MPa/m；
页岩气藏压裂水平井的稳态产量方程表达式如下：
其中
式中：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。
在本发明提供的实施例中，所述步骤S3中页岩气藏压裂水平井非稳态模型把页岩气的流动分为体积压裂区和压裂区到储层边界页岩基质内外两个渗流场，页岩气的流动为径向流，同时考虑页岩气的解吸，得到页岩气非稳态下的产量；
单位质量的页岩基质在压力发生改变时解吸出的页岩气体积可由Langmuir等温吸附解吸表达式得到，如下：
式中：Q—单位质量页岩基质解吸页岩气体积，m3/kg；pi—原始地层压力，MPa；
在地面标准状况下，页岩气藏基质中的游离气量为：
式中：Gm—地下储集页岩气体积折算到地面标准状况下的体积，m3；Vm—页岩基质储集空间的体积，m3；Bgi—页岩气在原始地层压力和温度下的体积系数，无量纲；
页岩气藏基质的体积为：
式中：Vt—页岩储层视体积，m3；Sw—束缚水饱和度，无量纲；
根据朗格缪尔等温吸附理论，页岩气藏的吸附气体积为：
Ga—页岩气基质中吸附气在地面标准状况下的体积，m3；ρs—页岩基质的密度，kg/m3；
根据式(71)、(72)和(73)求得：
随着页岩气的不断被采出，地层压力下降，同时束缚水膨胀，基质孔隙收缩，裂缝中气体膨胀，令裂缝中的游离气为Gf。则储层中剩余游离气量为：
基质中：
改造区内：
式中：G'm—基质剩余游离气在地面标准状况下的体积，m3；G'f—裂缝剩余游离气在地面标准状况下的体积，m3；Gf—裂缝中的游离气在标准状况下的体积，m3；Cm—页岩基质的压缩系数，MPa-1；Bg—当前压力下的页岩气体积系数，无量纲；
根据朗格缪尔等温吸附理论，可得到页岩气藏剩余的吸附气量：
根据物质平衡原理，原始地层压力下的基质中的游离气量、原始地层压力下裂缝中的游离气量、原始地层压力下的吸附气量之和等于采出气量、当前地层压力下基质中的游离气量、当前地层压力下裂缝中的游离气量、当前地层压力下的吸附气量之和，即：
Gm+Gf+Ga＝Gp+G'm+G'f+G′a (78)
得到页岩气藏物质平衡方程表达式如下：
式中：Gm—地下储集页岩气体积折算到地面标准状况下的体积，m3；Gf—裂缝中的游离气在标准状况下的体积，m3；Gp—页岩气累计产出量，m3；
—基质孔隙度，无量纲；Sw—束缚水饱和度，无量纲；Bgi—页岩气在原始地层压力和温度下的体积系数，无量纲；ρs—页岩基质的密度，kg/m3；pi—原始地层压力，MPa；VL为朗格缪尔压力，Pa；pL为朗格缪尔体积，m3/kg；
其中：
页岩气压裂水平井生产过程中的压力传播半径方程表达式如下：
式中：R—t时刻的压力传播半径，m；K—储层渗透率，m2；t—生产时间，d；
—基质孔隙度，无量纲；μ—气体黏度，mPa·s；Ct—储层的综合压缩系数，MPa-1；
所述连续拟稳定法步骤如下：求得初始产量，假定一个时间步长，认为在该时间步长内的渗流为稳定渗流，计算出累计产量，根据物质平衡方程，算出在该时间段内压力波所波及到的区域的平均地层压力，以该平均地层压力作为下一时间步长的边界压力，计算下一时间步长下的产量，以此类推，便可得到产量与时间的关系，即非稳态下的产量。
由于页岩气藏是致密储层，为获得商业气流往往需要大规模的体积压裂。所以，在页岩储层中，压力波传播的第一阶段又可以分为两个时段，即压力波先从井底传到体积压裂区的外边界，然后再从体积压裂区的外边界传到储层边界。压力波传播的距离即压力波及半径仅与时间、储层和流体的物性有关。因此，先设定一个时间步长Δt，可以得到压力波及半径R1，再根据稳态下产量公式式得到初始产量q1。由于页岩气藏生产初期所产气体是体积压裂区的游离气，所以时间步长Δt要取小，以免压力波及超出体积压裂区的外边界，造成较大误差。此后，再取相同时间步长Δtf，得到Δt+Δtf时刻压力波及R2，根据容积法可算出在当前压力波及半径内的地质储量，假定在Δtf时间内页岩气稳定产出，得到在Δtf时间内的累计产量，根据物质平衡方程计算当前压力波及半径内的地层平均压力p，然后根据稳态下页岩气产量公式计算Δt+Δtf时刻的产量q2，重复以上工作，得到页岩气非稳态的产量。
实施例：
某井场的页岩储层单根毛细管表观渗透率的基础数据如表1所示：
表1页岩气储层相关数据表
| 参数名 | 符号 | 单位 | 数值 |
| 气体类型 | CH<sub>4</sub> | — | — |
| 分子碰撞直径 | δ | m | 0.42×10<sup>-9</sup> |
| 气体分子直径 | d<sub>m</sub> | m | 3.8×10<sup>-10</sup> |
| 气体常数 | R | J/(mol·K) | 8.314 |
| 温度 | T | K | 423 |
| 气体摩尔质量 | M | g/mol | 16 |
| 气体黏度 | μ | Pa·s | 1.84×10<sup>-5</sup> |
| 压缩因子 | z | 无量纲 | 0.85 |
| 切向动量调节系数 | α | 无量纲 | 0.8 |
| 气体分子密度 | ρ | kg/m<sup>3</sup> | 0.655 |
| 原始地层压力,储层平均压力 | P,p<sub>avg</sub> | Pa | 30×10<sup>6</sup> |
| 表面最大浓度 | C<sub>smax</sub> | mol/m<sup>3</sup> | 25040 |
| Langmuir压力 | p<sub>L</sub> | Pa | 2.46×10<sup>-6</sup> |
| 表面扩散系数 | D<sub>s</sub> | m<sup>2</sup>/s | 2.89×10<sup>-10</sup> |
现场取回岩心分析得到基质孔隙度、束缚水饱和度及毛细管管径分布结果见表2。
表2页岩岩心基质孔隙度及毛细管孔径分布
根据页岩气井压裂现场微地震数据反演以及实验测试获得的分形维数及缝网基础参数结果见表3。
表3缝网区域基础参数
| 参数 | 单位 | 数值 |
| 分形维数 | - | 1.47 |
| 最大裂缝长度r<sub>max</sub> | m | 150 |
| 最小裂缝长度r<sub>min</sub> | m | 3 |
| 裂缝区域宽度 | m | 100 |
| 裂缝区域长度 | m | 100 |
| 裂缝面密度 | 条/m<sup>2</sup> | 3 |
| b<sub>x</sub>上的应力系数(α<sub>x</sub>，α<sub>y</sub>)<sup>*</sup> | - | 0.034，0.146 |
| b<sub>y</sub>上的应力系数(β<sub>x</sub>，β<sub>y</sub>)<sup>*</sup> | - | 0.101，0.071 |
| 剩余开度(b<sub>r</sub>) | μm | 50 |
| 最大开度(b<sub>max</sub>) | μm | 300 |
| 最大剪切开度(d<sub>max</sub>) | μm | 500 |
| 裂缝初始开度(b<sub>i</sub>) | μm | 300 |
| 摩擦角<sup>*</sup> | ° | 24.9 |
| 临界压力比 | 无因次 | 2.45 |
| 有效水平主应力(σ<sub>x</sub>) | MPa | 37.8～65 |
| 有效水平主应力(σ<sub>y</sub>) | MPa | 32.5～67 |
计算页岩气压裂水平井非稳态产量需要的基础数据见表4。
表4地层基础参数
根据以上参数，通过式(36)计算得到整个岩心的表观渗透率，如图1所示，图1表示页岩基质表观渗透率计算结果。
根据以上表3数据，得到页岩压裂水平井微地震反演裂缝分布图像，如图2所示；根据微地震反演出的裂缝分布图像进行灰度处理，识别出该区域内的裂缝分布，如图3所示；通过盒子法计算出裂缝区域的分形维数图，如图4所示。
根据以上数据，通过式(62)计算得到复杂裂缝网络的渗透率，分析得到x方向与y方向水平主应力相同时渗透率的变化情况，如图5所示；x方向与y方向水平主应力不同时渗透率的变化情况，如图6所示。
根据以上数据，通过连续拟稳定法得到页岩气藏压裂水平井复杂裂缝网络的非稳态产量，且本发明预测结果与Eclipse模型对比结果如图7所示，从图7可以看出本发明计算结果与该井的实际生产情况吻合，且误差较小；通过对计算结果分析，得到不同分形维数下页岩气压裂水平井累产的变化情况，如图8所示；得到不同压裂半径下，页岩气压裂水平井累产的变化情况，如图9所示；得到不同裂缝开度下，页岩气压裂水平井累产的变化情况，如图10所示。
以上所述，并非对本发明作任何形式上的限制，虽然本发明已通过实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容作出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。
Claims (4)
1.一种计算页岩气藏复杂裂缝网络非稳态产量的方法，其特征在于，包括以下步骤：
S1、收集页岩气藏地层、页岩气的基本参数，建立页岩基质多孔介质毛细管渗透率模型；
S2、收集页岩气藏微地震裂缝成像数据，根据微地震裂缝成像数据，识别出裂缝条数、位置及形态特征；运用盒子法获取表征岩石复杂裂缝网络的分形维数；根据求取的分形维数，利用蒙特卡洛随机建模理论生成二维离散裂缝网络模型；通过立方定律计算复杂裂缝网络渗透率；
S3、假设储层等厚、圆形边界、半径为re，体积压裂区域改造充分、圆形边界、半径为rf、与储层厚度相同；考虑页岩气藏压裂水平井外区为基质区域渗流区域、内区为复杂缝网区域流动区域，建立页岩气藏压裂水平井稳态产能模型；建立页岩气压力传播半径方程、页岩气藏物质平衡方程，运用连续拟稳定法建立页岩气藏压裂水平井非稳态模型，求解页岩气藏压裂水平井非稳态产量；
所述步骤S1中页岩基质多孔介质毛细管渗透率模型是根据毛细管大小的分布不同分别计算不同尺度下毛细管的渗透率，然后对不同尺度下的毛细管渗透率进行叠加，从而得到整个岩心的渗透率大小，其表达式如下：
式中：kapp-整个岩心的表观渗透率，m2；
-基质孔隙度，无量纲；kapp，i-不同尺度毛细管的表观渗透率，m2；λi-不同尺度毛细管的分布频率，无量纲；i，n-计数单位，无量纲；τ-迂曲度，无量纲，其中迂曲度计算公式如下：
式中：m为迂曲度拟合参数，取值0.77；
页岩储层中普遍具有一定的束缚水饱和度，不同页岩尺寸页岩毛细管半径对应的有效流动半径的计算公式如下：
式中：rie-不同尺寸毛细管对应的有效流动半径，mm；ri-不同尺寸毛细管对应的孔隙半径，mm；Sw-页岩孔隙束缚水饱和度，无量纲；i-计数符号，无量纲；
不同孔径尺度下的毛细管表观渗透率计算公式如下：
式中：kapp，i-不同尺寸毛细管的表观渗透率，m2；i-计数符号，无量纲；dm-气体分子直径，m；rie-不同尺寸毛细管对应的有效流动半径，m；p-储层压力，Pa；pL-Langmuir压力，Pa；Fi-不同尺寸毛细管有效流动半径对应的滑脱系数，无量纲；εi-不同尺寸毛细管有效流动半径对应的贡献系数，无量纲；Dki-不同尺寸毛细管有效流动半径对应的克努森扩散系数，m2/s；ρ—气体密度，kg/m3；μ-气体黏度，Pa·s；M-气体摩尔质量，g/mol；Ds-表面扩散系数，m2/s；Csmax-吸附气最大吸附度，mol/m3；p-储层压力，Pa；
所述步骤S2中复杂裂缝网络渗透率通过微地震数据地震成像获取分形维数，通过蒙特卡洛随机建模理论建立二维离散裂缝网络模型，并进一步通过立方定律来计算压裂施工后区域内复杂缝网的渗透率，其表达式如下：
式中：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。
2.根据权利要求1所述的一种计算页岩气藏复杂裂缝网络非稳态产量的方法，其特征在于，所述页岩气藏地层、页岩气的基本参数包括储层边界压力、储层温度、供给半径、基质孔隙度、裂缝孔隙度、体积压裂区半径、井眼半径、储层厚度、页岩气粘度、气体类型、分子碰撞直径、气体分子直径、气体常数、气体摩尔质量、切向动量调节系数、气体分子密度、平均压力、表面最大浓度、朗格缪尔压力和表面扩散系数。
3.根据权利要求1所述的一种计算页岩气藏复杂裂缝网络非稳态产量的方法，其特征在于，所述步骤S1中页岩基质多孔介质毛细管渗透率模型，根据达西定律求得页岩气藏整个岩心的表观渗透率的步骤是利用页岩气储层参数计算得到克努森系数，根据克努森系数判断毛细管中的流态，利用毛细管中的流态和运移机制建立气体质量运移方程，接下来建立统一的质量运移方程，利用统一的质量运移方程得到毛细管的表观渗透率，最后进行叠加求得整个岩心的表观渗透率；
其中克努森系数的计算公式如下：
式中：Kn-克努森系数，无因次；kB-玻尔兹曼常数，1.3805×10-23J/K；p-储层压力，Pa；T—页岩储层温度，K；π-常数，3.14；δ-气体分子碰撞直径，m；d-毛细管直径，m。
4.根据权利要求1所述的一种计算页岩气藏复杂裂缝网络非稳态产量的方法，其特征在于，所述步骤S2中复杂裂缝网络渗透率包括正常张开型裂缝渗透率、剪切型裂缝渗透率，计算渗透率过程中，正应力和剪应力分别进行计算，最后将正应力和剪应力诱导产生裂缝所具有渗透率进行叠加再复合得到，表达式如下：
kx＝knx+kdx
ky＝kny+kdy
式中：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。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810311636.6A CN108518212B (zh) | 2018-04-09 | 2018-04-09 | 一种计算页岩气藏复杂裂缝网络非稳态产量的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810311636.6A CN108518212B (zh) | 2018-04-09 | 2018-04-09 | 一种计算页岩气藏复杂裂缝网络非稳态产量的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108518212A CN108518212A (zh) | 2018-09-11 |
| CN108518212B true CN108518212B (zh) | 2020-10-16 |
Family
ID=63430771
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810311636.6A Active CN108518212B (zh) | 2018-04-09 | 2018-04-09 | 一种计算页岩气藏复杂裂缝网络非稳态产量的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108518212B (zh) |
Families Citing this family (35)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109268004B (zh) * | 2018-09-21 | 2020-06-26 | 北京科技大学 | 一种页岩气藏介质结构耦合及缝网形态判识方法 |
| CN109342287B (zh) * | 2018-10-25 | 2021-03-09 | 中国石油天然气股份有限公司 | 一种气水稳态渗流的判定方法 |
| CN111188610A (zh) * | 2018-10-29 | 2020-05-22 | 中国石油化工股份有限公司 | 一种致密气藏压裂气井产能确定方法及装置 |
| CN109658016B (zh) * | 2019-02-12 | 2022-11-22 | 中国石油大学(北京) | 致密气藏补给边界的识别方法及系统 |
| CN111626530B (zh) * | 2019-02-27 | 2023-08-25 | 中国石油天然气股份有限公司 | 确定压裂页岩气藏井控储量的方法及装置 |
| CN110018293B (zh) * | 2019-05-05 | 2021-09-21 | 西安石油大学 | 一种考虑水敏性的页岩黏土多相耦合甲烷含气量计算方法 |
| CN110206522B (zh) * | 2019-06-10 | 2019-12-20 | 西南石油大学 | 一种页岩气藏压裂水平井压裂液返排模拟方法 |
| CN110295885B (zh) * | 2019-06-27 | 2022-04-01 | 西南石油大学 | 一种考虑页岩多尺度效应的流态评价方法 |
| CN110671100B (zh) * | 2019-10-10 | 2022-08-30 | 东北石油大学 | 一种制造模拟岩石非均质性装置中棋盘状仿真体的方法 |
| CN111177946B (zh) * | 2020-01-10 | 2022-03-11 | 西南石油大学 | 高压致密气藏射孔总表皮和压实带表皮系数模拟计算方法 |
| CN111353227B (zh) * | 2020-02-28 | 2023-03-14 | 西安石油大学 | 基于跨尺度多流动空间气体传输机制的co2强化页岩气藏开发动态仿真方法 |
| CN111340298B (zh) * | 2020-02-28 | 2021-11-02 | 西安石油大学 | 页岩储层诱导非均质性动态演化预测方法 |
| CN111444612B (zh) | 2020-03-26 | 2021-04-16 | 北京科技大学 | 一种致密油藏水平井多级压裂流场形态模拟方法 |
| CN111650088B (zh) * | 2020-06-10 | 2021-07-27 | 河海大学 | 一种流态混凝土拌合物流变性能实时检测方法 |
| CN111832164B (zh) * | 2020-07-03 | 2023-05-12 | 中国石油大学(北京) | 页岩气产量预测方法、装置及设备 |
| CN111852431B (zh) * | 2020-08-19 | 2022-04-12 | 西南石油大学 | 一种缝网压裂中缝网结构参数的优化方法和装置 |
| CN112084689B (zh) * | 2020-08-25 | 2023-12-15 | 中海油田服务股份有限公司 | 天然气储层的非稳态渗流模拟方法及系统 |
| CN111983192B (zh) * | 2020-08-25 | 2022-03-11 | 西南石油大学 | 向斜背景页岩气大量逸散深度定量确定方法 |
| CN112001134B (zh) * | 2020-09-04 | 2021-09-07 | 北京科技大学 | 一种非常规气藏流场结构识别方法及装置 |
| CN112102487B (zh) * | 2020-09-17 | 2022-02-15 | 西南石油大学 | 一种基于多重混合分形的非常规储层三维渗透率确定方法 |
| CN112179815B (zh) * | 2020-09-21 | 2022-04-05 | 西南石油大学 | 一种基于孔隙网络模型的单相非稳态渗流模型建立方法 |
| CN112081583B (zh) * | 2020-09-25 | 2021-08-17 | 西南石油大学 | 非常规储层缝网压裂多尺度支撑产量计算方法及装置 |
| CN112434426B (zh) * | 2020-11-26 | 2022-09-27 | 北京科技大学 | 页岩气多级压裂水平井台阶梯度压降开发方法及装置 |
| CN112360422B (zh) * | 2020-12-08 | 2022-04-08 | 西南石油大学 | 一种页岩气藏压裂水平井产量预测方法及系统 |
| CN113034003B (zh) * | 2021-03-26 | 2022-03-15 | 西南石油大学 | 一种页岩气井产能快速评价方法 |
| CN115341894A (zh) * | 2021-04-27 | 2022-11-15 | 中国石油天然气股份有限公司 | 油气井产能评价方法和装置 |
| CN114183113B (zh) * | 2021-05-19 | 2024-06-14 | 中国石油天然气股份有限公司 | 一种简化表征页岩气井压裂裂缝形态的方法 |
| CN115680592B (zh) * | 2021-07-28 | 2024-06-25 | 中国石油天然气股份有限公司 | 压裂有效性参数的确定方法、装置、设备及存储介质 |
| CN113792932B (zh) * | 2021-09-18 | 2023-11-07 | 西南石油大学 | 一种利用微震-损伤-渗流关系的页岩气产量预测方法 |
| CN114048695B (zh) * | 2021-11-10 | 2022-09-09 | 西南石油大学 | 一种基于返排数据的页岩气有效缝网体积反演方法 |
| CN115270663B (zh) * | 2022-09-27 | 2022-12-09 | 中国石油大学(华东) | 一种描述页岩基质与天然裂缝间非稳态窜流的方法及系统 |
| CN115600425B (zh) * | 2022-11-07 | 2023-05-02 | 中国科学院武汉岩土力学研究所 | 一种基于蛭石热膨胀致裂的页岩可压性评价装置及方法 |
| CN116629165B (zh) * | 2023-07-24 | 2023-09-22 | 中国石油大学(华东) | 储层压裂改造区和未改造区参数反演方法、系统及设备 |
| CN117633409B (zh) * | 2024-01-25 | 2024-04-09 | 中国科学院地质与地球物理研究所 | 计算页岩油气储层压裂缝网渗流参数的方法、系统及设备 |
| CN118095015B (zh) * | 2024-04-23 | 2024-07-23 | 中国科学院武汉岩土力学研究所 | 一种裂隙岩体渗流性能预测方法、装置以及处理设备 |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2016028564A1 (en) * | 2014-08-22 | 2016-02-25 | Schlumberger Canada Limited | Methods for monitoring fluid flow and transport in shale gas reservoirs |
| CN105046006B (zh) * | 2015-07-29 | 2018-05-04 | 中国石油天然气股份有限公司 | 一种页岩气藏水平井多段压裂产能预测方法及装置 |
| CN107066674B (zh) * | 2017-01-18 | 2021-08-10 | 西南石油大学 | 计算页岩气藏体积压裂水平井非稳态产量的方法 |
-
2018
- 2018-04-09 CN CN201810311636.6A patent/CN108518212B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN108518212A (zh) | 2018-09-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108518212B (zh) | 一种计算页岩气藏复杂裂缝网络非稳态产量的方法 | |
| CN109441422B (zh) | 一种页岩气井间距优化开采方法 | |
| CN107301306B (zh) | 用于致密砂岩气藏压裂水平井的动态无阻流量预测方法 | |
| Wu et al. | A model for gas transport in micro fractures of shale and tight gas reservoirs | |
| Javadpour et al. | Nanoscale gas flow in shale gas sediments | |
| CN105089663B (zh) | 一种大斜度/水平井地层电阻率各向异性校正方法 | |
| Li et al. | Composition-transient analysis in shale-gas reservoirs with consideration of multicomponent adsorption | |
| CN107622328A (zh) | 一种页岩气藏分段多簇压裂水平井产量预测方法 | |
| CN106932323B (zh) | 一种页岩气藏气体有效孔隙度反演方法 | |
| CN111927417B (zh) | 一种页岩气分段压裂水平井组储量动用状况评价方法 | |
| CN110984973A (zh) | 缝洞型碳酸盐岩气藏单井控制储量的确定方法 | |
| CN103258091A (zh) | 非常规油气藏水平井段三维岩体力学模型建立的方法及装置 | |
| CN106054279B (zh) | 一种煤岩脆性指数的确定方法 | |
| CN107130959A (zh) | 一种煤层气产量预测方法 | |
| Ruth et al. | Combining electrical measurements and mercury porosimetry to predict permeability | |
| CN109184676A (zh) | 一种页岩气藏有效改造体积评价方法 | |
| CN113947041A (zh) | 一种页岩储层压裂缝网络扩展流动一体化模拟方法 | |
| Li et al. | Effect of adsorption and permeability correction on transient pressures in organic rich gas reservoirs: Vertical and hydraulically fractured horizontal wells | |
| Li et al. | A novel multiphysics multiscale multiporosity shale gas transport model for geomechanics/flow coupling in steady and transient states | |
| Yang et al. | Impact of unrecovered shale gas reserve on methane emissions from abandoned shale gas wells | |
| CN116502553B (zh) | 非常规油气藏裂缝堵塞表皮系数及裂缝参数的反演方法 | |
| CN117252054A (zh) | 基于双重介质模型的页岩气压裂水平井产量预测方法 | |
| Yildiz | Productivity of selectively perforated vertical wells | |
| CN113348458B (zh) | 评估非均质地层中的烃的方法和系统 | |
| CN116291416A (zh) | 一种考虑井间干扰的页岩气加密井产能预测方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |