CN112012731A - 一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法 - Google Patents
一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法 Download PDFInfo
- Publication number
- CN112012731A CN112012731A CN202010973613.9A CN202010973613A CN112012731A CN 112012731 A CN112012731 A CN 112012731A CN 202010973613 A CN202010973613 A CN 202010973613A CN 112012731 A CN112012731 A CN 112012731A
- Authority
- CN
- China
- Prior art keywords
- gas
- phase
- water
- equation
- pressure
- 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.)
- Granted
Links
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 210
- 238000000034 method Methods 0.000 title claims abstract description 88
- 230000005514 two-phase flow Effects 0.000 title claims abstract description 78
- 238000010276 construction Methods 0.000 title claims description 13
- 230000035699 permeability Effects 0.000 claims abstract description 55
- 239000005416 organic matter Substances 0.000 claims abstract description 51
- 239000012071 phase Substances 0.000 claims description 203
- 239000000126 substance Substances 0.000 claims description 56
- 238000004519 manufacturing process Methods 0.000 claims description 51
- 239000011435 rock Substances 0.000 claims description 39
- 239000011159 matrix material Substances 0.000 claims description 33
- 238000012937 correction Methods 0.000 claims description 10
- 239000008346 aqueous phase Substances 0.000 claims description 8
- 239000012530 fluid Substances 0.000 claims description 4
- 210000002615 epidermis Anatomy 0.000 claims description 3
- 239000011148 porous material Substances 0.000 abstract description 18
- 230000009471 action Effects 0.000 abstract description 11
- 230000000694 effects Effects 0.000 abstract description 10
- 238000009792 diffusion process Methods 0.000 abstract description 9
- 238000001179 sorption measurement Methods 0.000 abstract description 8
- 230000007246 mechanism Effects 0.000 abstract description 6
- 230000008859 change Effects 0.000 abstract description 5
- 238000003795 desorption Methods 0.000 abstract description 5
- 230000000704 physical effect Effects 0.000 abstract description 4
- 238000010521 absorption reaction Methods 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 180
- 238000004088 simulation Methods 0.000 description 10
- 238000012545 processing Methods 0.000 description 9
- 239000002734 clay mineral Substances 0.000 description 6
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000013178 mathematical model Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 238000000418 atomic force spectrum Methods 0.000 description 3
- 230000008878 coupling Effects 0.000 description 3
- 238000010168 coupling process Methods 0.000 description 3
- 238000005859 coupling reaction Methods 0.000 description 3
- 238000009795 derivation Methods 0.000 description 3
- 230000009977 dual effect Effects 0.000 description 3
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 230000005465 channeling Effects 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 230000006835 compression Effects 0.000 description 2
- 238000007906 compression Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 239000003345 natural gas Substances 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 150000001875 compounds Chemical class 0.000 description 1
- 230000008094 contradictory effect Effects 0.000 description 1
- 239000013078 crystal Substances 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000006073 displacement reaction Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 229910010272 inorganic material Inorganic materials 0.000 description 1
- 239000011147 inorganic material Substances 0.000 description 1
- 238000013508 migration Methods 0.000 description 1
- 230000000644 propagated effect Effects 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000002336 sorption--desorption measurement Methods 0.000 description 1
- 238000011144 upstream manufacturing Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH DRILLING; MINING
- E21B—EARTH DRILLING, e.g. DEEP DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B49/00—Testing the nature of borehole walls; Formation testing; Methods or apparatus for obtaining samples of soil or well fluids, specially adapted to earth drilling or wells
-
- E—FIXED CONSTRUCTIONS
- E21—EARTH DRILLING; MINING
- E21B—EARTH DRILLING, e.g. DEEP DRILLING; OBTAINING OIL, GAS, WATER, SOLUBLE OR MELTABLE MATERIALS OR A SLURRY OF MINERALS FROM WELLS
- E21B47/00—Survey of boreholes or wells
- E21B47/06—Measuring temperature or pressure
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- 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—Systems or methods specially adapted for specific business sectors, e.g. utilities or tourism
- G06Q50/02—Agriculture; Fishing; Mining
Abstract
本发明公开了一种考虑气水两相流动的页岩气藏三孔三渗模型及动态预测方法，考虑页岩气藏由物理性质独立的有机质系统、无机质系统和裂缝系统组成，同时考虑了无机质纳米级孔隙吸附水层的影响对渗透率进行了修正，综合考虑页岩气流动为压力差作用下的黏性流动、有机质和黏土矿物的吸附解吸、滑脱效应和克努森扩散多重机制的作用，推导得到了页岩气藏气水两相流动三孔三渗模型；采用IMPES方法，利用Matlab数学软件对页岩气藏多级压裂水平井两相三维三孔三渗模型进行编程求解，精准预测的产量压力动态变化。旨在解决现有技术中存在的页岩气藏气水两相双重介质模型进行产量压力动态预测不够精准的技术问题。
Description
技术领域
本发明涉及页岩气开发技术领域，尤其涉及一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法。
背景技术
页岩气作为我国的重要非常规能源，为了实现对页岩气藏经济开采常采用压裂技术，压裂后会出现气水两相流动现象。目前国内外学者对页岩气藏气水两相渗流规律的研究多是基于基质-裂缝双重介质模型，实现了页岩气藏气水两相流动的产气、产水预测，但欠缺对基质中亲油性有机质和亲水性无机质物理性质、吸附解吸能力差异性的考虑。在页岩储层中有机质通常以集合体形式与黏土矿物共生，有机质内发育大量吸附特性的纳米级孔隙，无机质内除了发育水湿性的纳米～微米级孔隙，同时无机质中黏土矿物因其特殊的晶体结构同样也具备一定的甲烷吸附能力。现有研究认为由于基质渗透率很低，整个基质中往往假设为单相气体流动，未考虑亲油性有机质和亲水性无机质的差异。事实上基质中的黏土矿物和脆性矿物除了发育大量纳米级孔隙外还存在一些微米级微裂缝，当地层为气水两相流动的时候，水会进入到无机质中甚至吸附在黏土矿物表面，而水通常不会进入到亲油性较强的有机质里。因此，有必要将基质中的有机质和无机质分开考虑进一步认识基质中气水两相流动规律，对更加精细描述页岩气藏气水两相流动问题具有重要意义。
页岩气藏渗流模型的研究正在逐步完善，在传统双重孔隙介质模型的研究基础上又发展了多重介质多重运移机制的渗流模型，但关于页岩气藏气水两相渗流规律的研究还处于探究阶段，特别是基质中气水两相流动规律认识还不够清楚。目前的页岩气藏气水两相渗流模型还存在以下问题：
(1)目前的模型多是基于基质-裂缝双重介质模型，未能体现亲油性有机质和亲水性无机质的差异；
(2)未能对有机质和无机质的气水两相流动规律分别考虑，未建立有机质 -无机质-裂缝共同形成三个系统的气水两相流动三孔三渗模型；
(3)未考虑吸附水层对无机质孔喉半径和孔隙度的影响，在渗透率修正中也未能体现吸附水层的影响。
因此，如何构建更精准的页岩气藏气水两相双重介质模型进行产量压力动态预测，是一个亟需解决的技术问题。
上述内容仅用于辅助理解本发明的技术方案，并不代表承认上述内容是现有技术。
发明内容
本发明的主要目的在于提供一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，旨在解决现有技术中存在的页岩气藏气水两相双重介质模型进行产量压力动态预测不够精准的技术问题。
为实现上述目的，本发明提出一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，所述三孔三渗模型构建及压力动态预测方法包括如下步骤：
构建气水两相流动三孔三渗模型：
S101：分别构建裂缝系统、无机质系统的气水两相流动方程、毛管力辅助方程和饱和度辅助方程以及有机质系统的气相流动方程；
S102：基于获取的裂缝系统、无机质系统和有机质系统的初始条件及内外边界条件，将饱和度辅助方程、毛管力辅助方程代入气水两相流动方程，分别得到裂缝系统、无机质系统和有机质系统下的气相压力方程；
S103：利用有限差分法得到气相压力的有限差分方程，通过IMPES方法对有限差分方程做线性化处理得到气相压力的线性代数方程组，将该线性代数方程组作为水平井气水两相流动三孔三渗模型；
压力动态预测：
S201：对存在于裂缝系统网络、无机质系统系统和有机质系统网络中的页岩气井进行点源或点汇数学处理，将网格内井的产量用拟稳态流公式表示得到井的生产指数；
S202：采用定流压方式生产，利用MATLAB对多级压裂水平井气水两相流动双孔双渗数值模型进行编程求解，获取该井以定井底流压生产时的产气量与产水量。
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，所述构建裂缝系统的气水两相流动方程、毛管辅助方程和饱和度辅助方程包括：
A1：将单位体积页岩无机质向裂缝窜流水量表达式带入裂缝网络中气相渗流微分方程，得到裂缝系统下的气相方程：
式中：▽为Hamilton算子；kf为裂缝系统渗透率，m2；kfrg为裂缝网络气相相对渗透率，无量纲；Bg为气体的体积系数，无量纲；μg为气体粘度，mPa·s；pfg为裂缝网络气相压力，MPa；ρfg为气体密度，kg/m3；D为气体所处的垂直深度，m；qsg为地面条件下单位体积页岩中采出的气量，sm3/s；α为形状因子，1/m2；ξc为页岩岩石无机质渗透率修正系数，无量纲；kcrg为无机质网络气相相对渗透率，无量纲；pcg为无机质系统气相压力，MPa；pfg为裂缝系统气相压力，MPa；φf为裂缝系统孔隙度，无量纲；sfg为裂缝网络气相饱和度，无量纲。
A2：将单位体积页岩无机质向裂缝窜流水量的表达式带入裂缝网络中水相渗流方程，得到裂缝系统下的水相方程：
式中：kfrw为裂缝网络水相相对渗透率，无量纲；Bw为水相的体积系数，无量纲；μw为水相粘度，mPa·s；pfw为裂缝系统水相压力，MPa；ρfw为水相密度， kg/m3；qsw为地面条件下单位体积页岩中采出的水量，sm3/s；kc为无机质系统渗透率，m2；kcrw为基质网络水相相对渗透率，无量纲；Bw为水相的体积系数，无量纲；μw为水相粘度，mPa·s；pcw为无机质网络水相压力，MPa；pfw为裂缝网络水相压力，MPa；sfw为裂缝网络水相饱和度，无量纲。
A3：获取裂缝网络中的气相饱和度和水相饱和度，得到饱和度辅助方程：
sfg+sfw＝1；
式中：sfg为裂缝网络气相饱和度，无量纲；sfw为裂缝网络水相饱和度，无量纲。
A4：获取裂缝网络中的毛管力，得到毛管辅助方程：
pfc(sfw)＝pfg-pfw；
式中：pfc为裂缝网络中的毛管力，MPa；pfg为裂缝网络气相压力，MPa；pfw为裂缝网络水相压力，MPa。
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，所述构建基质系统的气水两相流动方程、毛管辅助方程和饱和度辅助方程包括：
B1：将单位体积页岩无机质向裂缝窜流气量表达式带入无机质网络中气相渗流微分方程，得到基质系统下的气相方程：
式中：kc为无机质固有渗透率，m2；ζc为页岩岩石无机质渗透率修正系数，无量纲；kcrg为无机质网络气相相对渗透率，无量纲；pcg为页岩岩石无机质中气相压力，MPa；Smg为基质网络中气相饱和度，无量纲；φm为基质孔隙度；VL为基质的兰氏体积，sm3/m3；pL为有基质的兰氏压力，MPa。
B2：将单位体积页岩无机质向裂缝窜流水量的表达式带入无机质网络中水相渗流方程，得到基质系统下的水相方程：
式中：kc为页岩岩石无机质固有渗透率，m2；kcrw为页岩岩石无机质中水相相对渗透率，无量纲；pcw为页岩岩石无机质中水相压力，MPa；qscfw为地面条件下单位体积页岩中无机质向裂缝的供给水量，sm3/s；scw为页岩岩石无机质中水相饱和度，无量纲。
B3：获取无机质网络中的气相饱和度和水相饱和度，得到饱和度辅助方程：
scg+scw＝1；
B4：获取无机质网络中的毛管力，得到毛管辅助方程：
pcc(scw)＝pcg-pcw；
式中：pcc为页岩岩石无机质中的毛管力，MPa。
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，构建有机质系统的气相流动方程通过质量守恒定律推导得到有机质系统中气相渗流微分方程：
式中：ko为页岩岩石有机质固有渗透率，m2；ζo为页岩岩石有机质渗透率修正系数，无量纲；pog为页岩岩石有机质中气相压力，MPa；sog为页岩岩石有机质中气相饱和度，无量纲。
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，所述获取的裂缝系统、无机质系统和有机质系统的初始条件及内外边界条件包括：
C1：裂缝系统和基质系统下的初始气相压力、初始含水饱：
式中：
为裂缝系统中初始气相压力，MPa；
为无机质系统初始气相压力，MPa；
为有机质系统初始气相压力，MPa；
为裂缝系统初始含水饱和度，无量纲；
为无机质系统初始含水饱和度，无量纲；
C2：页岩气藏三孔三渗模型的外边界条件：
C3：页岩气藏三孔三渗模型的内边界条件：
其中：
式中：pgi,j,k为网格坐标为i,j,k的气体压力，MPa；δ(i,j,k)为Dirac函数，(存在源汇项时值为1反之为0)；pwf为井底流压，MPa；S为表皮系数，无量纲；qvg为产气量，sm3/s；qvw为产水量，sm3/s。
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，所述得到裂缝系统下的气相压力方程，通过将毛管力辅助方程和饱和度辅助方程代入气水两相渗流方程，得到只含有气相压力的方程：
D1：将毛管力辅助方程代入三孔三渗模型中裂缝系统的两相渗流方程之中：
式中：
D2：将饱和度辅助方程代入上述两相渗流方程式子，并分别乘以气体体积系数Bg和水体积系数Bw，相加得到裂缝系统气相的压力方程：
其中：
CGfg＝-▽·[λfg▽(ρfggD)]
CGfw＝-▽·[λfw▽(pfc+ρfwgD)]
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，利用有限差分法得到气相压力的有限差分方程，通过IMPES方法对有限差分方程做线性化处理得到气相压力的线性代数方程组：
对三孔三渗模型的裂缝系统压力隐式求解，得到压力线性方程组：
其中：
式中：N为求解网格中北向系数；S为求解网格中南向系数；W为求解网格中西向系数；E为求解网格中东向系数；C为求解网格中中心系数；A为求解网格中上表面系数；B为求解网格中下表面系数；Qi,j,k为源汇项；下标i,j,k表示网格坐标；上标n表示求解时间步。
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，所述得到无机质系统下的气相压力方程，通过将毛管力辅助方程和饱和度辅助方程代入气水两相渗流方程，得到只含有气相压力的方程：
E1：将毛管力辅助方程代入三孔三渗模型中无机质系统的两相渗流方程之中：
其中：
E2：将上述两相渗流方程式子代入饱和度辅助方程，并分别乘以气体体积系数Bg和水体积系数Bw，相加得到无机质系统气相的压力方程：
其中：
CGcg＝-▽·[λcg▽(ρcggD)]
CGcw＝-▽·[λcw▽(pcc+ρcwgD)]
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，利用有限差分法得到气相压力的有限差分方程，通过IMPES方法对有限差分方程做线性化处理得到气相压力的线性代数方程组：
对三孔三渗模型的无机质系统压力隐式求解，得到压力线性方程组：
其中：
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，所述利用有限差分法得到气相压力的有限差分方程，通过IMPES方法对有限差分方程做线性化处理得到气相压力的线性代数方程组：
对三孔三渗模型的有机质系统压力隐式求解，得到压力线性方程组：
其中：
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，所述对存在于裂缝系统网络、无机质系统网络和有机质系统网络中的页岩气井进行点源或点汇数学处理，将网格内井的产量用拟稳态流公式表示得到井的生产指数，其中：
拟稳态流公式为：
井的生产指数为：
式中：Δn为在n方向的网格步长，m；pwf为井底流压，MPa；PID为井的生产指数；krl为l相流体的相对渗透率；re为井点处网格块的等效半径，m；rw为井筒半径，m；S为表皮因子，无量纲；pf为井底流压，Pa；Δn为在n方向的网格步长，m。
优选的，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，所述采用定流压方式生产，利用MATLAB对多级压裂水平井气水两相流动双孔双渗数值模型进行编程求解，获取该井以定井底流压生产时的产气量与产水量，其中：
在定流压生产时，第k完井段的产气量可表示为：
第k完井段的产水量为：
本发明中，本发明提供一种考虑气水两相流动的页岩气藏三孔三渗模型及动态预测方法，考虑页岩气藏由物理性质独立的有机质系统、无机质系统和裂缝系统组成，同时考虑了无机质纳米级孔隙吸附水层的影响对渗透率进行了修正，综合考虑页岩气流动为压力差作用下的黏性流动、有机质和黏土矿物的吸附解吸、滑脱效应和克努森扩散多重机制的作用，推导得到了页岩气藏气水两相流动三孔三渗模型；采用IMPES方法，利用Matlab数学软件对页岩气藏多级压裂水平井两相三维三孔三渗模型进行编程求解，精准预测的产量压力动态变化。旨在解决现有技术中存在的页岩气藏气水两相双重介质模型进行产量压力动态预测不够精准的技术问题。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图示出的结构获得其他的附图。
图1为本发明提出的页岩气藏气-水两相流动双孔双渗物理模型示意图；
图2为吸附水层的无机质纳米孔隙模型示意图；
图3为本发明提出的一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法的利用计算机求解数值模型的步骤示意图；
图4为本发明提出的地质网络模型示意图；
图5为基质系统中毛管力曲线示意图；
图6为裂缝系统中毛管力曲线示意图；
图7为基质和裂缝系统相对渗透率曲线示意图；
图8为三孔三渗模型模拟生产第10天示意图；
图9为三孔三渗模型模拟生产第50天示意图；
图10为三孔三渗模型模拟生产第300天示意图；
图11为三孔三渗模型模拟生产第600天示意图。
本发明目的的实现、功能特点及优点将结合实施例，参照附图做进一步说明。
具体实施方式
应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
另外，各个实施例之间的技术方案可以相互结合，但是必须是以本领域普通技术人员能够实现为基础，当技术方案的结合出现相互矛盾或无法实现时应当人认为这种技术方案的结合不存在，也不在本发明要求的保护范围之内。
在本实施例中，一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，包括如下构建及预测步骤：
S1、建立页岩气藏气-水两相流动三孔三渗物理模型；
S2、建立裂缝系统的气水两相流动方程；
S3、建立无机质系统的气水两相流动方程；
S4、建立有机质系统的气相流动方程；
S5、确定气水两相流动三孔三渗模型的初始条件及内外边界条件；
S6、利用有限差分法对数学模型进行离散化求得有限差分方程组，通过 IMPES方法对其做线性化处理得到线性代数方程组，求解得到相应的数值模型；
S7、对页岩气井进行点源或点汇数学处理，采用定流压方式生产，利用Matlab对多级压裂水平井气水两相流动三孔三渗数值模型进行编程求解，分析压力动态变化。
步骤S1：
在本实施例中，建立页岩气藏气-水两相流动双孔双渗物理模型具体步骤为：
如图1所示，将有机质、无机质与裂缝视为三重介质分开考虑，同时无机质纳米孔隙表面有吸附水层的影响，页岩气藏气-水两相流动三孔三渗物理模型如图1所示。
页岩气藏气水两相流动三孔三渗模型具体基本假设条件如下：
(1)流动过程中忽略温度的变化；
(2)气相和水相互不相溶；
(3)裂缝中同时存在气-水两相流动且遵循达西定律；
(4)考虑压裂液在毛管力的作用下可渗吸进入无机质中，气体在无机质中的流动为压力差作用下的黏性流动、吸附解吸、滑脱以及克努森扩散多重机制的作用；
(5)有机质中不考虑水的作用，只存在单相气体渗流，气体在有机质纳米级孔隙中的流动为压力差作用下的黏性流动、吸附解吸、滑脱以及克努森扩散多重机制的作用；
(6)毛管力仅存在于无机质和裂缝系统中，有机质中忽略毛管力；
(7)考虑重力作用。
步骤S2：
在本实施例中，建立裂缝系统的气水两相流动方程具体步骤为：
(1)气相方程
根据渗流力学基本原理，可以得到裂缝网络中气相渗流方程：
式中，▽为Hamilton算子；kf为裂缝系统渗透率，m2；kfrg为裂缝网络气相相对渗透率，无量纲；Bg为气体的体积系数，无量纲；μg为气体粘度，mPa·s；pfg为裂缝网络气相压力，MPa；ρfg为气体密度，kg/m3；D为气体所处的垂直深度，m；qsg为地面条件下单位体积页岩中采出的气量，sm3/s；qscfg为地面条件下单位体积页岩中无机质向裂缝的供给气量，sm3/s；φf为裂缝系统孔隙度，无量纲；sfg为裂缝网络气相饱和度，无量纲。
将单位体积页岩中无机质向裂缝窜流气量的表达式代入上式可得：
(2)水相方程
同样，推导得到裂缝网络中水相渗流方程：
式中，kfrw为裂缝网络水相相对渗透率，无量纲；pfw为裂缝网络中水相压力，MPa；qscfw为地面条件下单位体积页岩中无机质向裂缝的供给水量，sm3/s；qsw为地面条件下单位体积页岩中采出的水量，sm3/s；sfw为裂缝网络水相饱和度，无量纲；
将单位体积页岩中无机质向裂缝窜流水量的表达式代入上式可得：
(3)辅助方程
裂缝网络中气相饱和度和水相饱和度满足以下方程：
式中：sfg为裂缝网络气相饱和度，无量纲；sfw为裂缝网络水相饱和度，无量纲。
裂缝网络中存在毛管力，满足以下方程：
式中：pfc为裂缝网络中的毛管力，MPa；pfg为裂缝网络气相压力，MPa；pfw为裂缝网络水相压力，MPa。
步骤S3：
在本实施例中，建立无机质系统的气水两相流动方程具体步骤为：
页岩岩石基质中的纳米级孔隙大小通常在几纳米到几百纳米之间，导致页岩气体的流动方式不再是线性达西流动，因此需对其渗透率进行修正。Beskok 建立了考虑克努森扩散和滑脱效应的流动方程，其对达西公式修正为：
式中，k0为为储层的固有渗透率，m2；μ为气体黏度，Pa·s；Kn为克努森数；α、b为分别为稀薄系数和滑脱因子。其中，克努森数可由下式进行计算：
式中，λ为为气体分子自由行程，m；rn为为纳米级孔隙的喉道半径，m；KB为为Boltzmann常数，J·K-1；δ为为流动气体分子有效碰撞直径，m。其中，稀薄系数可以可由下式进行计算：
因此无机质中考虑克努森扩散和滑脱效应的渗透率可以作如下修正：
当地层为气水两相流动的时候，在毛管力的作用下水进入到无机质中甚至水会吸附在无机质纳米孔隙表面对渗透率产生影响，因此无机质孔隙的气水两相流动还需考虑吸附水层的影响，本方法采用曹成等人模型对考虑吸附水层的无机质孔喉半径和孔隙度进行修正。考虑吸附水层的无机质纳米孔隙模型如图2 所示，单根无机质纳米孔道长度为L，初始孔隙直径为Dc,0，含水饱和度为Sw的有效孔隙直径为Dc：
无机质系统孔隙直径Dc的含水饱和度表达式为：
对上式进一步化简可得：
当含水饱和度为sw时，有效孔隙度的表达式为：
将式子(1-12)、(1-13)代入Beskok渗透率修正模型即可得到考虑无机质中同时考虑吸附水层影响、克努森扩散和滑脱效应的修正渗透率。
由质量守恒定律可推导得到无机质系统中气相渗流微分方程：
式中，kc为页岩岩石无机质固有渗透率，m2；ζc为页岩岩石无机质渗透率修正系数，无量纲；kcrg为页岩岩石无机质中气相相对渗透率，无量纲；pcg为页岩岩石无机质中气相压力，MPa；qscfg为地面条件下单位体积页岩中无机质向裂缝的供给气量，sm3/s；qsocg为地面条件下单位体积页岩中有机质向无机质的供给气量，sm3/s；scg为页岩岩石无机质中气相饱和度，无量纲。
将窜流气量代入上式可得：
同理，根据质量守恒定律，可以得到无机质系统中水相渗流微分方程：
式中，kc为页岩岩石无机质固有渗透率，m2；kcrw为页岩岩石无机质中水相相对渗透率，无量纲；pcw为页岩岩石无机质中水相压力，MPa；qscfw为地面条件下单位体积页岩中无机质向裂缝的供给水量，sm3/s；scw为页岩岩石无机质中水相饱和度，无量纲。
将窜流水量代入上式可得：
无机质中气相饱和度和水相饱和度满足以下方程：
无机质中存在毛管力，满足以下方程：
式中，pcc为页岩岩石无机质中的毛管力，MPa。
步骤S4：
在本实施例中，建立有机质系统的气相流动方程具体步骤为：
当地层为气水两相流动的时候，水通常不会进入到亲油性较强的有机质里，因此有机质中无需考虑吸附水层对渗透率的影响，只需考虑页岩有机质纳米级孔隙中发生克努森扩散和滑脱现象，与无机质系统渗透率修正方法相同，这里不再赘述。
由质量守恒定律可推导得到有机质系统中气相渗流微分方程：
式中，ko为页岩岩石有机质固有渗透率，m2；ζo为页岩岩石有机质渗透率修正系数，无量纲；pog为页岩岩石有机质中气相压力，MPa；sog为页岩岩石有机质中气相饱和度，无量纲。
步骤S5：
在本实施例中，确定气水两相流动三孔三渗模型的初始条件及内外边界条件具体步骤为：
(1)初始条件
对于页岩气藏气水两相流动三孔三渗模型，满足如下初始条件：
式中，
为裂缝系统中初始气相压力，MPa；
为无机质系统初始气相压力，MPa；
为有机质系统初始气相压力，MPa；
为裂缝系统初始含水饱和度，无量纲；
为无机质系统初始含水饱和度，无量纲。
(2)外边界条件
本方法模型假设外边界封闭，所以压力在边界处满足如下关系式：
(3)内边界条件
本方法模型假设内边界为定井底流压，产量可采用径向流公式表示为：
其中：
式中：pgi,j,k为网格坐标为i,j,k的气体压力，MPa；δ(i,j,k)为Dirac函数，(存在源汇项时值为1反之为0)；pwf为井底流压，MPa；S为表皮系数，无量纲；qvg为产气量，sm3/s；qvw为产水量，sm3/s。
其中有关渗透率的坐标变换参考表1-1：
表1–1坐标变换表
步骤S6：
在本实施例中，利用有限差分法对数学模型进行离散化求得有限差分方程组，通过IMPES方法对其做线性化处理得到线性代数方程组，求解得到相应的数值模型具体步骤为：
采用有限差分法对数学模型进行离散化求得有限差分方程组，即非线性方程组，再用IMPES方法对其做线性化处理得到线性代数方程组，建立相应的数值模型，对页岩气藏气水两相渗流模型的求解采用IMPES方法，首先对压力方程组进行隐式求解，然后再对饱和度进行显示求解。
(1)裂缝系统
1)压力方程的推导
将裂缝系统的毛管力方程代入裂缝的两相渗流方程：
其中：
式中：λfw为基质与裂缝间水相窜流系数，无量纲；λfg为基质与裂缝间气相窜流系数，无量纲。
式子(1-27)和(1-28)右边项可作如下变换：
将式子(1-30)乘上气体体积系数Bg，式子(1-31)乘上水的体积系数 Bw，相加得到：
其中：
在这里，定义综合压缩系数Cft如下：
同样，将式子(1-27)左边乘上Bg，式子(1-28)左边乘上Bw，相加得到裂缝系统气相的压力方程为：
式子(1-35)参数具体定义如下：
2)隐式求解压力
首先对压力方程(1-35)两端乘上单元六面体体积VB＝ΔxiΔyjΔzk，然后采用有限差分法对其线性化处理。差分方程展开项冗长，在此以第一项为例进行展开。
对上式化简得到：
定义传导参数Tfl如下：
为了简化差分方程，可作如下定义：
其中：
得到裂缝系统气相压力的差分方程为：
其中：
最终可得到裂缝系统气相压力的线性方程组形式如下：
其中：
式中：N为求解网格中北向系数；S为求解网格中南向系数；W为求解网格中西向系数；E为求解网格中东向系数；C为求解网格中中心系数；A为求解网格中上表面系数；B为求解网格中下表面系数；Qi,j,k为源汇项；下标i,j,k表示网格坐标；上标n表示求解时间步。
传导系数的求解过程中，对与时间无关的部分T进行调和平均计算：
对fl进行上游权处理：
3)显示求解饱和度
当求解得到裂缝系统气相压力后，可由下式计算裂缝系统的气相饱和度：
则裂缝系统的水相饱和度为：
(2)无机质系统
1)压力方程的推导
将无机质系统的毛管力辅助方程代入无机质的两相渗流方程之中：
式中：
式子(1-59)和(1-60)右边项可作如下变换：
式子(1-62)乘上气体体积系数Bg，式子(1-63)乘上水的体积系数Bw，相加得到：
其中：
在这里，定义无机质系统的综合压缩系数Cct如下：
同样，将式子(1-59)左边乘上Bg，式子(1-60)左边乘上Bw，相加得到无机质系统气相的压力方程为：
式子(1-67)参数具体定义如下：
2)隐式求解压力
首先对无机质系统的压力方程(1-67)两端乘上单元六面体体积VB＝Δxi ΔyjΔzk，然后采用有限差分法对其线性化处理。差分过程与裂缝系统差分过程相同，这里不再详细叙述，得到无机质系统的气相压力差分方程如下：
式中：
最终可得到无机质系统中气相压力的线性方程组形式如下：
其中：
3)显示求解饱和度
求得到无机质系统气相压力后，可采用下式求解出无机质系统的气相饱和度：
则无机质系统的水相饱和度为：
(3)有机质系统
根据三孔三渗模型假设条件，有机质系统中只存在单相气体流动，因此可直接对有机质系统气相压力进行隐式求解：
上述式子右边项可作如下变换：
将式子(1-77)代入式子(1-76)，有机质系统气相的压力方程为：
其中：
接下来对有机质系统的压力方程(1-67)两端乘上单元六面体体积VB＝Δ xiΔyjΔzk，然后采用有限差分法对其线性化处理，差分过程不再详细叙述，得到有机质系统压力的差分方程如下：
最终可得到有机质系统中气相压力的线性方程组形式如下：
其中：
通过对上式求解，则可得到有机质系统中的气相压力。
步骤S7：
在本实施例中对页岩气井进行点源或点汇数学处理，采用定流压方式生产，利用Matlab对多级压裂水平井气水两相流动三孔三渗数值模型进行编程求解，分析压力动态变化具体步骤为：
(1)井处理
网格内若有井存在，通常把井视为一点源或点汇进行数学处理，根据数模理论，网格内井的产量可以用拟稳态流公式表示：
则井的生产指数可定义为：
式中：Δn为在n方向的网格步长，m；pwf为井底流压，MPa；PID为井的生产指数；krl为l相流体的相对渗透率；re为井点处网格块的等效半径，m；rw为井筒半径，m；S为表皮因子；Δn为在n方向的网格步长，m。
(2)定流压生产
当井以定井底流压生产时，第k完井段的产气量可表示为：
则第k完井段的产水量为：
(3)计算机求解与模型结果分析
基于前述有限差分法求得的页岩气藏气水两相渗流数值模型，下面对模型进行求解。由于线性方程组的计算量庞大，因此本节将采用Matlab软件对其进行编程求解，程序详细求解流程如图3所示：
下面将进一步对本方法建立的页岩气藏气水两相渗流模型进行模拟验证，首先建立地质网格模型如图4所示：
参考实际页岩气藏地层参数，详细参数设置如表1所示：
表1模拟验证参数设置表
对于页岩基质系统和裂缝系统存在毛管力，本方法采用Brooks等提出的半经验公式对毛管力进行计算：
式中，sew为归一化的含水饱和度，无量纲；pe为排驱压力，MPa；λ为孔隙大小分布量度，无量纲。
归一化的含水饱和度表示为：
式中，swc为束缚水饱和度，无量纲。
根据上式则可绘制出基质和裂缝系统中气水两相流动的毛管力曲线，如图5 和图6所示：
本方法采用Li等人提出的方法对基质和裂缝系统的相对渗透率进行计算。
其中，水相相对渗透率可由下式计算：
同样，气相相对渗透率可由下式计算：
根据式子(1-89)、(1-90)可计算得到页岩气藏气水两相流动时的相渗曲线如图7所示：
根据上述模拟参数，通过Matlab软件编程输出三孔三渗模型网格气相压力的变化情况如图8-11所示。
由图8-图11可以观察到，当模拟至第10天时，主裂缝压力明显下降，在图中可观察到压降漏斗，同时可发现网格裂缝压力迅速下降，无机质压力下降缓慢，有机质压力下降最为缓慢。随着生产的进行，当模拟至第300天时，压力进一步传播到边界处，裂缝压力持续降低，同时无机质参与能量供给进而产生压力下降现象，有机质压力较之前有所下降但下降趋势较无机质压力缓慢。当模拟至第600天时，压力传播到边界处，所有网格裂缝压力均持续下降，同时边界处的无机质压力也出现了明显下降现象，有机质压力下降缓慢。通过模拟结果发现，早期裂缝中的天然气先被采出，且裂缝网络作为产量供给的主要通道，因此整个过程裂缝压力下降迅速，随着生产的进行，无机质中的游离天然气开始流向裂缝同时着吸附气进一步解吸为游离气参与流动，无机质中压力在中后期下降较明显，当生产压差进一步扩大，有机质中的游离气和吸附气也开始参与流动，因此中后期也能观察到有机质压力下降现象。通过分析生产过程中的压力动态变化，进一步验证了本方法所建立的页岩气藏气水两相三孔三渗模型的准确性。
在本实施例中，(1)本方法将有机质单独提出来视为一个渗流介质，考虑页岩气藏由物理性质独立的有机质系统、无机质系统和裂缝系统组成，建立了页岩气藏气水两相流动三孔三渗模型，该模型对有机质和无机质分开建立渗流方程，假设气体先从有机质系统流向无机质系统，再从无机质系统流向天然裂缝，无机质系统、天然裂缝系统中为气水两相流动，有机质系统为单相气体流动，同时考虑了无机质纳米级孔隙吸附水层的影响对渗透率进行了修正，综合考虑页岩气流动为压力差作用下的黏性流动、有机质和黏土矿物的吸附解吸、滑脱效应和克努森扩散多重机制的作用，推导得到了页岩气藏气水两相流动三孔三渗模型。
(2)采用有限差分法求解页岩气藏气水两相流动三孔三渗模型，通过IMPES 方法，结合Matlab软件对模型进行编程求解，最后结合实际页岩气田的地层参数和生产制度对压力和产气量进行了动态预测。
(3)本方法将页岩基质中有机质和无机质的气水两相流动规律分开考虑，能准确描述页岩气藏气水两相流动规律，更能准确预测页岩气井的压力动态变化和气井生产能力。
本发明所揭露的方法、系统和模块，可以通过其它的方式实现。例如，以上所描述的实施例仅是示意性的，例如，所述模块的划分，可以仅仅是一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个模块或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以说通过一些接口，系统或模块的间接耦合或通信连接，可以是电性，机械或其它的形式。
另外，在本发明各个实施例中的各功能模块可以集成在一个处理模块中，也可以是各个模块单独物理存在，也可以两个或两个以上模块集成在一个模块中。
以上所述仅是本发明的优选实施方式，应当理解本发明并非局限于本文所披露的形式，不应看作是对其他实施例的排除，而可用于各种其他组合、修改和环境，并能够在本文所述构想范围内，通过上述教导或相关领域的技术或知识进行改动。而本领域人员所进行的改动和变化不脱离本发明的精神和范围，则都应在本发明所附权利要求的保护范围内。
Claims (10)
1.一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，其特征在于，所述三孔三渗模型构建及压力动态预测方法包括如下步骤：
构建气水两相流动三孔三渗模型：
S101：分别构建裂缝系统、无机质系统的气水两相流动方程、毛管力辅助方程和饱和度辅助方程以及有机质系统的气相流动方程；
S102：基于获取的裂缝系统、无机质系统和有机质系统的初始条件及内外边界条件，将饱和度辅助方程、毛管力辅助方程代入气水两相流动方程，分别得到裂缝系统、无机质系统和有机质系统下的气相压力方程；
S103：利用有限差分法得到气相压力的有限差分方程，通过IMPES方法对有限差分方程做线性化处理得到气相压力的线性代数方程组，将该线性代数方程组作为水平井气水两相流动三孔三渗模型；
压力动态预测：
S201：对存在于裂缝系统网络、无机质系统系统和有机质系统网络中的页岩气井进行点源或点汇数学处理，将网格内井的产量用拟稳态流公式表示得到井的生产指数；
S202：采用定流压方式生产，利用MATLAB对多级压裂水平井气水两相流动双孔双渗数值模型进行编程求解，获取该井以定井底流压生产时的产气量与产水量。
2.如权利要求1所述的一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，其特征在于，所述构建裂缝系统的气水两相流动方程、毛管辅助方程和饱和度辅助方程包括：
A1：将单位体积页岩无机质向裂缝窜流水量表达式带入裂缝网络中气相渗流微分方程，得到裂缝系统下的气相方程：
式中：▽为Hamilton算子；kf为裂缝系统渗透率，m2；kfrg为裂缝网络气相相对渗透率，无量纲；Bg为气体的体积系数，无量纲；μg为气体粘度，mPa·s；pfg为裂缝网络气相压力，MPa；ρfg为气体密度，kg/m3；D为气体所处的垂直深度，m；qsg为地面条件下单位体积页岩中采出的气量，sm3/s；α为形状因子，1/m2；ξc为页岩岩石无机质渗透率修正系数，无量纲；kcrg为无机质网络气相相对渗透率，无量纲；pfg为裂缝系统气相压力，MPa；φf为裂缝系统孔隙度，无量纲；sfg为裂缝网络气相饱和度，无量纲；
A2：将单位体积页岩无机质向裂缝窜流水量的表达式带入裂缝网络中水相渗流方程，得到裂缝系统下的水相方程：
式中：kfrw为裂缝网络水相相对渗透率，无量纲；Bw为水相的体积系数，无量纲；μw为水相粘度，mPa·s；pfw为裂缝系统水相压力，MPa；ρfw为水相密度，kg/m3；qsw为地面条件下单位体积页岩中采出的水量，sm3/s；kc为无机质系统渗透率，m2；kcrw为基质网络水相相对渗透率，无量纲；Bw为水相的体积系数，无量纲；μw为水相粘度，mPa·s；pcw为无机质网络水相压力，MPa；pfw为裂缝网络水相压力，MPa；sfw为裂缝网络水相饱和度，无量纲。
A3：获取裂缝网络中的气相饱和度和水相饱和度，得到饱和度辅助方程：
sfg+sfw＝1；
式中：sfg为裂缝网络气相饱和度，无量纲；sfw为裂缝网络水相饱和度，无量纲。
A4：获取裂缝网络中的毛管力，得到毛管辅助方程：
pfc(sfw)＝pfg-pfw；
式中：pfc为裂缝网络中的毛管力，MPa；pfg为裂缝网络气相压力，MPa；pfw为裂缝网络水相压力，MPa。
3.如权利要求2所述的一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，其特征在于，所述构建基质系统的气水两相流动方程、毛管辅助方程和饱和度辅助方程包括：
B1：将单位体积页岩无机质向裂缝窜流气量表达式带入无机质网络中气相渗流微分方程，得到基质系统下的气相方程：
式中：kc为无机质固有渗透率，m2；ζc为页岩岩石无机质渗透率修正系数，无量纲；kcrg为无机质网络气相相对渗透率，无量纲；pcg为页岩岩石无机质中气相压力，MPa；Smg为基质网络中气相饱和度，无量纲；φm为基质孔隙度；VL为基质的兰氏体积，sm3/m3；pL为有基质的兰氏压力，MPa。
B2：将单位体积页岩无机质向裂缝窜流水量的表达式带入无机质网络中水相渗流方程，得到基质系统下的水相方程：
式中：kc为页岩岩石无机质固有渗透率，m2；kcrw为页岩岩石无机质中水相相对渗透率，无量纲；pcw为页岩岩石无机质中水相压力，MPa；qscfw为地面条件下单位体积页岩中无机质向裂缝的供给水量，sm3/s；scw为页岩岩石无机质中水相饱和度，无量纲；
B3：获取无机质网络中的气相饱和度和水相饱和度，得到饱和度辅助方程：
scg+scw＝1；
B4：获取无机质网络中的毛管力，得到毛管辅助方程：
pcc(scw)＝pcg-pcw；
式中：pcc为页岩岩石无机质中的毛管力，MPa。
4.如权利要求3所述的一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，其特征在于，构建有机质系统的气相流动方程通过质量守恒定律推导得到有机质系统中气相渗流微分方程：
式中：ko为页岩岩石有机质固有渗透率，m2；ζo为页岩岩石有机质渗透率修正系数，无量纲；pog为页岩岩石有机质中气相压力，MPa；sog为页岩岩石有机质中气相饱和度，无量纲。
5.如权利要求4所述的一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，其特征在于，所述获取的裂缝系统、无机质系统和有机质系统的初始条件及内外边界条件包括：
C1：裂缝系统和基质系统下的初始气相压力、初始含水饱合度：
式中：
为裂缝系统中初始气相压力，MPa；
为无机质系统初始气相压力，MPa；
为有机质系统初始气相压力，MPa；
为裂缝系统初始含水饱和度，无量纲；
为无机质系统初始含水饱和度，无量纲；
C2：页岩气藏三孔三渗模型的外边界条件：
C3：页岩气藏三孔三渗模型的内边界条件：
上式中：
式中：pgi,j,k为网格坐标为i,j,k的气体压力，MPa；δ(i,j,k)为Dirac函数，(存在源汇项时值为1反之为0)；pwf为井底流压，MPa；S为表皮系数，无量纲；qvg为产气量，sm3/s；qvw为产水量，sm3/s。
6.如权利要求5所述的一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，其特征在于，所述得到裂缝系统下的气相压力方程，通过将毛管力辅助方程和饱和度辅助方程代入气水两相渗流方程，得到只含有气相压力的方程：
D1：将毛管力辅助方程代入三孔三渗模型中裂缝系统的两相渗流方程之中：
其中：
D2：将饱和度辅助方程代入上述两相渗流方程式子，并分别乘以气体体积系数Bg和水体积系数Bw，相加得到裂缝系统气相的压力方程：
其中：
7.如权利要求6所述的一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，其特征在于，利用有限差分法得到气相压力的有限差分方程，通过IMPES方法对有限差分方程做线性化处理得到气相压力的线性代数方程组：
对三孔三渗模型的裂缝系统压力隐式求解，得到压力线性方程组：
其中：
式中：N为求解网格中北向系数；S为求解网格中南向系数；W为求解网格中西向系数；E为求解网格中东向系数；C为求解网格中中心系数；A为求解网格中上表面系数；B为求解网格中下表面系数；Qi,j,k为源汇项；下标i,j,k表示网格坐标；上标n表示求解时间步。
8.如权利要求7所述的一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，其特征在于，所述得到无机质系统下的气相压力方程，通过将毛管力辅助方程和饱和度辅助方程代入气水两相渗流方程，得到只含有气相压力的方程：
E1：将毛管力辅助方程代入三孔三渗模型中无机质系统的两相渗流方程之中：
其中：
E2：将上述两相渗流方程式子代入饱和度辅助方程，并分别乘以气体体积系数Bg和水体积系数Bw，相加得到无机质系统气相的压力方程：
其中：
9.如权利要求8所述的一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，其特征在于，利用有限差分法得到气相压力的有限差分方程，通过IMPES方法对有限差分方程做线性化处理得到气相压力的线性代数方程组：
对三孔三渗模型的无机质系统压力隐式求解，得到压力线性方程组：
其中：
10.如权利要求9所述的一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法，其特征在于，所述利用有限差分法得到气相压力的有限差分方程，通过IMPES方法对有限差分方程做线性化处理得到气相压力的线性代数方程组：
对三孔三渗模型的有机质系统压力隐式求解，得到压力线性方程组：
其中：
所述对存在于裂缝系统网络、无机质系统网络和有机质系统网络中的页岩气井进行点源或点汇数学处理，将网格内井的产量用拟稳态流公式表示得到井的生产指数，其中：
拟稳态流公式为：
井的生产指数为：
式中：Δn为在n方向的网格步长，m；pwf为井底流压，MPa；PID为井的生产指数；krl为l相流体的相对渗透率；re为井点处网格块的等效半径，m；rw为井筒半径，m；S为表皮因子；Δn为在n方向的网格步长，m。
所述采用定流压方式生产，利用MATLAB对多级压裂水平井气水两相流动双孔双渗数值模型进行编程求解，获取该井以定井底流压生产时的产气量与产水量，其中：
在定流压生产时，第k完井段的产气量可表示为：
第k完井段的产水量为：
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010973613.9A CN112012731B (zh) | 2020-09-16 | 2020-09-16 | 一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010973613.9A CN112012731B (zh) | 2020-09-16 | 2020-09-16 | 一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112012731A true CN112012731A (zh) | 2020-12-01 |
| CN112012731B CN112012731B (zh) | 2021-08-31 |
Family
ID=73522020
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010973613.9A Active CN112012731B (zh) | 2020-09-16 | 2020-09-16 | 一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112012731B (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116335654A (zh) * | 2023-05-17 | 2023-06-27 | 重庆科技学院 | 一种模拟页岩气特殊机理的压裂水平井产量预测方法 |
| US11905807B2 (en) | 2021-06-30 | 2024-02-20 | Saudi Arabian Oil Company | Workflow to model and condition water distribution in a hydrocarbon reservoir |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130346040A1 (en) * | 2012-06-20 | 2013-12-26 | Schlumberger Technology Corporation | Shale gas production forecasting |
| CN103792581A (zh) * | 2013-12-31 | 2014-05-14 | 刘致水 | 基于岩石物理模型的页岩气储层地应力测井预测方法 |
| CN106337677A (zh) * | 2016-10-17 | 2017-01-18 | 长江大学 | 页岩气压裂缝网气水两相导流能力测试系统及测试方法 |
| CN108133080A (zh) * | 2017-12-04 | 2018-06-08 | 陕西延长石油(集团)有限责任公司研究院 | 一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法 |
| CN109408859A (zh) * | 2018-09-05 | 2019-03-01 | 中国石油集团川庆钻探工程有限公司 | 页岩气藏压裂水平井二维三重介质数值模型建立方法 |
| CN109815516A (zh) * | 2018-09-10 | 2019-05-28 | 中国石油天然气股份有限公司 | 对页岩气井产能进行预测的方法及装置 |
-
2020
- 2020-09-16 CN CN202010973613.9A patent/CN112012731B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130346040A1 (en) * | 2012-06-20 | 2013-12-26 | Schlumberger Technology Corporation | Shale gas production forecasting |
| CN103792581A (zh) * | 2013-12-31 | 2014-05-14 | 刘致水 | 基于岩石物理模型的页岩气储层地应力测井预测方法 |
| CN106337677A (zh) * | 2016-10-17 | 2017-01-18 | 长江大学 | 页岩气压裂缝网气水两相导流能力测试系统及测试方法 |
| CN108133080A (zh) * | 2017-12-04 | 2018-06-08 | 陕西延长石油(集团)有限责任公司研究院 | 一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法 |
| CN109408859A (zh) * | 2018-09-05 | 2019-03-01 | 中国石油集团川庆钻探工程有限公司 | 页岩气藏压裂水平井二维三重介质数值模型建立方法 |
| CN109815516A (zh) * | 2018-09-10 | 2019-05-28 | 中国石油天然气股份有限公司 | 对页岩气井产能进行预测的方法及装置 |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11905807B2 (en) | 2021-06-30 | 2024-02-20 | Saudi Arabian Oil Company | Workflow to model and condition water distribution in a hydrocarbon reservoir |
| CN116335654A (zh) * | 2023-05-17 | 2023-06-27 | 重庆科技学院 | 一种模拟页岩气特殊机理的压裂水平井产量预测方法 |
| CN116335654B (zh) * | 2023-05-17 | 2024-01-30 | 重庆科技学院 | 一种模拟页岩气特殊机理的压裂水平井产量预测方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112012731B (zh) | 2021-08-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106484933B (zh) | 一种用于确定页岩气井井控动态储量的方法及系统 | |
| Wang et al. | Spontaneous imbibition analysis in shale reservoirs based on pore network modeling | |
| CN107622165A (zh) | 一种页岩气水平井重复压裂产能计算方法 | |
| CN108133080B (zh) | 一种考虑非稳态吸附的非均质裂缝页岩气藏数值模拟方法 | |
| Zhang et al. | Employing a quad-porosity numerical model to analyze the productivity of shale gas reservoir | |
| CN112012731B (zh) | 一种基于气水两相流动页岩气藏三孔三渗模型构建及压力动态预测方法 | |
| CN112084718B (zh) | 基于渗流差异的页岩气藏单相气体三孔三渗模型构建方法 | |
| Tian et al. | A compositional model for gas injection IOR/EOR in tight oil reservoirs under coupled nanopore confinement and geomechanics effects | |
| CN106547938B (zh) | 裂隙-孔隙结构双重介质煤储层气水两相流数值模拟方法 | |
| Cui et al. | Impact of shale matrix mechanical interactions on gas transport during production | |
| CN107622328A (zh) | 一种页岩气藏分段多簇压裂水平井产量预测方法 | |
| Li et al. | An analytical method for modeling and analysis gas-water relative permeability in nanoscale pores with interfacial effects | |
| CN106481332A (zh) | 用于确定页岩气多段压裂水平井内外区动态储量的方法 | |
| CN112360422B (zh) | 一种页岩气藏压裂水平井产量预测方法及系统 | |
| Wu et al. | Pressure transient analysis of multiple fractured horizontal well in composite shale gas reservoirs by boundary element method | |
| CN111980695B (zh) | 基于不同占比的页岩气藏单相气体双孔双渗模型构建方法 | |
| CN112081583A (zh) | 非常规储层缝网压裂多尺度支撑产量计算方法及装置 | |
| AU2020103953A4 (en) | Method and system for predicting production of fractured horizontal well in shale gas reservoir | |
| CN110472348A (zh) | 一种页岩气藏非稳态渗流模型的建立方法 | |
| CN114169263A (zh) | 基于压裂液返排数据的页岩油气藏裂缝参数反演方法 | |
| Tian et al. | Performance evaluation of gas production with consideration of dynamic capillary pressure in tight sandstone reservoirs | |
| Huang et al. | A semi-analytical model to evaluate productivity of shale gas wells with complex fracture networks | |
| CN111963163B (zh) | 一种基于气水两相流动页岩气藏双孔双渗模型构建及压力动态预测方法 | |
| Le et al. | Bridging between macroscopic behavior of shale gas reservoirs and confined fluids in nanopores | |
| Su et al. | Relative permeability estimation of oil− water two-phase flow in shale reservoir |
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 | ||
| OL01 | Intention to license declared | ||
| OL01 | Intention to license declared |