CN113128040B - 一种基于孔隙结构特征的页岩表观渗透率预测方法 - Google Patents

一种基于孔隙结构特征的页岩表观渗透率预测方法 Download PDF

Info

Publication number
CN113128040B
CN113128040B CN202110389015.1A CN202110389015A CN113128040B CN 113128040 B CN113128040 B CN 113128040B CN 202110389015 A CN202110389015 A CN 202110389015A CN 113128040 B CN113128040 B CN 113128040B
Authority
CN
China
Prior art keywords
apparent permeability
pore
gas
shale
porous medium
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
Application number
CN202110389015.1A
Other languages
English (en)
Other versions
CN113128040A (zh
Inventor
薛培
赵谦平
祁攀文
曹成
张建锋
史鹏
孙建博
尹锦涛
刘刚
王晖
徐杰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shaanxi Yanchang Petroleum Group Co Ltd
Original Assignee
Shaanxi Yanchang Petroleum Group Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Shaanxi Yanchang Petroleum Group Co Ltd filed Critical Shaanxi Yanchang Petroleum Group Co Ltd
Priority to CN202110389015.1A priority Critical patent/CN113128040B/zh
Publication of CN113128040A publication Critical patent/CN113128040A/zh
Application granted granted Critical
Publication of CN113128040B publication Critical patent/CN113128040B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • 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
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
    • G01N15/08—Investigating permeability, pore-volume, or surface area of porous materials
    • G01N15/082—Investigating permeability by forcing a fluid through a sample
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • 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
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N15/00—Investigating characteristics of particles; Investigating permeability, pore-volume or surface-area of porous materials
    • G01N15/08—Investigating permeability, pore-volume, or surface area of porous materials
    • G01N2015/0866—Sorption
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
    • Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Business, Economics & Management (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Algebra (AREA)
  • Chemical & Material Sciences (AREA)
  • Strategic Management (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Human Resources & Organizations (AREA)
  • Economics (AREA)
  • Evolutionary Biology (AREA)
  • Quality & Reliability (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Hardware Design (AREA)
  • Development Economics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Game Theory and Decision Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Geometry (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Marketing (AREA)
  • Evolutionary Computation (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Computing Systems (AREA)
  • Dispersion Chemistry (AREA)
  • Health & Medical Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)

Abstract

本发明涉及一种基于孔隙结构特征的页岩表观渗透率预测方法。一种基于孔隙结构特征的页岩表观渗透率预测方法，依次建立微孔道表观渗透率模型及建立多孔介质表观渗透率模型。本发明建立的一种基于孔隙结构特征的页岩表观渗透率预测方法可用于不同埋深、储层温压条件下，页岩气开发过程中的页岩表观渗透率的动态预测。考虑了随着孔隙压力改变，不同孔隙直径的微孔道内气体传输形式的变化特征，考虑了有效应力及基质收缩对页岩表观渗透率的影响，明确提出了多孔介质内微孔道表观渗透率分配系数的理论推导过程，实现了采用气体吸附法实验获取的页岩孔隙结构实测数据计算表观渗透率分配系数，进而实现了基于孔隙结构特征的页岩表观渗透率预测。

Description

一种基于孔隙结构特征的页岩表观渗透率预测方法
技术领域
本发明涉及一种基于孔隙结构特征的页岩表观渗透率预测方法。
背景技术
页岩气是一种潜在资源量巨大的非常规天然气，随着勘探开发的深入，对页岩孔隙结构特征及渗流机理研究日益受到重视，同时取得了很大的进步。Reed、Sondergeld、Curtis等发现页岩中存在大量的纳米级孔隙。Loucks等研究了Barnett页岩中的纳微米孔隙，按照孔隙尺度将孔隙划分为微米孔(直径大于0.75μm)和纳米孔(直径小于0.75μm)。聂海宽等认为页岩主要储集类型为裂缝和孔隙，孔裂隙尺度处于纳米及微米级别，巨、大型裂缝为毫米级。由于页岩孔隙结构的多尺度特征，导致页岩气传输方式与常规气藏不同，页岩中气体传输方式呈现多样性。Wang和Reed指出页岩天然裂缝及人工裂缝中气体传输方式为达西流，纳米级孔隙内为非达西流。Schepers等指出页岩气在纳米级孔隙和微裂缝中流动方式以非达西渗流为主。关富佳等指出页岩气藏孔隙结构的多尺度决定了渗流方式的多样性。
目前吴克柳、王瑞、Javadpour F、段永刚等分别建立了微孔道内气体传输表观渗透模型，同时分析了微孔道内气体传输能力，但此类模型只能表征具有特定孔隙直径的微孔道的气体传输能力，尚且不能预测具有一定孔隙结构的多孔介质表观渗透率，无法应用于实际页岩储层的表观渗透率预测。
发明内容
本发明旨在针对上述问题，提出一种基于孔隙结构特征的页岩表观渗透率预测方法。该方法在充分考虑不同气体传输方式的基础上，采用分段函数的形式，建立了微孔道表观渗透率预测方法，分析了微孔道内气体传输特征；同时结合多孔介质孔隙结构特征，建立了多孔介质表观渗透率预测方法。
本发明的技术方案在于：
(一)本发明提出一种微孔道表观渗透率模型
气体传输方式由孔道特征及气体特性共同决定，普遍使用克努森数作为气体传输方式的判别标准。克努森数定义为气体分子平均自由程与孔道直径的比值，其表达式为：
式中：Kn为克努森数，量纲1；
为气体分子平均自由程，m；
kb为玻尔兹曼常数，1.38×10-23J·K-1；
T为绝对温度，K；
σ为气体分子平均碰撞直径，m；
p为孔隙压力，Pa；
d为孔道直径，m；
气体分子平均自由程与温度呈正相关，与孔隙压力、气体分子平均碰撞直径呈负相关。表1所示为CH4、C2H6、CO2单组分气体及多组分气体分子平均碰撞直径。多组分气体的平均分子碰撞直径以单组分气体摩尔百分比为权重，对各单组分气体平均分子碰撞直径进行加权求和而得；
表1典型页岩气相关参数
采用克努森数为0.001、0.1、10作为气体传输方式的划分界限，将气体传输方式划分为达西流、滑脱流、过渡流和努森扩散，如表2所示；
表2基于克努森数的气体传输形式划分
1.1努森扩散
当克努森数大于10，气体分子在孔道内的传输方式为自由分子流即努森扩散。气体分子努森扩散摩尔通量为
式中：Jkn为努森扩散摩尔通量，mol/(m2·s)；
R为通用气体常数，8.314J/(mol·K)；
M为气体摩尔质量，kg/mol；
C为气体摩尔浓度，mol/m3；
l为气体传输方向的距离，m；
对于真实气体而言，其状态方程为
pVmol＝ZRT式(3)
式中：p为气体压力，Pa；
Vmol为气体摩尔体积，m3/mol；
T为绝对温度，K；
Z为真实气体压缩因子，量纲1，使用Hall-Yarborough法计算；
由此可知
则
1.2过渡流
当克努森数为0.1～10之间，气体分子在孔道内的传输方式为过渡流。过渡流域内气体传输能力可由努森扩散与滑脱流动综合表征，即
Jtr＝θJsl+(1)Jkn式(6)
式中：Jtr为过渡流摩尔通量，mol/(m2·s)；
Jsl为滑脱流动摩尔通量，mol/(m2·s)；
θ为过渡流中存在滑脱现象的分子数占总分子数的比例系数，量纲1；Kn0.、n为拟合常数，Kn0.＝4.5，n＝5；图1所示为θ与Kn的关系曲线。当Kn低于0.1，θ趋近于1，表明此Kn范围内，气体传输方式中不存在努森扩散；当Kn大于15，θ趋近于0，表明此Kn范围内，气体传输方式中不存在滑脱流、达西流；
1.3滑脱流
当克努森数为0.001～0.1之间，气体分子在孔道内的传输方式为滑脱流动。滑脱流动中气体分子与孔道壁面间的碰撞不可忽略，气体分子在孔道壁面运动速度不为零，需使用滑脱边界条件对Hagen-Poiseuille方程进行修正；
Jsl＝δκJda式(7)
式中：Jsl为滑脱流摩尔通量，mol/(m2·s)；
δκ为滑脱修正系数，量纲1；
式中：μ为气体粘度，Pa·s；
r为孔道半径，m；
k为切向动量适应系数，k＝0.8，量纲1；
由于气体分子平均自由程λ亦可表示为则式(8)可以改写为
气体粘度μ计算使用Lee-Gonzden-Enkin半经验法；
Jda为达西流摩尔通量，mol/(m2·s)；
1.4达西流
当克努森数小于0.001，气体分子在孔道内的传输方式为达西流动。气体分子做达西流的达西流摩尔通量表示为
式中：Jda为达西流摩尔通量，mol/(m2·s)；
1.5微孔道表观渗透率模型
孔道内气体传输方式由温度、压力、气体类型、孔道直径及孔隙压力决定，气体传输能力可由相应气体传输方式的摩尔通量表征。不同气体传输方式的摩尔通量表达式中已考虑了孔道内气体分子之间、气体分子与壁面的两种碰撞类型，若采用达西流摩尔通量与努森扩散摩尔通量的相应占比求和的方式表征微孔道气体传输摩尔通量则会造成误差。因此采用分段函数的形式表征一定温度、压力、气体类型及孔隙直径的孔道内气体传输能力，此分段函数的分段界限为克努森数，即
式中：Jd为具有特定孔径的微孔道内气体传输摩尔通量，mol/(m2·s)；
对于多孔介质而言，假设其迂曲度是均匀的，则多孔介质的微孔道内气体传输摩尔通量Jdpm可表示为
式中：Jdpm为多孔介质的微孔道内气体传输摩尔通量，mol/(m2·s)；
τ为迂曲度，量纲1，迂曲度τ与孔隙度满足/>由此可得，多孔介质的微孔道内表观渗透率表达式为
式中：Kapp,d表示多孔介质的微孔道内气体传输表观渗透率，m2；
ρ为气体密度，Kg/m3；
结合式(5)、式(6)、式(7)、式(9)、式(10)、式(11)、式(12)、式(13)可得多孔介质的微孔道内表观渗透率为
式中：
(二)本发明提出一种多孔介质表观渗透率模型
上述(一)中以上所建立的微孔道表观渗透率模型是针对具有特定孔隙直径的单个微孔道而言。对于多孔介质而言，假设气体传输通道是由若干具有一定迂曲度及孔径分布的微孔道组成，微孔道之间由不具气体通过能力的物质组成，即毛管束模型。具体假设多孔介质中微孔道的孔隙直径分别为d1，d2，……，di，相应孔隙直径对应的微孔道数目分别为N1，N2，……，Ni。以下讨论内容针对具单位质量的多孔介质而言；具有特定孔隙直径的微孔道表观渗透率表示为Kapp,di，则该微孔道气体流速可按照达西定律表示为
式中为具有特定孔隙直径di的微孔道内气体流速，m/s；
该微孔道气体流量可表示为
式中：为具有特定孔隙直径di的微孔道气体流量，m3/s；
为具有特定孔隙直径di的微孔道截面积，m2；
则具有相同孔隙直径di的微孔道集合的总气体流量为
式中：为具相同孔隙直径di的微孔道集合的总气体流量，m3/s；
Ni为具相同孔隙直径di的微孔道的数量，整数；
多孔介质气体流量为所有微孔道气体流量的综合，可表示为
式中：qt为多孔介质气体流量，m3/s；
对于多孔介质而言，其表观渗透率为Kapp，则多孔介质气体流量qt同时可以表示为
式中：Kapp为多孔介质表观渗透率，m2；
S为多孔介质截面积，m2；
联立式(15)，式(16)，式(17)，式(18)，式(19)，可得
式中：表示具有孔隙直径di的微孔道的总截面积，m2；
表示具有孔隙直径di的微孔道的总截面积与多孔介质截面积比，小数；
定义f(di)为多孔介质内微孔道表观渗透率分配系数，量纲1，令则多孔介质表观渗透率为
由于假设多孔介质的孔道为平行毛细管束结构，因此孔道的截面积占比亦可表示为孔道的体积占比，则fi)可表示为
式中：为多孔介质孔隙度，％；
为具孔隙直径di的微孔道体积，cm3；
表示多孔介质总孔体积，cm3；
式(22)中涉及的孔隙体积、孔隙直径数据，可通过气体吸附法实验获取，且测点数量越多，孔隙结构刻画越精细，相应的表观渗透率预测结果的准确度越高。同时由于在微孔道表观渗透率中已经考虑了多孔介质孔道的迂曲度等(见式12)，因此此处不再考虑。
(三)本发明提出一种考虑有效应力影响下的多孔介质表观渗透率模型对于实际储层而言，产气过程中，上覆地层压力保持恒定，随着孔隙压力的降低，储层所承受的有效应力随之增加，造成孔道压缩，储层渗透率下降。考虑有效应力条件下的多孔介质表观渗透率可表示为：
式中，σeff＝pc-αp，σeff为有效应力，MPa；
α为Biot常数，量纲1；
pc为上覆地层压力，MPa；
p1为多孔介质孔隙闭合时所承受的有效应力，MPa；当有效应力σeff达到p1，此时Kapp,＝0；
β为与孔道表面粗糙度有关的常数，小数；
K0为有效应力为0时的表观渗透率，m2；
由于式(21)所表征的多孔介质表观渗透率未考虑有效应力的影响，因此可视K0＝Kapp。则考虑有效应力影响的多孔介质表观渗透率可表示为
(四)本发明提出一种考虑有效应力及基质收缩影响下的多孔介质表观渗透率模型在页岩储层降压产气过程中，随着孔隙压力的降低，孔隙表面的吸附气将解吸为游离气，页岩基质发生收缩，孔道尺寸变大。基质收缩效应将提升孔道的气体传输能力。采用基质收缩渗透率修正系数作为基质收缩对多孔介质表观渗透率的修正系数。基质收缩渗透率修正系数为
式中：ωs为页岩有机质收缩效应渗透率修正系数，量纲1；
pint为初始孔隙压力，MPa；
εL为Langmuir应变，量纲1；
pL为Langmuir压力，MPa；
因此考虑有效应力及基质收缩影响的多孔介质表观渗透率可表示为
本发明的技术效果在于：
本发明建立的一种基于孔隙结构特征的页岩表观渗透率预测方法可用于不同埋深、储层温压条件下，页岩气开发过程中的页岩表观渗透率的动态预测。该方法充分考虑了随着孔隙压力改变，不同孔隙直径的微孔道内气体传输形式的变化特征，考虑了有效应力及基质收缩对页岩表观渗透率的影响，明确提出了多孔介质内微孔道表观渗透率分配系数的理论推导过程，实现了采用气体吸附法实验获取的页岩孔隙结构实测数据计算表观渗透率分配系数，进而实现了基于孔隙结构特征的页岩表观渗透率预测。经验证，该方法对开发过程中页岩表观渗透动态变化的预测精度高。
附图说明
图1为实施例1中θ与Kn的关系曲线。
图2为实施例1中微孔道内孔隙压力与表观渗透率关系曲线。
图3为实施例1中微孔道内克努森数与表观渗透率关系曲线。
图4为实施例2页岩孔径分布曲线。
图5为实施例2方程拟合结果。
图6为实施例3模拟结果。
具体实施方式
实施例1—微孔道内气体传输表观渗透率模型
气体传输方式由孔道特征及气体特性共同决定，普遍使用克努森数作为气体传输方式的判别标准。克努森数定义为气体分子平均自由程与孔道直径的比值，其表达式为：
式中：Kn为克努森数，量纲1；
为气体分子平均自由程，m；
kb为玻尔兹曼常数，1.38×10-23J·K-1；
T为绝对温度，K；
σ为气体分子平均碰撞直径，m；
p为孔隙压力，Pa；
d为孔道直径，m；
气体分子平均自由程与温度呈正相关，与孔隙压力、气体分子平均碰撞直径呈负相关。表1所示为CH4、C2H6、CO2单组分气体及多组分气体分子平均碰撞直径。多组分气体的平均分子碰撞直径以单组分气体摩尔百分比为权重，对各单组分气体平均分子碰撞直径进行加权求和而得；
表1典型页岩气相关参数
采用克努森数为0.001、0.1、10作为气体传输方式的划分界限，将气体传输方式划分为达西流、滑脱流、过渡流和努森扩散，如表2所示；
表2基于克努森数的气体传输形式划分
1.1努森扩散
当克努森数大于10，气体分子在孔道内的传输方式为自由分子流即努森扩散。气体分子努森扩散摩尔通量为
式中：Jkn为努森扩散摩尔通量，mol/(m2·s)；
R为通用气体常数，8.314J/(mol·K)；
M为气体摩尔质量，kg/mol；
C为气体摩尔浓度，mol/m3；
l为气体传输方向的距离，m；
对于真实气体而言，其状态方程为
pVmol＝ZRT式(3)
式中：p为气体压力，Pa；
Vmol为气体摩尔体积，m3/mol；
T为绝对温度，K；
Z为真实气体压缩因子，量纲1，使用Hall-Yarborough法计算；
由此可知
则
1.2过渡流
当克努森数为0.1～10之间，气体分子在孔道内的传输方式为过渡流。过渡流域内气体传输能力可由努森扩散与滑脱流动综合表征，即
Jtr＝θJsl+(1)Jkn式(6)
式中：Jtr为过渡流摩尔通量，mol/(m2·s)；
Jsl为滑脱流动摩尔通量，mol/(m2·s)；
θ为过渡流中存在滑脱现象的分子数占总分子数的比例系数，量纲1；Kn0.、n为拟合常数，Kn0.＝4.5，n＝5；图1所示为θ与Kn的关系曲线。当Kn低于0.1，θ趋近于1，表明此Kn范围内，气体传输方式中不存在努森扩散；当Kn大于15，θ趋近于0，表明此Kn范围内，气体传输方式中不存在滑脱流、达西流；
1.3滑脱流
当克努森数为0.001～0.1之间，气体分子在孔道内的传输方式为滑脱流动。滑脱流动中气体分子与孔道壁面间的碰撞不可忽略，气体分子在孔道壁面运动速度不为零，需使用滑脱边界条件对Hagen-Poiseuille方程进行修正；
Jsl＝δκJda式(7)
式中：Jsl为滑脱流摩尔通量，mol/(m2·s)；
δκ为滑脱修正系数，量纲1；
式中：μ为气体粘度，Pa·s；
r为孔道半径，m；
κ为切向动量适应系数，κ＝0.8，量纲1；
由于气体分子平均自由程λ亦可表示为则式(8)可以改写为
气体粘度μ计算使用Lee-Gonzden-Enkin半经验法；
Jda为达西流摩尔通量，mol/(m2·s)；
1.4达西流
当克努森数小于0.001，气体分子在孔道内的传输方式为达西流动。气体分子做达西流的达西流摩尔通量表示为
式中：Jda为达西流摩尔通量，mol/(m2·s)；
1.5微孔道表观渗透率模型
孔道内气体传输方式由温度、压力、气体类型、孔道直径及孔隙压力决定，气体传输能力可由相应气体传输方式的摩尔通量表征。不同气体传输方式的摩尔通量表达式中已考虑了孔道内气体分子之间、气体分子与壁面的两种碰撞类型，若采用达西流摩尔通量与努森扩散摩尔通量的相应占比求和的方式表征微孔道气体传输摩尔通量则会造成误差。因此采用分段函数的形式表征一定温度、压力、气体类型及孔隙直径的孔道内气体传输能力，此分段函数的分段界限为克努森数，即
式中：Jd为具有特定孔径的微孔道内气体传输摩尔通量，mol/(m2·s)；
对于多孔介质而言，假设其迂曲度是均匀的，则多孔介质的微孔道内气体传输摩尔通量Jdpm可表示为
式中：Jdpm为多孔介质的微孔道内气体传输摩尔通量，mol/(m2·s)；
τ为迂曲度，量纲1，迂曲度τ与孔隙度满足/>由此可得，多孔介质的微孔道内表观渗透率表达式为
式中：Kapp,表示示多孔介质的微孔道内气体传输表观渗透率，m2；
ρ为气体密度，Kg/m3；
结合式(5)、式(6)、式(7)、式(9)、式(10)、式(11)、式(12)、式(13)可得多孔介质的微孔道内表观渗透率为
式中：
模型表征
利用建立的微孔道表观渗透率模型，结合相关模拟参数(表3)，计算了纳米级孔道(1、5、10、20、50、100、250、500nm)内甲烷传输表观渗透率，计算结果如图2所示；
表3微孔道表观渗透率模拟参数列表
图2表示不同孔隙直径微孔道内气体传输表观渗透率随孔隙压力的变化。孔隙直径越小，微孔道表观渗透率越低，表观渗透率变化幅度越大。例如d＝1nm，Kapp,为(8.55×10-7-6.96×10-5)mD；d＝20nm，Kapp,为(1.55×10-4-7.34×10-4)mD。同时高孔隙压力(10-30MPa)阶段，微孔道表观渗透率随孔隙压力的变化幅度较小，低孔隙压力(0.5-5MPa)阶段，表观渗透率变化幅度大；
本专利建立的微孔道内气体传输表观渗透率模型，充分考虑了特定温压条件、孔隙直径及气体特性等因素下的气体传输方式的差异。图3表示微孔道内气体传输表观渗透率随克努森数的变化，微孔道内气体传输存在着传输形式的变化。例如d＝1nm孔道内，存在扩散流与过渡流的转变；d＝5、10、20、50、100nm孔道内存在过渡流与滑脱流的转变；d＝250、500nm孔道内存在滑脱流与达西流的转变。
实施例2—考虑有效应力效应的页岩表观渗透率模型
上述(一)中以上所建立的微孔道表观渗透率模型是针对具有特定孔隙直径的单个微孔道而言。对于多孔介质而言，假设气体传输通道是由若干具有一定迂曲度及孔径分布的微孔道组成，微孔道之间由不具气体通过能力的物质组成，即毛管束模型。具体假设多孔介质中微孔道的孔隙直径分别为d1，d2，……，di，相应孔隙直径对应的微孔道数目分别为N1，N2，……，Ni。以下讨论内容针对具单位质量的多孔介质而言；具有特定孔隙直径的微孔道表观渗透率表示为Kapp,di，则该微孔道气体流速可按照达西定律表示为
式中为具有特定孔隙直径di的微孔道内气体流速，m/s；
该微孔道气体流量可表示为
式中：为具有特定孔隙直径di的微孔道气体流量，m3/s；
为具有特定孔隙直径di的微孔道截面积，m2；
则具有相同孔隙直径di的微孔道集合的总气体流量为
式中：为具相同孔隙直径di的微孔道集合的总气体流量，m3/s；
Ni为具相同孔隙直径di的微孔道的数量，整数；
多孔介质气体流量为所有微孔道气体流量的综合，可表示为
式中：qt为多孔介质气体流量，m3/s；
对于多孔介质而言，其表观渗透率为Kapp，则多孔介质气体流量qt同时可以表示为
式中：Kapp为多孔介质表观渗透率，m2；
S为多孔介质截面积，m2；
联立式(15)，式(16)，式(17)，式(18)，式(19)，可得
式中：表示具有孔隙直径di的微孔道的总截面积，m2；
表示具有孔隙直径di的微孔道的总截面积与多孔介质截面积比，小数；
定义f(di)为多孔介质内微孔道表观渗透率分配系数，量纲1，令则多孔介质表观渗透率为
由于假设多孔介质的孔道为平行毛细管束结构，因此孔道的截面积占比亦可表示为孔道的体积占比，则f(di)可表示为
式中：为多孔介质孔隙度，％；
为具孔隙直径di的微孔道体积，cm3；
表示多孔介质总孔体积，cm3；
式(22)中涉及的孔隙体积、孔隙直径数据，可通过气体吸附法、压汞法等实验获取，且测点数量越多，孔隙结构刻画越精细，相应的表观渗透率预测结果的准确度越高。同时由于在微孔道表观渗透率中已经考虑了多孔介质孔道的迂曲度等(见式12)，因此此处不再考虑；
对于实际储层而言，产气过程中，上覆地层压力保持恒定，随着孔隙压力的降低，储层所承受的有效应力随之增加，造成孔道压缩，储层渗透率性下降。考虑有效应力条件下的多孔介质表观渗透率可表示为：
式中，σeff＝pc-αp，σeff为有效应力，MPa；
α为Biot常数，量纲1；
pc为上覆地层压力，MPa；
p1为多孔介质孔隙闭合时所承受的有效应力，MPa；当有效应力σeff达到p1，此时Kapp,eff＝0；
β为与孔道表面粗糙度有关的常数，小数；
K0为有效应力为0时的表观渗透率，m2；
由于式(21)所表征的多孔介质表观渗透率未考虑有效应力的影响，因此可视K0＝Kapp。则考虑有效应力影响的多孔介质表观渗透率可表示为
模型验证
选用四川盆地寒武统牛蹄塘组的2块页岩岩心的孔隙结构及渗透率应力敏感数据开展模型验证。岩性样品样品编号分别为MK64-1、MK64-2。岩心基础数据如表4所示；
表4页岩岩心基础数据
岩样编号 层位 长度/cm 直径/cm 孔隙度/％ TOC/％
MK64-1 牛蹄塘组 2.323 2.519 1.797 8.336
MK64-2 牛蹄塘组 3.999 2.534 1.349 3.732
采用PoroPDP-200型覆压孔渗仪，测定了不同有效应力下的页岩岩心渗透率，实验气体为氮气，实验采用定孔隙压力(7MPa)，改变围压的方式实现有效应力(2、5、10、15、25MPa)的变化。由于实验过程中孔隙压力是恒定的，不存在气体的解吸附，因此在此过程中不考虑基质收缩效应的影响；
各有效应力下页岩岩样的渗透率如表5所示；
表5不同有效应力下页岩岩心的渗透率
采用低温液氮吸附实验测定页岩岩心MK64-1、MK64-2的孔隙结构分布，测量结果如图4所示；
依据多孔介质内微孔道表观渗透率分配系数f(di)的定义，页岩样品MK64-1、MK64-2的f(di)值分别采用孔隙结构测定实验获取的不同孔径下的孔体积分布数据计算；
模型计算其余相关参数见表6所示；
表6模型计算相关参数
方程拟合结果如图5所示。不同有效应力下表观渗透率模型预测值偏低于实测值，模型预测表观渗透率应力敏感特征与实测值相符，均呈现随有效应力增加，表观渗透率下降的规律，经计算模型预测值相对于实验值的平均相对误差为12.60％。预测结果表明，本专利建立的基于孔隙结构特征，结合微孔道表观渗透率分段函数表征的页岩表观渗透率预测模型能够有效的预测页岩表观渗透率的动态变化。
实施例3-考虑有效应力及基质收缩作用影响的页岩表观渗透率模型上述(二)中以上所建立的考虑有效应力的页岩表观渗透率模型是针对孔隙压力恒定，有效应力变化的气体流动过程而言。在页岩储层降压产气过程中，随着孔隙压力的降低，孔隙表面的吸附气将解吸为游离气，页岩基质发生收缩，孔道尺寸变大。基质收缩效应将提升孔道的气体传输能力。采用基质收缩渗透率修正系数作为基质收缩对多孔介质表观渗透率的修正系数。基质收缩渗透率修正系数为
式中：ωs为页岩有机质收缩效应渗透率分系数，量纲1；
pint为初始孔隙压力，MPa；
εL为Langmuir应变，量纲1；
pL为Langmuir压力，MPa；
因此考虑有效应力及基质收缩影响的多孔介质表观渗透率可表示为
在实际储层的衰竭式开采中，储层上覆地层压力保持恒定，孔隙压力逐渐下降。本节利用所建立的多孔介质表观渗透率模型，结合四川盆地牛蹄塘组页岩样品MK64-1孔隙结构数据(如图4)，分别模拟了埋深1500m、2000m、3000m，对应围压Pc＝35MPa、45MPa、68MPa，对应初始孔隙压力15MPa、20MPa、30MPa条件下，页岩气传输过程中的页岩渗透率动态变化。
其中页岩气的气体组成及相关参数如表1所示。其余模型参数如表6所示。模拟结果如图6所示。图6所示为模拟不同围压及孔隙压力条件下，页岩表观渗透率的动态变化。
模拟结果显示，随着孔隙压力降低，表观渗透率下降。模拟的孔隙压力范围内，围压35MPa时，表观渗透下降幅度为71.33％，围压45MPa时，下降幅度为93.45％，围压68MPa时，当孔隙压力下降至13.5MPa，表观渗透率接近于0，此时有效应力接近于实验获得的页岩孔隙闭合时的有效应力。可见，围压与表观渗透率的下降幅度呈正相关关系。
由于低孔隙压力阶段气体传输方式向努森扩散方向转变，以及吸附气解吸基质收缩作用造成孔道扩张，二者均能够起到提高气体传输能力的作用，而图6显示，随着孔隙压力的降低，页岩表观渗透率持续降低，这表明相对于气体传输方式转变及基质收缩作用而言，有效应力作用在页岩表观渗透率动态变化过程中起主导作用。

Claims (3)

1.一种基于孔隙结构特征的页岩表观渗透率预测方法，其特征在于：该步骤如下：(1)建立微孔道表观渗透率模型如下：
式中：Kapp,d为多孔介质内具有孔隙直径的微孔道表观渗透率，m2；
d为孔道直径，m；
μ为气体粘度，Pa·s；
M为气体摩尔质量，kg/mol；
τ为迂曲度，量纲1；
Z为真实气体压缩因子，量纲1；
ρ为气体密度，Kg/m3；
R为通用气体常数，8.314J/(mol·K)；
T为气体温度，℃；
M为气体摩尔质量，kg/mol；
θ为过渡流中存在滑脱现象的分子数占总分子数的比例系数，量纲1；
δκ为滑脱修正系数，量纲1；
p为气体压力，Pa；
Kn为克努森数，量纲1；
κ为切向动量适应系数，量纲1；
Kn0.5、n为拟合常数，Kn0.5＝4.5，n＝5；
kb为玻尔兹曼常数，1.38×10-23J·K-1；
σ为气体分子平均碰撞直径，m；
(2)建立多孔介质表观渗透率模型如下：
式中：Kapp为多孔介质表观渗透率，m2；
f(di)为多孔介质内微孔道表观渗透率分配系数，量纲1；
Kapp,di为具有特定孔隙直径di的微孔道表观渗透率，m2；
式中：为多孔介质孔隙度，％；
Ni为具有相同孔隙直径di的微孔道的数量，整数；
为具孔隙直径di的微孔道体积，cm3；
表示多孔介质总孔体积，cm3。
2.根据权利要求1所述基于孔隙结构特征的页岩表观渗透率预测方法，其特征在于：还包括在有效应力影响下的多孔介质表观渗透率模型为
式中：Kapp,eff为考虑有效应力条件下的多孔介质表观渗透率，m2；
σeff＝pc-αp，σeff为有效应力，MPa；
pc为上覆地层压力，MPa；
α为Biot系数，量纲1；
p1为多孔介质孔隙闭合时所承受的有效应力，MPa；
β为与孔道表面粗糙度有关的常数，小数。
3.根据权利要求2所述基于孔隙结构特征的页岩表观渗透率预测方法，其特征在于：还包括在有效应力及基质收缩影响的多孔介质表观渗透率模型为
式中：ωs为页岩有机质收缩效应渗透率修正系数，量纲1；
为考虑有效应力及基质收缩影响的多孔介质表观渗透率；
式中：pint为初始孔隙压力，MPa；
L为Langmuir应变，量纲1；
pL为Langmuir压力，MPa。
CN202110389015.1A 2021-04-12 2021-04-12 一种基于孔隙结构特征的页岩表观渗透率预测方法 Active CN113128040B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110389015.1A CN113128040B (zh) 2021-04-12 2021-04-12 一种基于孔隙结构特征的页岩表观渗透率预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110389015.1A CN113128040B (zh) 2021-04-12 2021-04-12 一种基于孔隙结构特征的页岩表观渗透率预测方法

Publications (2)

Publication Number Publication Date
CN113128040A CN113128040A (zh) 2021-07-16
CN113128040B true CN113128040B (zh) 2024-02-20

Family

ID=76776295

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110389015.1A Active CN113128040B (zh) 2021-04-12 2021-04-12 一种基于孔隙结构特征的页岩表观渗透率预测方法

Country Status (1)

Country Link
CN (1) CN113128040B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114136862B (zh) * 2021-11-29 2024-04-09 西南石油大学 一种双润湿性页岩的液体表观渗透率计算方法
CN118067582A (zh) * 2022-11-22 2024-05-24 中国石油天然气股份有限公司 一种岩石原位转化动态储层物性变化评价方法及系统
CN116124667B (zh) * 2022-11-23 2024-02-09 西南石油大学 一种页岩气储层渗透率表征方法

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107461192A (zh) * 2017-06-01 2017-12-12 西南石油大学 储层条件下页岩动态表观渗透率的计算方法
CN108343433A (zh) * 2018-02-28 2018-07-31 西南石油大学 用于计算页岩微裂缝缝宽变化下气体质量传输参数的方法
CN108959789A (zh) * 2018-07-12 2018-12-07 西南石油大学 一种考虑应力敏感效应的页岩表观渗透率计算方法
CN109902918A (zh) * 2019-01-16 2019-06-18 西南石油大学 一种页岩压裂自支撑裂缝缝网渗透率计算方法
CN110210460A (zh) * 2019-06-26 2019-09-06 中国石油大学(华东) 一种考虑多重因素影响的页岩气表观渗透率计算方法
CN110567858A (zh) * 2019-10-17 2019-12-13 西南石油大学 一种基于分形理论的预测页岩纳米孔隙渗透率的方法
CN111208052A (zh) * 2020-02-28 2020-05-29 西安石油大学 基于改进Kozeny-Carman模型的页岩储层渗透率预测方法
CN111208051A (zh) * 2020-02-28 2020-05-29 西安石油大学 复杂机理协同影响下页岩储层渗透率动态预测方法
CN111340298A (zh) * 2020-02-28 2020-06-26 西安石油大学 页岩储层诱导非均质性动态演化预测方法

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107461192A (zh) * 2017-06-01 2017-12-12 西南石油大学 储层条件下页岩动态表观渗透率的计算方法
CN108343433A (zh) * 2018-02-28 2018-07-31 西南石油大学 用于计算页岩微裂缝缝宽变化下气体质量传输参数的方法
CN108959789A (zh) * 2018-07-12 2018-12-07 西南石油大学 一种考虑应力敏感效应的页岩表观渗透率计算方法
CN109902918A (zh) * 2019-01-16 2019-06-18 西南石油大学 一种页岩压裂自支撑裂缝缝网渗透率计算方法
CN110210460A (zh) * 2019-06-26 2019-09-06 中国石油大学(华东) 一种考虑多重因素影响的页岩气表观渗透率计算方法
CN110567858A (zh) * 2019-10-17 2019-12-13 西南石油大学 一种基于分形理论的预测页岩纳米孔隙渗透率的方法
CN111208052A (zh) * 2020-02-28 2020-05-29 西安石油大学 基于改进Kozeny-Carman模型的页岩储层渗透率预测方法
CN111208051A (zh) * 2020-02-28 2020-05-29 西安石油大学 复杂机理协同影响下页岩储层渗透率动态预测方法
CN111340298A (zh) * 2020-02-28 2020-06-26 西安石油大学 页岩储层诱导非均质性动态演化预测方法

Also Published As

Publication number Publication date
CN113128040A (zh) 2021-07-16

Similar Documents

Publication Publication Date Title
CN113128040B (zh) 一种基于孔隙结构特征的页岩表观渗透率预测方法
CN107461192B (zh) 储层条件下页岩动态表观渗透率的计算方法
CN108959789B (zh) 一种考虑应力敏感效应的页岩表观渗透率计算方法
Cai et al. Shale gas transport model in 3D fractal porous media with variable pore sizes
Singh et al. Nonempirical apparent permeability of shale
Barsotti et al. A review on capillary condensation in nanoporous media: Implications for hydrocarbon recovery from tight reservoirs
CN109100278B (zh) 一种考虑页岩孔径分布特征的表观渗透率计算方法
Wu et al. A model for gas transport in micro fractures of shale and tight gas reservoirs
Dumore et al. Drainage capillary-pressure functions and the influence of connate water
Kang et al. Carbon dioxide storage capacity of organic-rich shales
Yuan et al. An analytical model of apparent gas permeability for tight porous media
CN108229048B (zh) 一种考虑页岩基质差异的多场耦合渗透率计算方法
CN108710723B (zh) 一种页岩储层多孔介质表观渗透率的计算方法
CN110598167B (zh) 低渗透油藏油水相对渗透率实验数据的处理方法
Wua et al. Adsorbed gas surface diffusion and bulk gas transport in nanopores of shale reservoirs with real gas effect-adsorption-mechanical coupling
Xu et al. Real gas transport in tapered noncircular nanopores of shale rocks
Sheng et al. Pore-scale modeling and analysis of surface diffusion effects on shale-gas flow in Kerogen pores
CN110567858A (zh) 一种基于分形理论的预测页岩纳米孔隙渗透率的方法
CN104200105A (zh) 确定致密砂岩气充注物性下限的方法
Erfani Gahrooei et al. Wettability alteration of reservoir rocks to gas wetting condition: A comparative study
CN109085104B (zh) 一种致密气藏的气体相对渗透率的确定方法
CN112129682B (zh) 一种基于边界层动态变化的低渗透油藏驱替介质优选方法
Kuila et al. Application of Knudsen flow in modeling gas-flow in shale reservoirs
CN105675441A (zh) 一种径向流条件下不同水合物饱和度时的气水相对渗透率测定方法
San Manley et al. An integrated method for measuring gas permeability and diffusivity of porous solids

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