岩石力学与工程学报

Parameter	$p L1 m$ /MPa	$p L2 m$ /MPa	$ε L1 m$ /%	$ε L2 m$ /%	$ε Lf$ /%	$p c$ /MPa
Value range	1～2.3^[28]	4～7^[28]	0.08～0.12^[27]	0.12～0.2	0.04～0.08	7.2～8.0
Value	2	6	0.08	0.12	0.04	7.4

Parameter

p L1 m

/MPa

p L2 m

/MPa

ε L1 m

ε L2 m

ε Lf

p c

/MPa

Value range

1～2.3^[28]

4～7^[28]

0.08～0.12^[27]

0.12～0.2

0.04～0.08

7.2～8.0

Value

0.08

0.12

0.04

7.4

Parameter	$p L1 m$ /MPa	$p L2 m$ /MPa	$ε L1 m$ /%	$ε L2 m$ /%	$ε Lf$ /%	$p c$ /MPa
Value range	1～2.3^[28]	4～7^[28]	0.08～0.12^[27]	0.12～0.2	0.04～0.08	7.2～8.0
Value	2	6	0.08	0.12	0.04	7.4

Parameter

p L1 m

/MPa

p L2 m

/MPa

ε L1 m

ε L2 m

ε Lf

p c

/MPa

Value range

1～2.3^[28]

4～7^[28]

0.08～0.12^[27]

0.12～0.2

0.04～0.08

7.2～8.0

Value

0.08

0.12

0.04

7.4

Parameter	Degree of deterioration/%	Literature sources
Young's modulus	－24.60～26.79 －29.28～33.15 －20.23～37.39	Q. Lyu et al.^[18] H. Yin et al.^[31] Q. Lyu et al.^[33]
Poisson's ratio	29.41～58.82	Q. Lyu et al.^[33]

Parameter

Degree of deterioration/%

Literature sources

Young's modulus

－24.60～26.79
－29.28～33.15
－20.23～37.39

Q. Lyu et al.^[18]
H. Yin et al.^[31]
Q. Lyu et al.^[33]

Poisson's ratio

29.41～58.82

Q. Lyu et al.^[33]

Parameter	Degree of deterioration/%	Literature sources
Young's modulus	－24.60～26.79 －29.28～33.15 －20.23～37.39	Q. Lyu et al.^[18] H. Yin et al.^[31] Q. Lyu et al.^[33]
Poisson's ratio	29.41～58.82	Q. Lyu et al.^[33]

Parameter

Degree of deterioration/%

Literature sources

Young's modulus

－24.60～26.79
－29.28～33.15
－20.23～37.39

Q. Lyu et al.^[18]
H. Yin et al.^[31]
Q. Lyu et al.^[33]

Poisson's ratio

29.41～58.82

Q. Lyu et al.^[33]

Initial porosity of matrix $ϕ mo$	Initial porosity of fracture $ϕ fo$	Maximum adsorption capacity of matrix part $n max m$ /(m³·kg)	Maximum adsorption capacity of fracture part $n max f$ /(m³·kg)	Biot coefficient of matrix $α$	Biot coefficient of fracture $β$	Initial permeability of matrix part $k m0$ /m²	Initial permeability of fracture part $k f0$ /m²	Sensitivity exponent M	Injected target pressure $p in$ /MPa
0.04	0.04	0.12	0.05	0.65	0.8	5×10^－19	1×10^－18	3	14

Initial porosity of matrix

ϕ mo

Initial porosity of fracture

ϕ fo

Maximum adsorption capacity of matrix part

n max m

/(m³·kg)

Maximum adsorption capacity of fracture part

n max f

/(m³·kg)

Biot coefficient of matrix

α

Biot coefficient of fracture

β

Initial permeability of matrix part

k m0

/m²

Initial permeability of fracture part

k f0

/m²

Sensitivity exponent M

Injected target pressure

p in

/MPa

0.04

0.12

0.05

0.65

0.8

5×10^－19

1×10^－18

Initial porosity of matrix $ϕ mo$	Initial porosity of fracture $ϕ fo$	Maximum adsorption capacity of matrix part $n max m$ /(m³·kg)	Maximum adsorption capacity of fracture part $n max f$ /(m³·kg)	Biot coefficient of matrix $α$	Biot coefficient of fracture $β$	Initial permeability of matrix part $k m0$ /m²	Initial permeability of fracture part $k f0$ /m²	Sensitivity exponent M	Injected target pressure $p in$ /MPa
0.04	0.04	0.12	0.05	0.65	0.8	5×10^－19	1×10^－18	3	14

Initial porosity of matrix

ϕ mo

Initial porosity of fracture

ϕ fo

Maximum adsorption capacity of matrix part

n max m

/(m³·kg)

Maximum adsorption capacity of fracture part

n max f

/(m³·kg)

Biot coefficient of matrix

α

Biot coefficient of fracture

β

Initial permeability of matrix part

k m0

/m²

Initial permeability of fracture part

k f0

/m²

Sensitivity exponent M

Injected target pressure

p in

/MPa

0.04

0.12

0.05

0.65

0.8

5×10^－19

1×10^－18

Basic expression of model	Parameter interpretation	Literature sources
$k k 0 = exp (− 3 C f Δ σ − α biot Δ p)$	$Δ σ$ is the stress change; $C f$ is the fracture compression coefficient,1/MPa	J. Q. Shi and S. Durucan^[49]
$k k 0 = (1 + α biot ϕ 0 Δ ε e) 3 f (K n)$	$Δ ε e$ is the effective strain change, $f (K n)$ is the slip correction coefficient related to the Knudsen number	M. Y. Wei et al.^[50]
$k k 0 = exp [3 (1 K − 1 K f) (Δ σ − Δ p) + Δ ε lfs] (l l 0) 2 f (K n) f (K n0)$	L is the fracture length; $ε lfs$ is the fracture strain determined by the volume fraction of adsorbed gas invading the shale, gas diffusion coefficient, and diffusion distance, defined as a piecewise function	J. Zeng et al.^[51]
$k k 0 = {exp (− 3 C f Δ σ e) (Elastic stage) exp (− 3 D λ D C f Δ σ e) (Plastic stage)$	$Δ σ e$ is the effective stress change; $λ D$ is the influence coefficient of rock damage on permeability; $D$ is the damage coefficient, defined as a piecewise function based on tensile and shear strains	C. Y. Wang et al.^[52]
$k k 0 = (ϕ ϕ 0) 3 exp (α k D)$	$α k$ is the fitted value of the influence coefficient of rock mass damage on permeability	W. C. Zhu et al.^[47]
$k k 0 = (ϕ ϕ 0) M + 1 exp (η z 0 − η z)$	M is the stress sensitivity constant obtained from experimental fitting,z is the initial coordination number of pores and fractures,and $η$ is a constant related to the fractal dimension and M.	D. Davudov and R. G. Moghanloo^[41]

Basic expression of model

Parameter interpretation

Literature sources

k k 0 = exp (− 3 C f Δ σ − α biot Δ p)

Δ σ

is the stress change;

C f

is the fracture compression coefficient,1/MPa

J. Q. Shi and S. Durucan^[49]

k k 0 = (1 + α biot ϕ 0 Δ ε e) 3 f (K n)

Δ ε e

is the effective strain change,

f (K n)

is the slip correction coefficient related to the Knudsen number

M. Y. Wei et al.^[50]

k k 0 = exp [3 (1 K − 1 K f) (Δ σ − Δ p) + Δ ε lfs] (l l 0) 2 f (K n) f (K n0)

L is the fracture length;

ε lfs

is the fracture strain determined by the volume fraction of adsorbed gas invading the shale, gas diffusion coefficient, and diffusion distance, defined as a piecewise function

J. Zeng et al.^[51]

k k 0 = {exp (− 3 C f Δ σ e) (Elastic stage) exp (− 3 D λ D C f Δ σ e) (Plastic stage)

Δ σ e

is the effective stress change;

λ D

is the influence coefficient of rock damage on permeability;

D

is the damage coefficient, defined as a piecewise function based on tensile and shear strains

C. Y. Wang et al.^[52]

k k 0 = (ϕ ϕ 0) 3 exp (α k D)

α k

is the fitted value of the influence coefficient of rock mass damage on permeability

W. C. Zhu et al.^[47]

k k 0 = (ϕ ϕ 0) M + 1 exp (η z 0 − η z)

M is the stress sensitivity constant obtained from experimental fitting,z is the initial coordination number of pores and fractures,and

η

is a constant related to the fractal dimension and M.

D. Davudov and R. G. Moghanloo^[41]

Basic expression of model	Parameter interpretation	Literature sources
$k k 0 = exp (− 3 C f Δ σ − α biot Δ p)$	$Δ σ$ is the stress change; $C f$ is the fracture compression coefficient,1/MPa	J. Q. Shi and S. Durucan^[49]
$k k 0 = (1 + α biot ϕ 0 Δ ε e) 3 f (K n)$	$Δ ε e$ is the effective strain change, $f (K n)$ is the slip correction coefficient related to the Knudsen number	M. Y. Wei et al.^[50]
$k k 0 = exp [3 (1 K − 1 K f) (Δ σ − Δ p) + Δ ε lfs] (l l 0) 2 f (K n) f (K n0)$	L is the fracture length; $ε lfs$ is the fracture strain determined by the volume fraction of adsorbed gas invading the shale, gas diffusion coefficient, and diffusion distance, defined as a piecewise function	J. Zeng et al.^[51]
$k k 0 = {exp (− 3 C f Δ σ e) (Elastic stage) exp (− 3 D λ D C f Δ σ e) (Plastic stage)$	$Δ σ e$ is the effective stress change; $λ D$ is the influence coefficient of rock damage on permeability; $D$ is the damage coefficient, defined as a piecewise function based on tensile and shear strains	C. Y. Wang et al.^[52]
$k k 0 = (ϕ ϕ 0) 3 exp (α k D)$	$α k$ is the fitted value of the influence coefficient of rock mass damage on permeability	W. C. Zhu et al.^[47]
$k k 0 = (ϕ ϕ 0) M + 1 exp (η z 0 − η z)$	M is the stress sensitivity constant obtained from experimental fitting,z is the initial coordination number of pores and fractures,and $η$ is a constant related to the fractal dimension and M.	D. Davudov and R. G. Moghanloo^[41]

Basic expression of model

Parameter interpretation

Literature sources

k k 0 = exp (− 3 C f Δ σ − α biot Δ p)

Δ σ

is the stress change;

C f

is the fracture compression coefficient,1/MPa

J. Q. Shi and S. Durucan^[49]

k k 0 = (1 + α biot ϕ 0 Δ ε e) 3 f (K n)

Δ ε e

is the effective strain change,

f (K n)

is the slip correction coefficient related to the Knudsen number

M. Y. Wei et al.^[50]

k k 0 = exp [3 (1 K − 1 K f) (Δ σ − Δ p) + Δ ε lfs] (l l 0) 2 f (K n) f (K n0)

L is the fracture length;

ε lfs

is the fracture strain determined by the volume fraction of adsorbed gas invading the shale, gas diffusion coefficient, and diffusion distance, defined as a piecewise function

J. Zeng et al.^[51]

k k 0 = {exp (− 3 C f Δ σ e) (Elastic stage) exp (− 3 D λ D C f Δ σ e) (Plastic stage)

Δ σ e

is the effective stress change;

λ D

is the influence coefficient of rock damage on permeability;

D

is the damage coefficient, defined as a piecewise function based on tensile and shear strains

C. Y. Wang et al.^[52]

k k 0 = (ϕ ϕ 0) 3 exp (α k D)

α k

is the fitted value of the influence coefficient of rock mass damage on permeability

W. C. Zhu et al.^[47]

k k 0 = (ϕ ϕ 0) M + 1 exp (η z 0 − η z)

M is the stress sensitivity constant obtained from experimental fitting,z is the initial coordination number of pores and fractures,and

η

is a constant related to the fractal dimension and M.

D. Davudov and R. G. Moghanloo^[41]

页岩注入CO₂渗透率多阶段演化特性的模型解析

PDF下载

李文睿 ¹^,²^,³ , 王登科 ¹^,²^,³^,^* , 魏建平 ¹^,²^,³ , 余金昊 ¹^,² , 赵文霖 ¹^,² , 温志辉 ¹^,²^,³

岩石力学与工程学报 | 理论与试验研究 2026,45(2): 432-448

收起

岩石力学与工程学报 | 理论与试验研究 2026, 45(2): 432-448

页岩注入CO₂渗透率多阶段演化特性的模型解析

全屏

李文睿¹^,²^,³, 王登科¹^,²^,³^,^*, 魏建平¹^,²^,³, 余金昊¹^,², 赵文霖¹^,², 温志辉¹^,²^,³

作者信息

^1.河南理工大学安全科学工程学院，河南焦作　454000

^2.河南理工大学河南省瓦斯地质与瓦斯治理重点实验室–省部共建国家重点实验室培育基地，河南焦作　454000

^3.河南理工大学煤矿瓦斯灾害防治应急管理部重点实验室，河南焦作　454000

LI Wenrui (1994–), assistant professor, is engaged in research on unconventional natural gas development and carbon dioxide geological storage. E-mail：li3wenrui123@163.com

李文睿（1994–），现任讲师，主要从事非常规天然气开发及二氧化碳地质封存方面的研究工作。E-mail：li3wenrui123@163.com

通讯作者:

^* 王登科（1980–），现任教授，主要从事安全科学与工程方面的研究工作。E-mail：wdk@hpu.edu.cn

Model interpretation of the multi-stage evolution characteristics of shale permeability during the injection of CO₂

Wenrui LI¹^,²^,³, Dengke WANG¹^,²^,³^,^*, Jianping WEI¹^,²^,³, Jinhao YU¹^,², Wenlin ZHAO¹^,², Zhihui WEN¹^,²^,³

Affiliations

^1.School of Safety Science and Engineering, Henan Polytechnic University, Jiaozuo, Henan 454000, China

^2.State Key Laboratory Cultivation Base for Gas Geology and Gas Control, Henan Polytechnic University, Jiaozuo, Henan 454000, China

^3.Key Laboratory of Coal Mine Disasters Prevention, Henan Polytechnic University, Jiaozuo, Henan 454000, China

出版时间: 2026-02-01 doi: 10.3724/1000-6915.jrme.2025.0505

文章导航

摘要

收起

精准预测CO₂注入页岩储层的渗透率动态演化规律对于碳封存和增强页岩气开采具有重要意义。针对传统渗透率模型难以准确描述CO₂由低压气态到超临界态全压注入页岩过程中的渗透率全过程演化问题，基于基质–裂隙双弹性介质系统，构建由组分渗透率权重决定的页岩渗透率演化模型；通过引入基质力学劣化、二次吸附，以及应变滞后效应等关键因素，建立起多效应耦合控制的渗透率演化方程。基于双弹性介质体重叠结构开展并行交叉耦合数值求解，实现对CO₂全压注入过程中页岩渗透率非线性演化过程的准确描述。进一步，影响效应解耦研究结果表明：CO₂对基质材料力学参数的劣化界定渗透率波动范围的边界阈值，机械应变与吸附应变之间的非同步响应显著放大演化阶段差异，形成具有显著辨识度的阶段转换过程；应变滞后效应则显著延长演化时间跨度；气体吸附作用与力学响应共同调控各演化阶段的转变节点，尤其是二次吸附引发的额外膨胀应变明显增强了演化过程中的阶段划分程度。此外，本研究还深入分析流–固耦合渗透率建模基本框架和不同数值模拟方法的特性，不仅可以加深对CO₂注入过程中页岩渗透率演化过程的理解，还为流体地质封存渗透率理论建模及数值模拟方法提供借鉴意义。

关键词

岩石力学 / 碳封存 / 页岩 / 渗透率演化 / 超临界态 / 二次吸附

Abstract

收起

Accurately predicting the dynamic evolution of permeability during CO₂ injection into shale reservoirs is crucial for carbon sequestration and enhanced shale gas recovery. However, traditional permeability models often fail to comprehensively describe the full-range evolution of permeability throughout the entire CO₂ injection process in shale—from the low-pressure gaseous state to the supercritical state. To address this limitation, this study develops a shale permeability evolution model based on a dual-elastic system comprising both the matrix and fractures, determined by component permeability weighting. By incorporating key factors such as mechanical degradation of the matrix, secondary adsorption, and strain hysteresis effects, we establish a governing equation for permeability evolution under multi-effect coupling. Utilizing an overlapping dual-elastic medium structure, we perform parallel cross-coupling numerical solutions, achieving an accurate representation of the nonlinear permeability evolution during full-pressure CO₂ injection. Furthermore, a decoupled analysis of influencing effects reveals that the degradation of mechanical parameters of the matrix material due to CO₂ defines the boundary thresholds for permeability fluctuation ranges. The asynchronous response between mechanical strain and adsorption strain significantly amplifies differences across evolutionary stages, leading to clearly distinguishable phase transitions. Additionally, the strain hysteresis effect prolongs the duration of evolution. Gas adsorption and mechanical responses jointly regulate the transition points between evolutionary stages, with the secondary adsorption-induced swelling strain particularly enhancing phase differentiation throughout the evolution process. This study also provides an in-depth analysis of the fundamental framework of fluid-solid coupled permeability modeling and explores the characteristics of different numerical simulation methods. The findings not only deepen the understanding of shale permeability evolution during CO₂ injection but also offer valuable insights for theoretical modeling and numerical simulation of permeability in geological fluid sequestration.

Key words

rock mechanics / carbon sequestration / shale / permeability evolution / supercritical state / secondary adsorption

引用本文

李文睿, 王登科, 魏建平, 余金昊, 赵文霖, 温志辉. 页岩注入CO₂渗透率多阶段演化特性的模型解析. 岩石力学与工程学报, 2026 , 45 (2) : 432 -448 . DOI: 10.3724/1000-6915.jrme.2025.0505

Wenrui LI, Dengke WANG, Jianping WEI, Jinhao YU, Wenlin ZHAO, Zhihui WEN. Model interpretation of the multi-stage evolution characteristics of shale permeability during the injection of CO₂[J]. Chinese Journal of Rock Mechanics and Engineering, 2026 , 45 (2) : 432 -448 . DOI: 10.3724/1000-6915.jrme.2025.0505

正文

收起

1　引言

收起

页岩渗透率演化的精准预测是评估储层封存性能和强化页岩气开采的关键科学问题^[1]。向深部储层中注入气体的过程，可将其抽象为恒定围压(constant confining pressure，CCP)边界条件模型，众多室内物理模拟试验围绕于此展开试验研究。稳态法测试渗透率结果表明，氮气、甲烷和CO₂等吸附性气体注入，页岩渗透率普遍呈现“U型”演化特征，控制其演化的核心机制可归纳为：(1)孔隙压力驱动的有效应力重分布；(2)气体吸附引发的基质–裂隙变形竞争；(3)化学腐蚀导致的力学性质劣化；(4)气体多机制传输效应^[2-3]。理想状态下，页岩渗透率演化范围上限由有效应力控制的多孔弹性变形决定，下限则由滑脱效应决定^[4-6]，有效应力和滑脱效应的共同作用划定了实际渗透率演化区间^[7]。随着研究不断深入，发现存在相态转化的CO₂注入，页岩渗透率呈现出多阶段演化特征，非简单的“U型”演化规律^[8]。

为深入揭示页岩渗透率演化规律，最初理论模型将页岩渗透率直接视为裂隙渗透率，由于忽略了基质–裂隙相互作用，导致模型计算结果与试验数据不完全匹配^[9]。实际上，基质–裂隙相互作用受整体变形和裂隙开度变化共同控制^[10]。基于这一认知，学者们引入了应变分配模型，用来区分吸附应变对基质膨胀和裂隙开度的独立影响，数学定义了应变分割系数^[11]、内部膨胀因子^[12]和气体非平衡因子^[13]。研究结果表明，引入应变分配概念的双重孔隙介质渗透率模型，能很好地解释常围压条件下，渗透率上、下限包络区间内存在的不同“U型”演化趋势。W. X. Cheng等^[14]还提出了一个多层纳米壳模型，通过模拟气体吸附过程中纳米管壁和基质外壳的竞争变形作用，发现了“U型”渗透率后期回弹不明显的特异性情况。然而，在全压CO₂注入范围(低压气态到高压超临界态)内，CO₂–页岩相互作用易导致矿物溶解、孔隙结构变化和基质力学性质劣化等现象^[15-18]，强扩散性超临界CO₂尤其对微孔结构损伤明显^[19]。此外，气楔效应驱动裂隙网络扩展形成优势渗流通道，诱发岩体局部应力集中与弱面损伤演化，易导致渗透率出现突变特征^[20]。试验研究发现，孔隙比表面积和孔隙率会因矿物溶解和碳酸盐形成有所增加；在超临界CO₂相变增能与孔隙匹配的机制下，CO₂分子吸附空间向微孔演化，进而通过多层吸附与微孔填充等途径，实现了吸附能力的再次强化^[21]。L. Hou和D. Elsworth^[8]将这种现象称之为“二次吸附”，并通过含相态转化的CO₂注入试验，推断二次吸附及其诱发的间歇性膨胀应变正是导致渗透率出现异常演化现象的根本原因。

当前，无论是基于应变分配的渗透率模型^[8，15]、定义局部非均匀应变的模型^[14]，还是强调杨氏模量衰减的双孔介质模型^[22]，虽都与“U型”渗透率试验数据吻合良好，但均无法合理解释全压CO₂注入过程中页岩渗透率呈现出类“W型”的多阶段演化特征。S. W. Zhang等^[23]认为，基质和裂隙都有其独立的渗透率演化特征，仅依靠裂隙渗透率演化来代表页岩渗透率演化，结果可能不准确。L. Fan和S. M. Liu^[24]考虑到外骨架和吸附矿物的体积模量差异，提出了一个概念性的双材料模型，该模型很好地展示了气体吸附诱发的微结构非均匀变形。M. Y. Wei等^[25-26]综合考虑吸附诱发的非均匀应变滞后效应构建了煤渗透率演化模型，该模型不但阐明了控制渗透率“U型”演化的主要机制，还从理论上阐述了注气初期渗透率短暂增大的现象，对于页岩渗透率建模具有一定的借鉴意义。

综上所述，传统渗透率模型无法准确描述页岩注入CO₂过程中的渗透率非线性演化问题。本研究综合考虑基质–裂隙跨尺度流动、力学性质劣化、基质二次吸附膨胀应变，以及应变滞后效应等关键因素；构建了基质–裂隙平行竞争耦合渗透率演化模型，应用重叠连续介质数值模拟方法，揭示了全压CO₂注入过程中的渗透率多阶段演化基本模式，并厘清了各阶段演化的主控机制，从而为页岩储层封存效率优化与渗流稳定性评估提供理论工具。

2　双弹性孔隙介质材料耦合模型

收起

在页岩储层多物理场耦合分析中，本研究构建了裂隙–基质并联的双介质材料模型，其中，软质裂隙系统呈现低吸附性特征，而硬质基质系统具有强吸附性特征。如图1(a)，(b)所示，该模型将页岩实体离散为分布式物理节点网络，每个节点均包含裂隙系统与基质系统的交互作用单元。

如图1(c)所示，采用重叠连续介质方法，将每个物理节点解耦为4个数学节点，分别对应裂隙/基质变形方程、裂隙/基质体内流动方程，建立包含4个相互作用物理场的控制方程组，并通过3组交叉耦合关系实现全耦合求解：基质–裂隙界面处的局部力平衡条件、变形协调条件以及气体质量交换项。图1(d)为页岩岩心轴向剖面图，尺寸为5 cm×10 cm，p_in为注气压力，p_out为出口压力；p_i为注气压力变量，下标i取m和f时分别表示基质系统和裂隙系统；p_a为大气压力；p为孔隙压力；n为边界的单位法向量；u为位移向量；

Ω

指代几何边界，

Ω 1

～

Ω 4

分别为下、左、右和上边界。

数值模拟在CCP条件下进行，CO₂从顶部注入，左右两侧施加恒定外部应力限制变形，底部无流出，计算域内初始孔隙压力设定为0.1 MPa。

3　多物理场数值耦合关系

收起

裂隙和基质系统中体积变形与气体流动的控制方程基于以下假设：(1)页岩由基质和裂隙连续介质材料组成，且均为孔隙介质弹性系统；(2)应变属于微应变；(3)计算域等温且气体完全饱和；(4)不考虑页岩的含水情况。

3.1　力学平衡控制方程

弹性介质的应变–位移满足如下基本关系：

ε i j = (u i ， j + u j ， i) / 2

(1)

式中：

ε i j

为总应变张量，

u i ， j

和

u j ， i

均为位移分量。

固体材料的应力平衡基本方程为

σ i j ， j + f i = 0

(2)

式中：

σ i j ， j

为应力张量，

f i

为体积力。

基于多孔弹性线性本构关系，将体积应变分解为气体吸附应变与局部应变，并考虑基质–裂隙间压力差的相互作用，根据连续介质理论可定义页岩基质与裂缝的变形本构关系^[23]分别为

ε m i j = σ m i j 2 G m − (1 6 G m − 1 9 K m) σ k k δ i j + (α p m − Δ p 3 K m + ε fs 3) δ i j

(3)

ε f i j = σ f i j 2 G f − (1 6 G f − 1 9 K f) σ k k δ i j + (β p f − Δ p 3 K f + ε fs 3) δ i j

(4)

式中：

δ i j

为克罗内克符号；

ε m i j

和

ε f i j

分别为基质与裂隙吸附应变；

Δ p

为基质内孔隙压力

p m

与裂隙内孔隙压力

p f

之差，

Δ p = p f − p m

；

K i

和

G i

分别为弹性体的体积模量和剪切模量，

K i = E i / [3 (1 − 2 ν i)]

、

G i =

E i / [2 (1 + ν i)]

，

E i

和

ν i

分别为杨氏模量和泊松比；

α

和

β

分别为基质和裂隙系统的Biot系数；

σ k k

为弹性体所受的总应力，具有方向性。参数下标i为m和f，分别代表基质系统、裂隙系统。

联立式(1)～(4)，分别得到如下所示基质和裂隙变形的纳维斯托克斯型方程：

G m u i ， k k + 2 G m 1 − 2 ν m u k ， k i − α p m ， i − K m ε ms ， i + Δ p i + F i = 0

(5)

G f u i ， k k + 2 G f 1 − 2 ν f u k ， k i − β p f ， i − K f ε fs ， i − Δ p i + F i = 0

(6)

(1)全压力范围内CO₂吸附及吸附应变

采用结合了吸附密度修正因子的超临界Dubinin- Astakhov(SDA)模型来量化全压范围内CO₂吸附量变化，数学表达式^[21]如下：

n ex = n max exp {− γ [ln (ρ a ρ g)] 2} (1 − ρ g ρ a)

(7)

式中：

n ex

为过剩吸附量(m³/kg)，

n max

为最大吸附量，(m³/kg)，

ρ a

为CO₂真实密度(kg/m³)，

γ

为特征能量和吸附势拟合常数。其中，

n max exp {− γ [ln (ρ a / ρ g)] 2}

为微孔填充吸附模型(SDA)的绝对吸附量n_ab，

1 − ρ g / ρ a

为密度修正项。

将SDA模型分别应用于具有不同吸附特性的基质系统和裂隙系统，得到以下结果：

G m u i ， k k + 2 G m 1 − 2 v m u k ， k i − α p m ， i − K m ε ms ， i + Δ p i + F i = 0

(8)

G f u i ， k k + 2 G f 1 − 2 ν f u k ， k i − β p f ， i − K f ε fs ， i − Δ p i + F i = 0

(9)

在低压气态向超临界态高压等温注入过程中，CO₂密度和动力黏度会发生阶跃性突变(见图2)，本文采用内置Span-Wagner算法的REFPROP软件，联合COMSOL数值，分别计算基质系统和裂隙系统内的密度

ρ gm

和

ρ gf

。

如图2所示，CO₂吸附应变可划分为气态区和超临界区2个阶段。S. Zhi等^[27]试验研究表明，各自区域内的吸附应变均遵循Langmuir型函数：

ε s = ε L p /

(p + p L)

，式中，

ε s

为吸附应变，

ε L

为Langmuir应变常数，p_L为Langmuir压力常数。因此，在数值模拟CO₂全压注入过程中，基质系统和裂隙系统的吸附应变可表示为

ε i s = {ε L 1 i (p i p i + p L 1 i) (p i ≤ p c) ε L 2 i (p i p i + p L 2 i) (p c ＜ p i)

(10)

式中：

p L 1

和

ε L 1 i

分别为气态阶段的Langmuir压力和应变常数；

p L 2

和

ε L 2 i

分别为超临界阶段的Langmuir压力和应变常数；

p c

为CO₂相变临界压力，也是力学参数出现显著劣化的特征压力值。鉴于基质和裂隙吸附能力差异特性，二次吸附应变仅发生于基质系统。因此，

p L1 m ≠

p L2 m

，

p L1 f = p L2 f

，

ε L1 m ≠ ε L2 m

，

ε L1 f = ε L2 f

。为简化计算，裂隙Langmuir应变和压力分别统一改写为

ε Lf

和

p Lf

。表1列出了数值模拟参数取值范围，作为后续数值模拟参数取值的依据。

(2)基质力学参数性能劣化

三轴加载条件下实施强吸附性气体注入，低渗煤试样复合损伤因子变化呈现出典型的3个阶段过程^[29]，如图3所示。倪红坚等^[30]通过超临界态CO₂浸泡页岩试验发现，宏观力学损伤变化与页岩复合损伤因子呈现的3个阶段过程相吻合。

假设页岩–CO₂耦合过程中基质与岩体损伤程度一致，则弹性模量存在如下劣化过程^[29]：

E i = E 0 (1 − D)

(11)

式中：D为参数的劣化损伤因子，取值0～1；

E i

和

E 0

分别为弹性模量参数的动态值和初始值。式(11)是弹性损伤理论里弹性模量变量的核心定义，即通过损伤变量D量化材料微缺陷累积，以及对弹性模量的线性衰减效应。在经过广泛试验验证后成为连接微观损伤机制与宏观性能的关键桥梁。

考虑到有限元数值模拟过程中连续介质不存在无限制性损伤，因此以渐进式线性损伤定义参数变化，

E i

可设定为分段函数：

E i = {E 0 (p in ≤ p c) E s − E 0 p d − p c (p − p c) + E 0 (p c ＜ p in ≤ p d) E s (p d ＜ p in)

(12)

式中：

p c

为CO₂相变临界压力，也是力学参数出现显著劣化的特征压力值，

p c

= 7.4 MPa；

p d

为力学参数劣化达到稳定时的气体特征压力，

p d

= 12 MPa^[31]；

E 0

为未发生损伤的初始弹性模量；

E s

为损伤稳定后的弹性模量。无论弹性模量是否变化，泊松比增加都会扩大渗透率波动范围^[32]。因此，劣化泊松比也定义为

ν i = {ν 0 (p i ≤ p c) ν s − ν 0 p d − p c (p − p c) + ν 0 (p c ＜ p in ≤ p d) ν s (p d ＜ p in)

(13)

式中：

ν 0

为未发生损伤的初始泊松比，

ν s

为损伤稳定后的泊松比。对于劣化参数的变化范围，笔者做了相应的数据收集，如表2所示。需要说明，p_c为注入岩体气体达到压力平衡的目标值，而不是瞬态气压值。

3.2　气体流动控制方程

Knudsen数作为页岩气滑脱流动和多种气体扩散机制的分类依据，其表达式为

K n i = K B T d i 2 π δ 2 p i

(14)

式中：

K n i

为Knudsen数；

K B

为Boltzmann常数，

K B

= 1.38×10^－23 J/K；

d i

为多孔介质本质渗透率

k i ∞

和初始孔隙度

ϕ i 0

决定的等效水力直径，

d i =

32 k i ∞ / ϕ i 0

；T为储层温度(K)；

δ

为分子直径(m)；

p i

为孔隙压力(Pa)。

涵盖多尺度孔裂隙结构的页岩，在其等效水力直径d内存在气体的连续和非连续流动传输，两类流动控制方程可分别表示为

J c = f (K n) J D J dc = − D e ∇ (M C d)}

(15)

式中：

J c

为气体连续流动通量；

J D

为达西流动通量；M为气体相对分子质量；

J dc

为由时变有效扩散系数D_e控制的非连续流动质量通量，

C d

为气体浓度(g/m³)。对于管状孔隙，

f (K n) = 1 + 8 K n

。当

K n

<0.001时，连续流动等同于达西流，此时

J c = J D =

d 2 ρ g /

(32 μ ∇ p)

，kg/(m²·s)^[34]，

ρ g

为气体密度(kg/m³)，

μ

为气体动力黏度(Pa/s)，

∇ p

为孔隙压力梯度(Pa/m)。

气体时效扩散系数受气体浓度、分子位置、孔隙结构和扩散时间等多因素影响^[35]。综合考虑多种时间依赖型扩散方程，采用W. Zhao等^[36]所提出的有效解吸扩散系数计算模型：

D e (t) = D 0 τ 1 − e − ξ t t − (1 − e − ξ t) / ξ

(16)

式中：

D e (t)

为有效解吸扩散系数(m²/s)；

D 0

为初始扩散系数(m²/s)；

τ

为解吸时间参数；

ξ

为与孔径、气体密度和黏度相关的综合因子，当

ξ

趋向于无穷大时，

D e (t) = D 0 τ / t

。

气体分子在受限微管中的一次完整碰撞运动中，通过引入分子碰撞权重，可实现连续、非连续流动总质量通量的加权叠加^[7]：

J a = w s J c + (1 − w s) J dc

(17a)

式中：

w s

为连续流动的碰撞权重，计算式为

w s = K n / (K n + 2)

(17b)

进一步基于量纲等效转换原则：

k i a = − J i a μ /

(ρ g i ∇ p i)

(式中，

k ma

和

k fa

分别为基质系统和裂隙系统的等效视渗透率；

J ma

和

J fa

分别为基质系统和裂隙系统的连续、非连续流动总质量通量的加权叠加)，联立式(15)～(17)可得到基质、裂隙系统各自质量通量的等效视渗透率分别为

k ma = w s m d m 2 f (K n m) 32 + (1 − w s m) D 0 τ μ t p m

(18)

k fa = w s f d f 2 f (K n f) 32 + (1 − w s f) D 0 τ μ t p f

(19)

根据质量守恒定律，

∂ m ∂ t + ∇ (− k μ ρ g ∇ p) = Q s

，基质/裂隙系统中气体流动控制方程可分别表示为

∂ ∂ t (m m) + ∇ (− k ma μ ρ gm ∇ p m) = − Q s

(20a)

∂ ∂ t (m f) + ∇ (− k fa μ ρ gf ∇ p f) = Q s

(20b)

式中：

Q s

为基质与裂隙的质量交换量(kg/(m³·s))；，

m m

和

m f

分别为基质和裂隙系统中游离气和吸附气的质量总量：

m m = ρ g ϕ m + ρ ga ρ sm n ex m

(21a)

m f = ρ g ϕ f + ρ ga ρ sf n ex f

(21b)

式中：

ρ ga

为标准工况条件下气体密度(kg/m³)，

ρ sm

和

ρ sf

分别为基质系统和裂隙系统的质量密度，

ϕ m

和

ϕ f

分别为基质系统和裂隙系统的孔隙率，

n ex m

和

n ex f

分别为基质系统和裂隙系统的过剩吸附量。

3.3　基质–裂隙系统交叉耦合方程

考虑非线性流体效应和气体压差驱动，基质–裂隙系统间的气体质量交换为

Q s = τ ρ ga k ma μ (p m − p f)

(22)

式中：

τ

为多孔介质的传递形状因子，在REV尺度上模拟流体质量交换，取值为1 m^－2^[37]。

压差作用下具有不同力学性质的基质、裂隙，通过局部作用力交叉耦合产生的局部应变

Δ ε md

和

Δ ε fd

分别为

Δ ε md = p f − p m K m Δ ε fd = p m − p f K f}

(23)

当基质和裂隙达到压力平衡，即

p m = p f

，局部应变将消失，此时

Δ ε md = Δ ε fd = 0

。忽略气体吸附，基质和裂隙各自体积应变的变化量

Δ ε mv

和

Δ ε fv

可分别表示为

Δ ε mv = − 1 K m (σ ¯ − α p m)

(24)

Δ ε fv = − 1 K f (σ ¯ − β p f)

(25)

式中：

σ ¯

为岩体承受的3个主应力方向上压缩载荷

σ k k

的平均值，

σ ¯ = − σ k k / 3

。

吸附应变增量

Δ ε i s

是气压

p i

时膨胀应变与参考压力

p i 0

时的膨胀应变差，利用Langmuir型应变函数，基质和裂隙系统的吸附应变增量

Δ ε m s

和

Δ ε f s

分别为

Δ ε m s = {ε L1 m (p m p m + p L 1 − p m0 p m0 + p L1) (p m ≤ p c) ε L2 m (p m p m + p L2 − p m0 p m0 + p L2) (p c ＜ p m)

(26)

Δ ε f s = ε Lf (p f p f + p Lf − p f 0 p f 0 + p Lf)

(27)

应用伪一级动力学方程，考虑应变滞后效应的机械应变

ε i v

和吸附应变

ε i s

可分别表示^[38]为

ε i v = ε i v ∞ [1 − exp (− t / t m)]

(28)

ε i s = ε i ∞ s [1 − exp (− t / t s)]

(29)

式中：

ε i v ∞

和

ε i ∞ s

分别为极限机械应变和极限吸附应变，

t m

和

t s

分别为对应达到极限应变平衡所需时间。

基质系统的有效体积应变通过机械应变与裂隙–基质的局部应变之和，再减去吸附应变得到：

Δ ε me = Δ ε mv [1 − exp (− t t m)] − Δ ε m s [1 − exp (− t t s)] + Δ ε md = α p m − σ ¯ K m [1 − exp (− t t m)] − Δ ε m s ⋅ [1 − exp (− t t s)] − p f − p m K m

(30a)

式中：

1 − exp (− t t m)

为机械应变；

p f − p m K m

为基质局部应变；

Δ ε m s

为气体吸附诱发应变，计算式为

Δ ε m s = {ε L 1 m (p m) − ε L 1 m (p m 0) (p m ≤ p c) ε L 2 m (p m) − ε L 2 m (p m 0) (p c ＜ p m)

(30b)

与基质系统类似，裂隙系统的有效体积应变为

Δ ε fe = Δ ε fv [1 − exp (− t t m)] − Δ ε f s [1 − exp (− t t s)] + Δ ε fd = β p f − σ ¯ K f [1 − exp (− t t m)] − [ε Lf (p f) − ε Lf (p f 0)] [1 − exp (− t t s)] − p m − p f K f

(31)

式(31)中，等号右边第1项为机械应变，第2项为气体吸附诱发应变，第3项为裂隙局部应变。

多孔介质弹性体的孔隙度演化本质上随有效应变变化，满足如下表达式^[39]：

ϕ i ϕ i 0 = 1 + a biot ϕ i 0 Δ ε i e

(32)

式中：

a biot

为弹性体比奥系数，

Δ ε i e

为有效应变，

ϕ i

和

ϕ i 0

分别为动态孔隙率、初始孔隙率。

根据孔隙度–渗透率指数关系：

k i / k i 0 =

(ϕ i / ϕ i 0) M

，M为受岩石孔隙连通性、颗粒大小分布和微观结构等因素影响的孔隙度敏感指数，当用REV结构的几何参数定义孔隙度比值表征渗透率比值时，通常M = 3^[40]，但真实页岩心试验结果表明，M取值为25～35范围，某些情况下高达50^[41]。将式(30)，(31)分别代入式(32)，可得到基质和裂隙系统的渗透率比值

k m / k m 0

和

k f / k f 0

分别为

k m k m 0 = {1 − α (σ ¯ − α p m) [1 − exp (− t t m)] + p f − p m K m ϕ m 0 + α Δ ε m s [1 − exp (− t t s)] ϕ m 0} M

(33)

k f k f 0 = {1 − β (σ ¯ − β p f) [1 − exp (− t t m)] + p m − p f K f ϕ f 0 + β [ε Lf (p f) − ε Lf (p f 0)] [1 − exp (− t t s)] ϕ f 0} M

(34)

L. F. Fenton^[42]研究表明，整体渗透率可被定义为孔隙状和裂隙状渗透率的统计集合。假设基质和裂隙服从对数正态概率分布，且不存在控制渗流的主导流动路径，H. Singh等^[43]提出基质组分、裂隙组分和总体渗透率(

k bulk

)遵循如下数学关系：

ln (k bulk) = A ln (k m) + (1 − A) ln (k f)

(35)

式中：A为基质渗透率在整体渗透率中的贡献权重，取值0～1。从岩石的体积组成来看，A可解释为给定单位体积内基质的体积比例系数。将式(35)简化得到：

k bulk = k m A k f 1 − A

(36)

将式(33)和(34)代入式(36)，得到页岩整体渗透率演化模型为

k bulk k bulk ， 0 = k m A k f 1 − A k m0 A k f 0 1 − A = (k m k m 0) A (k f k f 0) 1 − A

(37)

D. Davudov和R. G. Moghanloo^[41]研究表明有效应力作用下孔隙连通性与页岩动态渗透率存在指数关系。将考虑连通性的修正项引入渗透率演化模型，得到如下表达式：

k bulk k bulk ， 0 = (k m k m 0) A (k f k f 0) 1 − A e L σ ¯

(38)

式中：L为连通性损失系数，在岩体弹性损伤范围内应用修正项

e L σ ¯

可提升模型的预测精度。当基质成分完全填充页岩时，A = 1，表示页岩渗透率比等同于基质系统渗透率比。事实上，很难分别量化基质和裂隙对整体渗透率的贡献，基于组分空间体积比例的方法，是一种折中方案。

4　数值模拟与试验匹配分析

收起

图4展示了全压条件下CO₂注入页岩的流–固耦合数值求解框架：以COMSOL自定义求解模块为核心，基于REFPROP精确提取不同相态CO₂关键物性参数；通过MATLAB构建数据交互接口，实现CO₂密度、黏度等物性参数与多物理场的协同耦合求解。其中，COMSOL平台主导完成2项关键计算：(1)基质–裂隙双系统流–固耦合，实现CO₂物性参数更新和固体变形的动态交互机制；(2) 更新受载岩体力学参数变化，对模型进行动态求解。

为确保模型具备较好收敛性，定义注气目标压力p_in为随时间变化的函数，即

p in = (p in − p a) + p in [1 − exp (− t / t d)]

(39)

式中：

p a

为参考压力(大气压)，

t d

为控制气体注入速率的特征时间。如图5所示，注气速率越快，基质–裂隙之间的压差增长速率越快，压差极值

Δ p max

越大。

4.1　气态CO₂/N₂注入

在气态CO₂注入条件下，吸附诱导的页岩膨胀应变和吸附量描述均遵循Langmuir型公式：

ε i s = ε L i p i / (p i + p L) n i s = n L i p i / (p i + p L)}

(40)

对于CO₂长时化学作用导致的应变滞后效应暂不做考虑，令

exp (− t / t m) = 0 exp (− t / t s) = 0}

(41)

式(41)表示数值模拟过程中，对吸附应变过程不做额外的应变平衡时间干预。由于CO₂与N₂注入过程中基质–裂隙系统的惯性效应及裂缝变形机制具有相似性，基于N₂注入页岩试验的渗透率动态响应数据^[44]，也可用于模型验证，模型相关参数取值如表3所示。

在实验室动态监测过程中应明确：当注气达到目标压力点时，不同时刻测得渗透率值可能呈现离散性。为建立试验数据与模拟结果的匹配对应关系，笔者采用颜色编码策略，将相同目标压力、不同模拟时间所得渗透率模拟值使用投影垂直条构建动态变化区间；将同等条件试验数据绘制于投影区间内，通过瞬时时间匹配关系预测试验所测时间点。

从图6中数值模拟结果可以看到：无论N₂还是气态CO₂，在注气压力增大到目标值过程中，尽管变化幅度不同，但页岩整体渗透率随时间始终呈现出类“U型”演化；将渗透率–时间演化数据转化为渗透率–目标压力投影面，发现试验数据恰好落在模型预测的范围内。将不同注入目标压力达成后的平衡渗透率模拟值与试验数据对比，发现除了图6(a)中N₂组高于6.8 MPa注入目标压力的渗透率出现明显回弹外，其余低压条件下，模拟与试验数据结果基本一致——都呈现逐渐下降趋势。不过，图6(b)中CO₂组试验数据的下降并没模拟结果显著，这说明岩心内部未达到对CO₂气体的吸附平衡，试验所测渗透率值并不代表最终实际平衡值。

对于渗透率演化趋势分析，笔者与H. H. Liu和J. L. Zhang^[45]的观点一致，即因低渗页岩在低压阶段由非线性流动效应主导，导致渗透率快速衰减；随后机械变形效应占据主导，促使渗透率出现回升现象。笔者更深入一步的研究表明，机械变形效应的强弱还受到气体吸附效应的影响，对CO₂吸附量–膨胀变形能力都更强的页岩显然回弹所需压力更大，如图6(b)所示。例如，同样4 MPa的目标注气压力，相比初始渗透率，N₂组的页岩动态渗透率最大损失25%，明显弱于CO₂组的最大损失55%。

4.2　超临界态CO₂注入

如图7(a)所示，超临界CO₂注入过程中页岩渗透率阶段变化比气态注入更明显，呈现出典型的“W”型演化模式；随着注入目标压力的增大，这种阶段性演化愈发明显。在更高目标压力(>10 MPa)条件下，非平衡阶段渗透率出现的回弹趋势更明显，这与L. Hou和D. Elsworth^[8]的3组试样试验结果相同；但渗透率平衡值呈现的却是缓慢下降趋势，这与图7(b)中L. Hou和D. Elsworth^[8]将气态和超临界态CO₂各自独立注入多组试样所测渗透率演化结果的平均值一致(L. Hou和D. Elsworth^[8]将渗透率演化划分为气态区“U型”变化，相态转换区“V型”变化，超临界区“U型”变化)。图7(b)中散点数据由L. Hou和D. Elsworth^[8]针对1～3号页岩样品，在气态、相变区和超临界态3个阶段分别进行独立CO₂注入试验后，测得的渗透率平均值构成。

将不同时刻数值模拟渗透率演化趋势与全压注入过程的页岩试样数据进行匹配，发现1号页岩注入压力低于6 MPa时，渗透率数据在200 s左右匹配，7 MPa时在150 s匹配，8 MPa时在1 000～1 300 s范围匹配；2和3号页岩的渗透率演化趋势与模拟时间105 s时结果相似，但各自达到注气目标压力的时间存在显著差异；2和3号页岩所测部分数据大致在200和1 500 s时得到匹配，这表明不同目标压力注入过程中实测渗透率数据非岩心完全达到平衡后的稳定值，而是处于非平衡状态下的瞬态值。同时，可以看到10⁵ s模拟所得渗透率“W型”演化曲线与L. Hou和D. Elsworth^[8]试验所测“双U＋V型”相似。但必须谨慎说明，模型与试验数据匹配的时间点(以s为单位)受多参数影响，其精确性不能完全等同于实际试验，仅可作为时间尺度上的预测值。

4.3　全压范围CO₂注入页岩渗透率演化模式

常围压条件下，CO₂由气态到超临界态注入页岩过程中，结合多种应变机制演化情况，页岩的渗透率演化规律和驱动因素可描述为：

阶段I：外部气体由低压逐渐注入致密低渗岩体中，低气压气体滑脱效应占主导，视渗透率增大，导致渗透率比值呈现下降趋势；但同时气体迅速进入渗透性更强的裂隙系统，软性裂隙率先做出开度增大响应，其渗透率增加；而硬性基质受裂隙挤压影响，其渗透率暂时降低。前期注气过程中，基质–裂隙的力学变形对页岩整体渗透率影响十分有限，因此受气体滑脱效应影响，页岩整体渗透率比值呈现显著降低趋势。

阶段II：两孔隙介质系统内孔隙压力逐渐上升，由于基质具有发达的微/纳米孔和强大的吸附能力，其显著的体积膨胀将挤压裂隙使其收缩闭合，基质组分渗透率显著增加，而裂隙组分渗透率下降；考虑到基质组分体积占比优势，以及体积膨胀对孔隙连通性的连通强化，从而导致整体渗透率表现出与基质组分相似的阶段性增大。

阶段III：随注入压力持续上升，基质组分吸附逐渐达到饱和，但基质和裂隙之间压差诱发的局部应变显著增大，导致基质、裂隙各自的有效应变降低，从而使整体的有效应变降低，最终导致整体渗透率出现一定程度的降低。

阶段IV：虽然压差导致局部应变逐渐消失，但CO₂与页岩矿物组分的持续化学反应，导致基质微孔变化并出现了短暂的二次吸附现象；在整体力学性质显著劣化情况下，二次吸附膨胀应变在抵消了一部分机械压缩应变后保持稳定，最终平衡渗透率虽有所增加，但未能恢复到初始值大小。

4.4　影响效应关键参数敏感性分析

压差传递相关的特征时间

t d

影响基质/裂隙局部应变竞争程度，平衡时间

t s

影响膨胀和机械应变之间的相互作用，对渗透率演化都有影响。因此，针对基质初始杨氏模量和泊松比(这里用E₀，

ν s0

代指E_m0，

ν m0

)，压差传递特征时间

t d

和应变平衡时间

t s

分别开展模型参数敏感性分析。基于表1所列实测参数劣化范围，模拟杨氏模量20%，30%和40%劣化损失，泊松比20%，40%和60%劣化增大，进行力学参数劣化量化分析，其余输入参数如表4中所列。

如图8(a)，(b)所示，基质力学参数劣化的加剧明显促进了阶段III整体渗透率的下降，并抑制了阶段IV渗透率的反弹；基质组分和裂隙组分各自渗透率比值的变化范围表明，力学参数经过劣化作用后，相同有效应力下软性裂隙对力学响应更为敏感。图8(c)中，

t d

越小意味着越大的气体注入速率，高注气速率明显放大了阶段I～III的渗透率变化范围；图8(d)中，

t s

越大表明吸附应变耗时越长，

t m

不变情况下应变平衡时间的延长明显扩大了阶段II到III的渗透率变化范围。

由图5可知，注入压力达到临界压力大约需要1 000 s，以这一节点为分界线，观察图8中的子图，发现超临界CO₂区域渗透率表现出持续时间更长、更显著的下降趋势。显然，超临界CO₂–页岩相互作用显著强化了渗透率应力敏感性，这与J. P. Zhou等^[19]相关试验研究结论是完全一致的。

4.5　渗透率模型及数值模拟方法分析

为何实验室高压CO₂被注入岩心试验过程中，有些试验所测渗透率演化后期会出现回弹增大和减小2种情况呢？笔者认为，需要对注入岩心的CO₂相态和流–固耦合作用机制进行区分。如图9(a)所示，在弹性变形假设前提下，外部应力和注气补能所诱发的裂隙–基质竞争变形具有可逆性，当高压气态CO₂被注入时，在有效应力降低和基质显著膨胀变形作用下，达到气体平衡的页岩渗透率值会显著增大。当超临界态CO₂被长期注入，一方面，矿化反应和气体吸附作用会引起页岩力学性质劣化、产生不可逆的塑性变形以及孔裂隙等微结构损伤，受孔裂隙连通性变化渗透率可能出现明显突增现象，甚至出现类似于压裂的渗透率阶跃行为，导致渗透率平衡值出现显著增大；另一方面，矿化反应也会引发缓慢矿物溶蚀沉淀并堵塞流体运移通道，再加上塑性劣化岩体渐进式的软化压实，也可能导致渗透率平衡值降低，但这一过程可能相当漫长。

如图9(b)所示，模块化定义相异组分材料的重叠连续非均质岩石，渗透率数学模型通常为连续函数，通过参数设定确保流–固耦合求解结果稳定，其多物理场耦合求解表示了基质与裂隙之间的竞争性变形关系，自洽适应性使整个岩体在内部组分差异变形情况下也处于相对稳定状态。

然而，对于离散型介质，其多物理场耦合求解过程具有连接性和扩展性，当应变能积累至临界阈值时，岩体中预先存在的弱相(如裂隙)作为应力集中区首先发生局部破裂，这些局部破裂扩展通过应力场的动态调整形成级联效应。当释放的应变能量传递路径与岩体结构网络中的“高连通性区域”相遇时，就会触发指数级裂隙扩展，导致整体宏观破裂和结构失稳。为确保求解的稳定，渗透率数学模型通常定义为分段函数，从而能够更灵活地反映渗透率演化的阶段性变化，如：注气阶段渐进性蠕变分数阶损伤渗透率模型^[46]，或结构破坏失效突变渗透率模型^[47]等。

对比来看，双连续介质方法在隐式处理裂缝–基质相互作用方面，比当前离散型裂缝网络(discrete fracture network，DFN)和离散裂缝矩阵(discrete fracture matrix，DFM)方法更为便捷，但在裂缝作为目标区域主要流动通道的情况下，无法显式捕捉裂缝几何形态及变化^[48]。因此，耦合地质力学与多相流问题仍是数学模型、数值模型及离散化方案领域的难点，要准确捕捉裂缝–基质相互作用下的流体动力学行为，需要根据研究对象建立合适的理论模型并选用恰当的求解方案。表5为笔者参考相关文献总结的几类流–固耦合渗透率数学模型，其模型特征直观体现了渗透率建模的基本思路。

5　影响效应独立解耦分析

收起

节4中，结合基质力学参数劣化、二次吸附及应变滞后效应构建的模型，与试验数据吻合良好。为了研究单一因素对渗透率演化的影响，笔者进行了3种情况解耦计算：(1)没有力学参数劣化的情况；(2) 没有二次吸附及分阶段应变的情况；(3) 没有应变滞后效应的情况。

5.1　解耦基质力学参数劣化

如图10(a)所示，不考虑基质力学性质劣化，整体渗透率仍呈现出四阶段演化规律。但是解耦基质力学参数劣化的渗透率比值，其变化范围(1～0.8)比图8(b)中耦合情况下的变化范围(1～0.98)更大，这表明基质的力学劣化并未放大渗透率的波动范围，反而对其起到了缓冲作用；图10(c)中的位移梯度云图则表明，基质与裂隙之间仍然存在显著的竞争应变作用。

5.2　解耦基质系统两阶段吸附应变

(1)Langmuir型吸附模型应用

如图11(a)，(b)所示，当利用Langmuir吸附模型来描述吸附能力时，页岩整体渗透率虽经历了四阶段演化，但基质和裂隙的有效应变均呈现出“U形”变化趋势。图11(c)中位移梯度云图与图10(c)对比可以明显看出：由Langmuir吸附函数描述涵盖低压到高压注入场景过程中的真实吸附量，将低估CO₂对基质微孔、中孔大小的改性作用，使得高压条件下被Langmuir模型简化的基质吸附膨胀作用被弱化，从而导致基质–裂隙系统间因吸附–变形机制差异导致的应变竞争效应显著削弱。

(2)解耦两阶段吸附膨胀变形

如图12(a)所示，页岩整体渗透率仍然呈现出四阶段演化趋势，最明显的变化是渗透率终值对基质p_Lm十分敏感。如图12(b)所示，单一阶段Langmuir吸附应变过程，p_Lm的值决定了气体注入后最终渗透率是否会大于初始渗透率。试验所测p_Lm总是倾向于小于6 MPa^[53]，这意味着气体注入达到平衡后，页岩渗透率终值会大于初始值。显然，这与试验所测结果相反，表明单Langmuir型方程描述吸附应变不符合实际情况。

相比之下，使用两阶段应变函数可以更准确地预测渗透率演化的整个过程，这也反映出CO₂相变对于吸附应变的确存在显著影响。此外，图12(c)中位移梯度随时间变化并不显著，表明两阶段吸附应变对基质–裂隙系统间的竞争应变影响显著。

5.3　解耦应变滞后效应

如图13(a)所示，当忽略应变滞后效应时，渗透率呈现多阶段演化特征。于阶段IV缓慢增大之后，又出现下降–稳定的过程，这可能是因为吸附应变已达极限后，力学性质劣化的页岩仍在经历缓慢的压缩应变累积，这一情况通过图13(b)中长时间显著的位移梯度变化得到佐证。同时，由于吸附变形比机械应变响应更早，因此基质、裂隙以及页岩整体的渗透率均存在短期上升现象。相比于图9(b)全耦合，无应变滞后效应的基质–裂隙竞争应变在注气较短时间内表现的更为显著(见图13(c))。

6　结论

收起

针对当前渗透率模型难以准确描述CO₂由低压气态到超临界态全压注入页岩过程中的渗透率演化过程，运用理论建模与数值模拟方法开展CO₂全压注入页岩渗透率演化过程研究，并揭示控制其阶段演化的基本机制，得出如下主要结论：

(1)受基质渗透率和裂隙渗透率演化之间复杂相互作用影响，页岩整体渗透率是一个耦合参数。在双孔隙均质材料模型中预测的渗透率演化反映的是组分渗透率，而在双弹性介质材料并行耦合框架内，页岩渗透率可表征为各组分渗透率加权后的耦合渗透率。基于组分渗透率加权确定的整体渗透率演化模型，能够有效刻画CO₂注入过程中渗透率基本演化模式。

(2)在弹性力学与连续介质理论框架下，基于力–化耦合作用机制构建的渗透率演化模型，系统阐明了CO₂相变诱发的渗透率四阶段演化过程，各阶段的渗透率演化规律及主控因素为：阶段I，渗透率比值下降，主要由气体滑脱效应占据主导，基质–裂隙的力学变形影响有限；阶段II，渗透率增加，主要由吸附膨胀控制；阶段III，渗透率下降，主要由竞争性局部应变驱动；阶段IV，渗透率恢复增加，主要归因于二次吸附膨胀以及基质力学性能劣化。

(3)裂隙与基质的力学差异显著扩大了渗透率演化范围。杨氏模量降低会减弱渗透率变化幅度，延缓演化过程，而泊松比的增加则与之相反；机械应变和吸附应变同步响应影响渗透率的演化模式和渗透率平衡终值，应变滞后效应仅延长了各阶段演化持续时间，使渗透率呈现出更为直观的多阶段演化特征；应用Langmuir吸附模型对渗透率演化趋势分析影响不大，但会低估平衡后渗透率终值；双Langmuir型吸附应变函数应用，显著提高了超临界CO₂区域渗透率演化预测的准确性。

当前，通过弹塑性损伤本构模型定义页岩力学损伤过程仍存在局限性，未来需考虑连续–非连续介质的耦合求解，以期更精准表征页岩复杂破坏机制。此外，笔者正构建基于原位CT可视化的多物理场耦合试验平台，致力于解构渗透率动态演化各阶段对应的微观结构变形规律，从而揭示页岩多尺度力学响应对渗透性演化的控制机制。

基金

收起

国家自然科学基金面上项目(52574238; 52274191)
河南理工大学河南省瓦斯地质与瓦斯治理重点实验室开放课题(WS2024B10)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

CHANG

，LIN

S S

，YANG

C H

，et al. A critical review of ScCO₂-enhanced gas recovery and geologic storage in shale reservoirs[J]. Gas Science and Engineering，2024，125：205317.

[2]

曲鸿雁，彭岩，刘继山，等. 气体吸附对页岩裂缝表观渗透率和页岩气采收率的影响[J]. 中国科学：技术科学，2018，48(8)：891–900.

(QU

Hongyan

，PENG

Yan

，LIU

Jishan

，et al. Impact of gas adsorption on apparent permeability of shale fracture and shale gas recovery rate[J]. Scientia Sinica Technologica，2018，48(8)：891–900.(in Chinese))

[3]

GAO

，LIU

J S

，LEONG

Y K

，et al. A review of swelling effect on shale permeability：assessments and perspectives[J]. Energy and Fuels，2023，37(5)：3 488–3 500.

[4]

CUI

G L

，LIU

J S

，WEI

M Y

，et al. Why shale permeability changes under variable effective stresses：New insights[J]. Fuel，2018，213：55–71.

[5]

陈帅，李波波，李建华，等. 不同应力条件下页岩表观渗透率模型试验研究[J]. 煤炭科学技术，2021，49(7)：169–178.

(CHEN

Shuai

，LI

Bobo

，LI

Jianhua

，et al. Experimental study on shale apparent permeability model under different stress conditions[J]. Coal Science and Technology，2021，49(7)：169–178.

[6]

SHI

，LIU

J S

，WANG

X M

，et al. A critical analysis of shale laboratory permeability evolution data[J]. Energy，2021，236：121405.

[7]

W R

，JU

，ZHOU

H W

，et al. Novel shale permeability model coupling sorption-induced differential deformation and multiple flow regimes[J]. Energy and Fuels，2022，36(19)：11 777–11 786.

[8]

HOU

，ELSWORTH

. Mechanisms of tripartite permeability evolution for supercritical CO₂ in propped shale fractures[J]. Fuel，2021，292：120188.

[9]

张宏学，刘卫群. 非平衡解吸状态下页岩气储层渗透率演化机制[J]. 岩土力学，2021，42(10)：2 696–2 704.

(ZHANG

Hongxue

，LIU

Weiqun

. Permeability evolution mechanism of shale gas reservoir in non-equilibrium desorption state[J]. Rock and Soil Mechanics，2021，42(10)：2 696–2 704.(in Chinese))

[10]

肖智勇，王长盛，王刚，等. 基质–裂隙相互作用对渗透率演化的影响：考虑基质变形和应力修正[J]. 岩土工程学报，2021，43(12)：2 209–2 219.

(XIAO

Zhiyong

，WANG

Changsheng

，WANG

Gang

，et al. Influences of matrix-fracture interaction on permeability evolution：considering matrix deformation and stress correction[J]. Chinese Journal of Geotechnical Engineering，2021，43(12)：2 209–2 219.

[11]

RONG

T L

，LIU

K L

，GUAN

. Evolution of internal swelling coefficient under different geometric conditions and its influence on gas drainage[J]. Natural Resources Research，2022，31(2)：1 027–1 044.

[12]

XIE

J Q

，XIONG

，TAN

Y L

，et al. Effects of anisotropic permeability evolution on shale gas production：an internal swelling factor model[J]. Energy and fuels，2022，36(2)：771–785.

[13]

肖智勇，孙小翔，王刚，等. 气体压差影响下的煤渗透率非平衡演化全过程模型[J]. 岩土工程学报，2024，47(2)：355–364.

(XIAO

Zhiyong

，SUN

Xiaoxiang

，WANG

Gang

，et al. Model for non-equilibrium evolution coal permeability of whole process under influences of gas pressure difference[J]. Chinese Journal of Geotechnical Engineering，2024，47(2)：355–364.(in Chinese))

[14]

CHENG

W X

，CUI

G L

，TAN

Y L

，et al. A multi-layer nanocased model to explain the U-shaped evolution of shale gas permeability at constant confining pressure[J]. Fuel，2024，359：130478.

[15]

HOU

，ELSWORTH

，GENG

X Y

. Swelling and embedment induced by sub- and supercritical-CO₂ on the permeability of propped fractures in shale[J]. International Journal of Coal Geology，2020，225：103496.

[16]

汤积仁，卢义玉，陈钰婷，等. 超临界CO₂作用下页岩力学特性损伤的试验研究[J]. 岩土力学，2018，39(3)：797–802.

(TANG

Jiren

，LU

Yiyu

，CHEN

Yuting

，et al. Experimental study of damage of shale mechanical properties under supercritical CO₂[J]. Rock and Soil Mechanics，2018，39(3)：797–802.(in Chinese))

[17]

YANG

Y F

，LI

Y W

，YAO

，et al. Dynamic pore-scale dissolution by CO₂-saturated brine in carbonates：Impact of homogeneous versus fractured versus vuggy pore structure[J]. Water Resources Research，2020，56(4).

https://doi.org/10.1029/2019WR026112

[18]

LYU

，WANG

K X

，HU

C E

，et al. Effects of supercritical CO₂/water imbibition under dynamic pressures on shale mechanics and acoustic emission characteristics[J]. Fuel，2022，321：124087.

[19]

ZHOU

J P

，TIAN

S F

，ZHOU

，et al. Experimental investigation on the influence of sub- and super-critical CO₂ saturation time on the permeability of fractured shale[J]. Energy，2020，191：116574.

[20]

FENG

，KANG

，SUN

Z d

，et al. Effects of supercritical CO₂ adsorption on the mechanical characteristics and failure mechanisms of shale[J]. Energy，2019，173：870–882.

[21]

WANG

L L

，LONG

Z J

，SONG

，et al. Supercritical CO₂ adsorption and desorption characteristics and pore structure controlling mechanism of tectonically deformed coals[J]. Fuel，2022，317：123485.

[22]

M Y

，LIU

B G

，CHU

Z F

，et al. Permeability，deformation characteristics，and damage constitutive model of shale under triaxial hydromechanical coupling[J]. Bulletin of Engineering Geology and the Environment，2022，81(3)，https://doi.org/10.1007/s10064-022-02584-4.

[23]

ZHANG

S W

，LIU

J S

，WEI

M Y

，et al. Coal permeability maps under the influence of multiple coupled processes[J]. International Journal of Coal Geology，2018，187：71–82.

[24]

FAN

，LIU

S M

. Numerical prediction of in situ horizontal stress evolution in coalbed methane reservoirs by considering both poroelastic and sorption induced strain effects[J]. International Journal of Rock Mechanics and Mining Sciences，2018，104：156–164.

[25]

WEI

M Y

，LIU

J S

，ELSWORTH

，et al. Impact of equilibration time lag between matrix and fractures on the evolution of coal permeability[J]. Fuel，2021，290：120029.

[26]

WEI

M Y

，LIU

J S

，ELSWORTH

，et al. Influence of gas adsorption induced non-uniform deformation on the evolution of coal permeability[J]. International Journal of Rock Mechanics and Mining Sciences，2019，114：71–78.

[27]

ZHI

，ELSWORTH

，LIU

L Y

. W-shaped permeability evolution of coal with supercritical CO₂ phase transition[J]. International Journal of Coal Geology，2019，211：103221.

[28]

ZHOU

S W

，XUE

H Q

，NING

，et al. Experimental study of supercritical methane adsorption in Longmaxi shale：Insights into the density of adsorbed methane[J]. Fuel，2018，211：140–148.

[29]

白鑫，王登科，田富超，等. 三轴应力加卸载作用下损伤煤岩渗透率模型研究[J]. 岩石力学与工程学报，2021，40(8)：1 536–1 546.

(BAI

Xin

，WANG

Dengke

，TIAN

Fuchao

，et al. Permeability model of damaged coal under triaxial stress loading-unloading[J]. Chinese Journal of Rock Mechanics and Engineering，2021，40(8)：1 536–1 546.(in Chinese))

[30]

倪红坚，郭兴，丁璐，等. 超临界二氧化碳浸泡对页岩力学性质影响的试验[J]. 中国石油大学学报：自然科学版，2019，43(2)：77–84.

(NI

Hongjian

，GUO

Xing

，DING

，et al. Experiment on mechanical properties of shale soaked with supercritical carbon dioxide[J]. Journal of China University of Petroleum：Natural Science，2019，43(2)：77–84.(in Chinese))

[31]

YIN

，ZHOU

J P

，XIAN

X F

，et al. Experimental study of the effects of sub- and super-critical CO₂ saturation on the mechanical characteristics of organic-rich shales[J]. Energy，2017，132：84–95.

[32]

WANG

，LIU

S M

，ZHAO

Y X

. Modeling of permeability for ultra-tight coal and shale matrix：A multi-mechanistic flow approach[J]. Fuel，2018，232：60–70.

[33]

LYU

，LONG

，RANJITH

P G

，et al. Experimental investigation on the mechanical properties of a low-clay shale with different adsorption times in sub-/super-critical CO₂[J]. Energy，2018，147：1 288–1 298.

[34]

LIU

，WANG

J G

，GAO

，et al. Impact of micro- and macro-scale consistent flows on well performance in fractured shale gas reservoirs[J]. Journal of natural gas science and engineering，2016，36：1 239–1 252.

[35]

赵伟，王凯，周爱桃，等. 扩散主控型煤层定义、特征及瓦斯促抽路径分析[J]. 煤炭学报，2025，50(7)：3 280–3 292.

(ZHAO

Wei

，WANG

Kai

，ZHOU

Aitao

，et al. Definition，characteristics and analysis of gas extraction improvement method in diffusion-dominated coal seams[J]. Journal of China Coal Society，2025，50(7)：3 280–3 292.(in Chinese))

[36]

ZHAO

，CHENG

Y P

，JIANG

H N

，et al. Modeling and experiments for transient diffusion coefficients in the desorption of methane through coal powders[J]. International Journal of Heat and Mass Transfer，2017，110：845–854.

[37]

MOHAGHEGHIAN

，HASSANZADEH

，CHEN

. Matrix-fracture transfer shape factor for modeling multimechanisitc multicomponent shale gas flow[J]. International Journal of Heat and Mass Transfer，2020，158：120022.

[38]

ZHAO

，WANG

，LIU

S M

，et al. Asynchronous difference in dynamic characteristics of adsorption swelling and mechanical compression of coal：Modeling and experiments[J]. International Journal of Rock Mechanics and Mining Sciences，2020，135：104498.

[39]

LIU

J S

，CHEN

Z W

，ELSWORTH

，et al. Interactions of multiple processes during CBM extraction：A critical review[J]. International Journal of Coal Geology，2011，87(3/4)：175–189.

[40]

ZHENG

J T

，ZHENG

L G

，LIU

H H

，et al. Relationships between permeability，porosity and effective stress for low-permeability sedimentary rock[J]. International Journal of Rock Mechanics and Mining Sciences，2015，78：304–318.

[41]

DAVUDOV

，MOGHANLOO

R G

. Impact of pore compressibility and connectivity loss on shale permeability[J]. International Journal of Coal Geology，2018，187：98–113.

[42]

FENTON

L F

. The sum of log-normal probability distributions in scatter transmission systems[J]. Communications Systems Ire Transactions on，1960，8(1)：57–67.

[43]

SINGH

，JAVADPOUR

，ETTEHADTAVAKKOL

，et al. Nonempirical apparent permeability of shale[J]. SPE Reservoir Evaluation and Engineering，2014，17(3)：414–424.

[44]

ZHAI

W B

，WU

，LIU

X Y

，et al. The seepage model for CO₂ in shale considering dynamic slippage，effective stress and gas adsorption[J]. Scientific Reports，2024，14(1)：30 623–30 697.

[45]

LIU

H H

，ZHANG

J L

. An efficient laboratory method to measure the combined effects of Knudsen diffusion and mechanical deformation on shale permeability[J]. Journal of Contaminant Hydrology，2020，232：103652.

[46]

ZHOU

H W

，ZHANG

，WANG

X Y

，et al. Effects of matrix-fracture interaction and creep deformation on permeability evolution of deep coal[J]. International Journal of Rock Mechanics and Mining Sciences，2020，127：104236.

[47]

ZHU

W C

，WEI

C H

，LI

，et al. Numerical modeling on destress blasting in coal seam for enhancing gas drainage[J]. International Journal of Rock Mechanics and Mining Sciences，2013，59：179–190.

[48]

HAWEZ

H K

，SANAEE

，FAISAL

N H

. A critical review on coupled geomechanics and fluid flow in naturally fractured reservoirs[J]. Journal of Natural Gas Science and Engineering，2021，95：104150.

[49]

SHI

J Q

，DURUCAN

. Drawdown induced changes in permeability of coalbeds：a new interpretation of the reservoir response to primary recovery[J]. Transport in Porous Media，2004，56(1)：1–16.

[50]

WEI

M Y

，LIU

J S

，FENG

X T

，et al. Evolution of shale apparent permeability from stress-controlled to displacement-controlled conditions[J]. Journal of Natural Gas Science and Engineering，2016，34：1 453–1 460.

[51]

ZENG

，LIU

J S

，LI

，et al. Evolution of shale permeability under the influence of gas diffusion from the fracture wall into the matrix[J]. Energy and fuels，2020，34(4)：4 393–4 406.

[52]

WANG

C Y

，ZHANG

D M

，TAN

D J

，et al. Study on the evolution of limestone damage and permeability under the tunnel driving stress path[J]. Construction and Building Materials，2024，425：136050.

[53]

HUI

，PAN

，LUO

P Y

，et al. Effect of supercritical CO₂ exposure on the high-pressure CO₂ adsorption performance of shales[J]. Fuel，2019，247：57–66.

2026年第45卷第2期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.3724/1000-6915.jrme.2025.0505

接收时间：2025-07-17
首发时间：2026-06-18
出版时间：2026-02-01

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2025-07-17
修回日期：2025-10-13

基金

National Natural Science Foundation of China(52574238; 52274191)

国家自然科学基金面上项目(52574238; 52274191)

Open Project of Henan Key Laboratory for Gas Geology and Gas Control, Henan Polytechnic University(WS2024B10)

河南理工大学河南省瓦斯地质与瓦斯治理重点实验室开放课题(WS2024B10)

作者信息

^1.河南理工大学安全科学工程学院，河南焦作　454000

^2.河南理工大学河南省瓦斯地质与瓦斯治理重点实验室–省部共建国家重点实验室培育基地，河南焦作　454000

^3.河南理工大学煤矿瓦斯灾害防治应急管理部重点实验室，河南焦作　454000

通讯作者:

^* 王登科（1980–），现任教授，主要从事安全科学与工程方面的研究工作。E-mail：wdk@hpu.edu.cn

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/yslxygcxb/CN/10.3724/1000-6915.jrme.2025.0505

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT

Parameter	$p L1 m$ /MPa	$p L2 m$ /MPa	$ε L1 m$ /%	$ε L2 m$ /%	$ε Lf$ /%	$p c$ /MPa
Value range	1～2.3^[28]	4～7^[28]	0.08～0.12^[27]	0.12～0.2	0.04～0.08	7.2～8.0
Value	2	6	0.08	0.12	0.04	7.4

Parameter

p L1 m

/MPa

p L2 m

/MPa

ε L1 m

ε L2 m

ε Lf

p c

/MPa

Value range

1～2.3^[28]

4～7^[28]

0.08～0.12^[27]

0.12～0.2

0.04～0.08

7.2～8.0

Value

0.08

0.12

0.04

7.4

Parameter	Degree of deterioration/%	Literature sources
Young's modulus	－24.60～26.79 －29.28～33.15 －20.23～37.39	Q. Lyu et al.^[18] H. Yin et al.^[31] Q. Lyu et al.^[33]
Poisson's ratio	29.41～58.82	Q. Lyu et al.^[33]

Parameter

Degree of deterioration/%

Literature sources

Young's modulus

－24.60～26.79
－29.28～33.15
－20.23～37.39

Q. Lyu et al.^[18]
H. Yin et al.^[31]
Q. Lyu et al.^[33]

Poisson's ratio

29.41～58.82

Q. Lyu et al.^[33]

Initial porosity of matrix $ϕ mo$	Initial porosity of fracture $ϕ fo$	Maximum adsorption capacity of matrix part $n max m$ /(m³·kg)	Maximum adsorption capacity of fracture part $n max f$ /(m³·kg)	Biot coefficient of matrix $α$	Biot coefficient of fracture $β$	Initial permeability of matrix part $k m0$ /m²	Initial permeability of fracture part $k f0$ /m²	Sensitivity exponent M	Injected target pressure $p in$ /MPa
0.04	0.04	0.12	0.05	0.65	0.8	5×10^－19	1×10^－18	3	14

Initial porosity of matrix

ϕ mo

Initial porosity of fracture

ϕ fo

Maximum adsorption capacity of matrix part

n max m

/(m³·kg)

Maximum adsorption capacity of fracture part

n max f

/(m³·kg)

Biot coefficient of matrix

α

Biot coefficient of fracture

β

Initial permeability of matrix part

k m0

/m²

Initial permeability of fracture part

k f0

/m²

Sensitivity exponent M

Injected target pressure

p in

/MPa

0.04

0.12

0.05

0.65

0.8

5×10^－19

1×10^－18

Basic expression of model	Parameter interpretation	Literature sources
$k k 0 = exp (− 3 C f Δ σ − α biot Δ p)$	$Δ σ$ is the stress change; $C f$ is the fracture compression coefficient,1/MPa	J. Q. Shi and S. Durucan^[49]
$k k 0 = (1 + α biot ϕ 0 Δ ε e) 3 f (K n)$	$Δ ε e$ is the effective strain change, $f (K n)$ is the slip correction coefficient related to the Knudsen number	M. Y. Wei et al.^[50]
$k k 0 = exp [3 (1 K − 1 K f) (Δ σ − Δ p) + Δ ε lfs] (l l 0) 2 f (K n) f (K n0)$	L is the fracture length; $ε lfs$ is the fracture strain determined by the volume fraction of adsorbed gas invading the shale, gas diffusion coefficient, and diffusion distance, defined as a piecewise function	J. Zeng et al.^[51]
$k k 0 = {exp (− 3 C f Δ σ e) (Elastic stage) exp (− 3 D λ D C f Δ σ e) (Plastic stage)$	$Δ σ e$ is the effective stress change; $λ D$ is the influence coefficient of rock damage on permeability; $D$ is the damage coefficient, defined as a piecewise function based on tensile and shear strains	C. Y. Wang et al.^[52]
$k k 0 = (ϕ ϕ 0) 3 exp (α k D)$	$α k$ is the fitted value of the influence coefficient of rock mass damage on permeability	W. C. Zhu et al.^[47]
$k k 0 = (ϕ ϕ 0) M + 1 exp (η z 0 − η z)$	M is the stress sensitivity constant obtained from experimental fitting,z is the initial coordination number of pores and fractures,and $η$ is a constant related to the fractal dimension and M.	D. Davudov and R. G. Moghanloo^[41]

Basic expression of model

Parameter interpretation

Literature sources

k k 0 = exp (− 3 C f Δ σ − α biot Δ p)

Δ σ

is the stress change;

C f

is the fracture compression coefficient,1/MPa

J. Q. Shi and S. Durucan^[49]

k k 0 = (1 + α biot ϕ 0 Δ ε e) 3 f (K n)

Δ ε e

is the effective strain change,

f (K n)

is the slip correction coefficient related to the Knudsen number

M. Y. Wei et al.^[50]

k k 0 = exp [3 (1 K − 1 K f) (Δ σ − Δ p) + Δ ε lfs] (l l 0) 2 f (K n) f (K n0)

L is the fracture length;

ε lfs

is the fracture strain determined by the volume fraction of adsorbed gas invading the shale, gas diffusion coefficient, and diffusion distance, defined as a piecewise function

J. Zeng et al.^[51]

k k 0 = {exp (− 3 C f Δ σ e) (Elastic stage) exp (− 3 D λ D C f Δ σ e) (Plastic stage)

Δ σ e

is the effective stress change;

λ D

is the influence coefficient of rock damage on permeability;

D

is the damage coefficient, defined as a piecewise function based on tensile and shear strains

C. Y. Wang et al.^[52]

k k 0 = (ϕ ϕ 0) 3 exp (α k D)

α k

is the fitted value of the influence coefficient of rock mass damage on permeability

W. C. Zhu et al.^[47]

k k 0 = (ϕ ϕ 0) M + 1 exp (η z 0 − η z)

M is the stress sensitivity constant obtained from experimental fitting,z is the initial coordination number of pores and fractures,and

η

is a constant related to the fractal dimension and M.

D. Davudov and R. G. Moghanloo^[41]

Experiment type	Environment temperature T/K	Initial porosity of matrix $ϕ m0$	Initial porosity of fracture $ϕ f0$	Maximum adsorption capacity of matrix part $n max m$ /(m³·kg)	Maximum adsorption capacity of fracture part $n max f$ /(m³·kg)	Poisson's ratio $ν$	Young's modulus of matrix part $E m$ /GPa	Young's modulus of fracture part $E f$ /GPa
N₂ group^[45]	298.15	0.04	0.04	0.015	0.005 0	0.15	5	2
CO₂ group^[44]	313.00	0.05	0.05	0.080	0.008 7	0.25	8	6
Supercritical CO₂ group^[8]	320.00	0.04	0.04	0.030	0.015 0	0.25	10	8
Experiment type	Biot coefficient of matrix $α$	Biot coefficient of fracture $β$	Initial permeability of matrixi $k m 0$ /m²	Initial permeability of fracture $k f0$ /m²	Langmuir strain constant of matrix part $ε Lm$ /%	Langmuir strain constant of fracture part $ε Lf$ /%	Langmuir pressure constant of fracture part $p Lf$ /MPa	Langmuir pressure constant of matrix part $p Lm$ /MPa
N₂ group^[45]	0.30	0.60	3×10^－18	1×10^－17	0.4	0.2	4	2
CO₂ group^[44]	0.65	0.85	1×10^－19	1×10^－17	0.8	0.4	8	2
Supercritical CO₂ group^[8]	0.70	0.90	5×10^－19	1×10^－17	$0.4 (ε L 1 m) ， 1.2 (ε L 2 m)$	$0.4 (ε L f)$	$4 (p L f)$	$2 (p L 1 m) ， 8 (p L 2 m)$

Experiment type	Environment temperature T/K	Initial porosity of matrix $ϕ m0$	Initial porosity of fracture $ϕ f0$	Maximum adsorption capacity of matrix part $n max m$ /(m³·kg)	Maximum adsorption capacity of fracture part $n max f$ /(m³·kg)	Poisson's ratio $ν$	Young's modulus of matrix part $E m$ /GPa	Young's modulus of fracture part $E f$ /GPa
N₂ group^[45]	298.15	0.04	0.04	0.015	0.005 0	0.15	5	2
CO₂ group^[44]	313.00	0.05	0.05	0.080	0.008 7	0.25	8	6
Supercritical CO₂ group^[8]	320.00	0.04	0.04	0.030	0.015 0	0.25	10	8
Experiment type	Biot coefficient of matrix $α$	Biot coefficient of fracture $β$	Initial permeability of matrixi $k m 0$ /m²	Initial permeability of fracture $k f0$ /m²	Langmuir strain constant of matrix part $ε Lm$ /%	Langmuir strain constant of fracture part $ε Lf$ /%	Langmuir pressure constant of fracture part $p Lf$ /MPa	Langmuir pressure constant of matrix part $p Lm$ /MPa
N₂ group^[45]	0.30	0.60	3×10^－18	1×10^－17	0.4	0.2	4	2
CO₂ group^[44]	0.65	0.85	1×10^－19	1×10^－17	0.8	0.4	8	2
Supercritical CO₂ group^[8]	0.70	0.90	5×10^－19	1×10^－17	$0.4 (ε L 1 m) ， 1.2 (ε L 2 m)$	$0.4 (ε L f)$	$4 (p L f)$	$2 (p L 1 m) ， 8 (p L 2 m)$

Experiment type	Environment temperature T/K	Initial porosity of matrix $ϕ m0$	Initial porosity of fracture $ϕ f0$	Maximum adsorption capacity of matrix part $n max m$ /(m³·kg)	Maximum adsorption capacity of fracture part $n max f$ /(m³·kg)	Poisson's ratio $ν$	Young's modulus of matrix part $E m$ /GPa	Young's modulus of fracture part $E f$ /GPa
N₂ group^[45]	298.15	0.04	0.04	0.015	0.005 0	0.15	5	2
CO₂ group^[44]	313.00	0.05	0.05	0.080	0.008 7	0.25	8	6
Supercritical CO₂ group^[8]	320.00	0.04	0.04	0.030	0.015 0	0.25	10	8
Experiment type	Biot coefficient of matrix $α$	Biot coefficient of fracture $β$	Initial permeability of matrixi $k m 0$ /m²	Initial permeability of fracture $k f0$ /m²	Langmuir strain constant of matrix part $ε Lm$ /%	Langmuir strain constant of fracture part $ε Lf$ /%	Langmuir pressure constant of fracture part $p Lf$ /MPa	Langmuir pressure constant of matrix part $p Lm$ /MPa
N₂ group^[45]	0.30	0.60	3×10^－18	1×10^－17	0.4	0.2	4	2
CO₂ group^[44]	0.65	0.85	1×10^－19	1×10^－17	0.8	0.4	8	2
Supercritical CO₂ group^[8]	0.70	0.90	5×10^－19	1×10^－17	$0.4 (ε L 1 m) ， 1.2 (ε L 2 m)$	$0.4 (ε L f)$	$4 (p L f)$	$2 (p L 1 m) ， 8 (p L 2 m)$