Model for non-equilibrium evolution coal permeability of whole process under influences of gas pressure difference
-
摘要: 煤层气开采过程中基质与裂隙之间由于渗透性能的差异会长期存在气体压差,使煤岩处于非平衡态,非平衡态下的基质-裂隙相互作用是不断变化的。目前大多数研究对基质吸附变形与裂隙开度变化之间的关系没有充分的解释。为此,定义了一个新的与气体压差相关的非平衡因子,以描述基质吸附变形在不同时刻对裂隙开度变化与煤体变形的影响程度,建立了煤岩渗透率非平衡演化理论模型,该模型包括了从初始平衡到最终平衡状态的全过程。进一步地,采用试验数据对模型进行了验证,用数值模拟分析了气体抽注过程的气体压力时空分布规律与煤渗透率演化规律。结果表明:①基质孔隙压力变化滞后于裂隙压力变化,气体压差会长期存在,其大小先从0迅速增大后缓慢降低恢复至0;②渗透率演化可分为5个阶段,首尾两阶段表示平衡状态,中间3个阶段表示非平衡状态,非平衡演化过程中渗透率发生非单调变化;③基质块与基质桥的宽度之比会影响渗透率的变化幅度,用于煤储层的不同边界条件会影响渗透率的变化幅度和最终数值。研究结果对合理开采煤层气具有一定指导意义。Abstract: During the process of coalbed methane (CBM) extraction, gas pressure difference will exist for a long time due to different permeabilities between fracture and matrix, which makes coal in non-equilibrium state. The interaction between matrix and fracture in non-equilibrium state is constantly changing. At present, most studies do not fully explain the relationship between matrix adsorption deformation and change of fracture aperture. Therefore, a new non-equilibrium factor related to gas pressure difference is defined to describe the influence of matrix adsorption deformation on change of fracture aperture and coal deformation at different time. A theoretical model for non-equilibrium evolution of coal permeability is established, which includes the whole process from the initial equilibrium to the final one. Furthermore, the model is verified through the experimental data, and the spatial and temporal distribution of gas pressure and the evolution of coal permeability during gas injection and extraction are analyzed by the numerical simulation. The results show that: (1) The change of matrix pore pressure lags behind that of fracture pressure. The gas pressure difference exists for a long time, and its size increases rapidly from 0 and then slowly decreases to 0. (2) The permeability evolution can be divided into five stages. The first and the last stages represent the equilibrium state, the middle three stages represent the non-equilibrium state, and the permeability changes non-monotonically during the non-equilibrium evolution process. (3) The ratio of matrix width to matrix bridge width will affect the variation amplitude of permeability, and the different boundary conditions used in coal reservoir will affect the variation amplitude and the final value of permeability. The results are of certain guiding significance for CBM extraction reasonably.
-
0. 引言
煤层气作为一种高效且清洁的能源受到越来越多学者的关注[1]。在煤层气开采过程中会引起煤岩的力学变形和吸附解吸变形[2]。煤岩渗透率变化是这两种作用的竞争结果,裂隙与基质孔隙之间的气体压差将加剧这两种过程的复杂性,从而增加渗透率预测的难度。因此研究基质与裂隙内气体压差对正确认识储层内渗透率的动态演化具有重要意义。
常用的煤岩渗透率理论模型可分为单孔隙模型与双重孔隙模型两大类。由于煤裂隙系统的渗透率远大于基质系统的渗透率[2],早期许多研究中将煤岩视为单孔隙模型[3-5],不存在气体压差的概念。双重孔隙介质模型对煤岩孔隙结构有很好的描述,区分了裂隙压力和基质孔隙压力[6-9],然而对孔隙压力的理解仅停留在不同介质系统具有不同的孔隙压力上,忽略了气体压差对煤岩变形的影响。这些传统模型未能解释应力控制边界条件下的实验室测量结果,在匹配现场数据方面也仅有有限的成功[10]。
上述研究中使用的煤岩理想化模型有火柴杆模型与立方体模型,模型中假设裂隙将基质完全隔开。而在实际情况中裂隙往往不能贯穿各基质块,基质之间通过基质桥连接,因此需要考虑基质与裂隙的相互作用[11]。这种相互作用使得基质吸附变形只有一部分影响裂隙开度变化,另一部分则影响煤整体变形。为了解决这一问题,一些研究通过引入一个值为0~1之间的内部膨胀系数f来表征基质吸附变形对裂隙开度变化的贡献[11-14],并通过数据拟合的方式将系数f确定为一个常量。然而,由于煤岩基质孔隙的致密性,抽注时裂隙压力与基质孔隙压力无法瞬间达到平衡。吸附性气体在裂隙与基质之间的长期扩散过程中,煤岩经历初始平衡、非平衡和最终平衡3个状态[15-16],在低压条件下甚至永远无法达到最终平衡[17]。非平衡态演化过程中基质吸附变形对裂隙开度变化的贡献是随时间变化的,演化初期裂隙附近的基质非均匀变形(局部变形)几乎完全影响裂隙开度变化(f≈1),而在最终平衡阶段,基质均匀变形(整体变形)几乎不影响裂隙开度变化(f≈0)。因此f被看作一个常量是不合理的。
过去渗透率测量的实验室观察均在短期内完成,得到的试验数据无法反映煤岩的长期非平衡演化过程。Wei等[18]对煤岩试样进行了超过80 d的CO2注气试验,使用改进的脉冲衰减技术进行了渗透率连续测量,研究了煤岩从初始平衡到最终平衡的裂隙渗透率演化。试验结果显示,煤岩渗透率随时间先迅速增大后降低,最后再缓慢恢复。为了解释这一结果,近期许多研究基于非平衡理论进行了建模。Wei等[19-20]通过区分裂隙应变的方式描述了局部变形与整体变形之间的关系,又通过划分基质区域的方式区分了不同位置处的基质变形对裂隙开度变化和煤体变形的影响。Zeng等[21-22]引入气体侵入/排出体积与基质整体体积的比值来量化了局部膨胀/收缩到全局膨胀/收缩的过程。Huang等[23]假设最靠近裂隙的基质点压力与裂隙压力相同,使用该点膨胀应变与整体膨胀应变的比值来表示非平衡因子,以描述基质的非均匀变形。Liu等[24]根据通过薄膜的非稳态扩散解析解计算了基质中气体浓度的分布,定义了新的非平衡因子。这些研究都推出了更符合长期试验数据的理论成果。
上述关于非平衡的研究在确定空间或时间变化因子方面存在很高的难度,或者对基质裂隙相互作用的力学解释不够充分。裂隙与基质的气体压差除了造成煤岩长期处于非平衡的吸附状态之外,还会造成基质和裂隙系统内有效应力的变化[25]。因此,本文考虑了气体压差引起的有效应力变化,并定义了一个新的与气体压差相关的非平衡因子,以描述非平衡演化过程中基质吸附变形对裂隙开度变化与煤体变形的影响。同时,本文考虑了基质桥的吸附变形,对基质块与基质桥吸附变形对裂隙开度变化的共同作用做出了解释。在此基础上,建立了一个气体压差影响下双重孔隙煤岩非平衡演化过程的渗透率理论模型,并进行了该理论模型控制下的煤岩试样注气与原位煤储层采气的数值模拟。研究结果对煤层气开采工程具有一定理论指导作用。
1. 煤岩变形控制方程
一般情况下,气体在基质中的储存形式包含游离态和吸附态,而裂隙中只有游离态。即气体吸附仅发生在基质系统中,不发生在裂隙系统中。假设吸附引起的应变在各个方向的分量相等,双重孔隙煤岩变形弹性力学本构关系可以定义为[6]
εij=12Gbσij−(16Gb−19Kb)σkkδij+α3Kbpmδij+β3Kbpfδij+εbs3δij。 (1) 式中:下标m,f,b分别代表基质、裂隙与煤体;G为剪切模量,G=E/2(1+ν),E为杨氏弹性模量(MPa);ν为泊松比;K为体积模量,K=E/3(1-2ν);α与β分别为基质与裂隙的Biot系数,表达式分别为α=1-Km/ Ks与β=1-Kb/Km;Ks为煤基质骨架(煤颗粒)的体积模量;p为孔隙压力(MPa);总应力σkk=σ11+σ22+σ33,MPa;δij为Kronecker符号;εbs为基质吸附变形产生的煤体应变。根据Langmuir等温吸附方程,基质的吸附应变εms可以定义为
εms=εLmpmpm+PLm。 (2) 式中:εLm与PLm分别为煤基质的Langmuir应变常数与压力常数。
将式(1)与弹性力学平衡方程、几何方程联立,可得基质与裂隙系统的Navier型变形控制方程:
G_{\mathrm{b}} \boldsymbol{u}_{i, k k}+\frac{G_{\mathrm{b}}}{1-2 v} \boldsymbol{u}_{k, k i}-\alpha P_{\mathrm{m}, i}-\beta P_{\mathrm{f}, i}-K_{\mathrm{b}} \varepsilon_{\mathrm{bs}, i}+\boldsymbol{F}_i=0 。 (3) 式中:u为位移矢量,F为体积力。
2. 煤岩孔隙率与渗透率模型
2.1 非平衡因子定义
向煤岩注气后,气体侵入裂隙壁附近的基质区域,引起基质局部膨胀,进而压缩裂隙,渗透率减小。随着气体向基质内部扩散,膨胀变形逐渐均匀,局部变形对裂隙开度的影响减弱,渗透率逐渐恢复。采气过程则与之相反。因此,引入一个数值介于0~1的非平衡因子Ue来量化基质吸附变形对裂隙开度变化的贡献程度。局部应变是由气体压差引起的基体与裂隙之间的应变差产生的[19],因此将该系数与气体压差建立起联系。
U_{\mathrm{e}}=\frac{\left|p_{\mathrm{f}}-p_{\mathrm{m}}\right|}{\max \left\{p_{\mathrm{f}}, p_{\mathrm{m}}\right\}}= \begin{cases}\frac{p_{\mathrm{f}}-p_{\mathrm{m}}}{p_{\mathrm{f}}} & \left(p_{\mathrm{f}}>p_{\mathrm{m}}\right) \\ \frac{p_{\mathrm{m}}-p_{\mathrm{f}}}{p_{\mathrm{m}}} & \left(p_{\mathrm{m}}>p_{\mathrm{f}}\right)\end{cases}。 (4) Ue从0增大至1后再减小最终回到0,当Ue=0时,基质的吸附变形不会对裂隙开度变化产生影响,当Ue=1时,基质的吸附变形完全影响裂隙开度变化。 {p_{\text{f}}} > {p_{\text{m}}} 适用于注气情况; {p_{\text{m}}} > {p_{\text{f}}} 适用于采气情况。
2.2 孔隙率与渗透率模型建立
煤裂隙是气体的主要流动通道,而煤基质则为气体提供主要的储存空间,裂隙渗透率远大于基质渗透率[6-9],因此,此处只研究裂隙渗透率模型。煤裂隙孔隙率的表达式为
\phi_{\mathrm{f}}=\frac{V_{\mathrm{f}}}{V_{\mathrm{b}}}。 (5) 式中:Vf为裂隙体积;Vb为煤体体积。对式(5)两边取微分,可得裂隙孔隙率的改变量为
{\text{d}}{\phi _{\text{f}}} = {\text{d}}\left( {\frac{{{V_{\text{f}}}}}{{{V_{\text{b}}}}}} \right) = {\phi _{\text{f}}}\left( {\frac{{{\text{d}}{V_{\text{f}}}}}{{{V_{\text{f}}}}} - \frac{{{\text{d}}{V_{\text{b}}}}}{{{V_{\text{b}}}}}} \right)。 (6) 根据式(1)可得裂隙体应变增量与煤岩体应变增量分别为[26]
\frac{\text{d}{V}_{\text{b}}}{{V}_{\text{b}}}=-\frac{1}{{K}_{\text{b}}}\left[\text{d}\overline{\sigma }-\alpha \text{d}{p}_{\text{m}}-\beta \text{d}{p}_{\text{f}}\right]+\text{d}{\varepsilon }_{\text{bs}}\text{ }。 (7) 式中:\overline \sigma 为平均应力,\overline \sigma = - {\sigma _{kk}}/3。煤体、裂隙的变形均由两部分组成:一部分由有效应力变化引起,另一部分由气体吸附引起[27],类似的,裂隙应变增量则可表示为[25]
\frac{\text{d}{V}_{\text{f}}}{{V}_{\text{f}}}=-\frac{1}{{K}_{\text{f}}}\left[\text{d}\overline{\sigma }-\beta \text{d}{p}_{\text{f}}-\text{d}p(t)\right]+\text{d}{\varepsilon }_{\text{fs}}\text{ }。 (8) 式中:p(t)为裂隙与基质之间的压差表达式,其值随注气或采气过程的时间而变化,p(t)= {p_{\text{f}}} - {p_{\text{m}}} ; {\varepsilon _{{\text{fs}}}} 为由基质吸附变形与基质桥吸附变形共同产生的裂隙应变。方程(7),(8)等式右侧第一项为有效应力变化产生的应变增量,第二项为吸附应变增量。
将式(7),(8)代入式(6)可得
\begin{aligned} \frac{\mathrm{d} \phi_{\mathrm{f}}}{\phi_{\mathrm{f}}}=\left(\frac{1}{K_{\mathrm{b}}}-\right. & \left.\frac{1}{K_{\mathrm{f}}}\right)\left(\mathrm{d} \bar{\sigma}-\beta \mathrm{d} p_{\mathrm{f}}\right)- \\ & \frac{\alpha}{K_{\mathrm{b}}} \mathrm{~d} p_{\mathrm{m}}+\frac{1}{K_{\mathrm{f}}} \mathrm{~d} p(t)+\mathrm{d} \varepsilon_{\mathrm{fs}}-\mathrm{d} \varepsilon_{\mathrm{bs}} \end{aligned}。 (9) 根据Betti-Maxwell互易定理可得Kf为[28]
{K}_{\text{f}}=\frac{{\varphi }_{\text{f}}}{\beta }{K}_{\text{b}}\text{ }。 (10) 考虑基质-裂隙相互作用时,常用的煤岩理想化模型为基质桥模型[11]。图 1为基质桥模型代表性单元,向煤岩注气时,基质吸附变形与基质桥吸附变形共同影响裂隙开度变化。图 1中深色区域表示吸附膨胀变形前的体积,浅色区域为吸附膨胀变形产生的体积差,红色箭头表示基质或基质桥的膨胀应变方向。基质桥与裂隙之间的接触比表面积较大,且基质桥自身体积较小,因此可以假设在注气后裂隙与基质桥之间的传质速度很快,即基质桥内的气体压力能在很短的时间内与裂隙压力达到平衡。基质由于体积较大,在与裂隙的传质过程中存在长期的非平衡态演化。因此,推导中假设基质桥中的气体吸附压力为与裂隙压力相等,基质中的气体压力则与裂隙压力存在压差。
基质桥是由于裂隙没有完全贯穿基质而产生的连接各基质块的部分,因此可以认为基质桥属于基质的一部分,其吸附特性与基质相同,可采用相同的Langmuir应变常数与压力常数。因此,吸附导致的基质桥膨胀应变为
{\varepsilon }_{\text{ms}}^{\text{r}}={\varepsilon }_{\text{Lm}}\frac{{p}_{\text{f}}}{{p}_{\text{f}}+{P}_{\text{Lm}}}\text{ }。 (11) 在注气过程中,基质桥吸附膨胀会增大裂隙开度,与基质吸附膨胀而压缩裂隙起到相反的作用。基于各向同性假设,由吸附膨胀导致的各个方向线应变为{\varepsilon _{{\text{s}}x}} = {\varepsilon _{{\text{s}}y}} = {\varepsilon _{{\text{s}}z}} = {\varepsilon _{\text{s}}}/3,以裂隙压缩方向为负,在水平方向,由基质吸附变形与基质桥吸附变形相互作用导致的裂隙开度变化为
\Delta {b}_{\text{s}}=b\cdot \frac{\Delta {\varepsilon }_{\text{ms}}^{\text{r}}}{3}-{U}_{\text{e}}\cdot a\cdot \frac{\Delta {\varepsilon }_{\text{ms}}}{3}\text{ }。 (12) 式中:b为裂隙开度(m);a为基质宽度(m)。水平方向上吸附导致的裂隙线应变增量为
\Delta {\varepsilon }_{\text{fs}}^{x}\text{=}\frac{\Delta {b}_{\text{s}}}{b}=\frac{\Delta {\varepsilon }_{\text{ms}}^{\text{r}}}{3}-{U}_{\text{e}}\cdot \frac{a}{b}\cdot \frac{\Delta {\varepsilon }_{\text{ms}}}{3}\text{ }。 (13) 因此,基质和基质桥吸附变形的共同影响的裂隙应变增量为
\Delta {\varepsilon }_{\text{fs}}\text{=}3{\varepsilon }_{\text{fs}}^{x}=\Delta {\varepsilon }_{\text{ms}}^{\text{r}}-{U}_{\text{e}}\cdot \frac{a}{b}\cdot \Delta {\varepsilon }_{\text{ms}}\text{ }。 (14) 由于基质吸附变形一部分导致裂隙开度变化,另一部分导致煤体变形,因此同理式(14),可得基质吸附变形导致的煤体应变增量为
\Delta {\varepsilon _{{\text{bs}}}} = {\text{(}}1 - {U_{\text{e}}}{\text{)}} \cdot \frac{a}{{a + b}} \cdot {\varepsilon _{{\text{ms}}}}。 (15) 由于裂隙开度b远小于基质宽度a,即a±b≈a,式(15)可简化为
\Delta {\varepsilon }_{\text{bs}}=(1-{U}_{\text{e}})\cdot {\varepsilon }_{\text{ms}}\text{ }。 (16) 联立式(2),(11),(14),(16)代入式(9),并对等式两侧进行积分、简化,得裂隙孔隙度模型为
\begin{gathered} \frac{\phi_{\mathrm{f}}}{\phi_{\mathrm{f} 0}}=\exp \left\{\left(\frac{1}{K_{\mathrm{b}}}-\frac{1}{K_{\mathrm{f}}}\right)\left[\left(\bar{\sigma}-\bar{\sigma}_0\right)-\beta\left(p_{\mathrm{f}}-p_{\mathrm{f} 0}\right)\right]-\frac{\alpha}{K_{\mathrm{b}}}\left(p_{\mathrm{m}}-p_{\mathrm{m} 0}\right)+\right. \\ \frac{1}{K_{\mathrm{f}}}\left[p(t)-p\left(t_0\right)\right]+\varepsilon_{\mathrm{Lm}} P_{\mathrm{Lm}} \cdot\left[\frac{p_{\mathrm{f}}-p_{\mathrm{f} 0}}{\left(p_{\mathrm{f}}+P_{\mathrm{Lm}}\right)\left(p_{\mathrm{f} 0}+P_{\mathrm{Lm}}\right)}-\right. \\ \left.\left.\quad\left(U_{\mathrm{e}} \cdot \frac{a}{b}+1\right) \cdot \frac{p_{\mathrm{m}}-p_{\mathrm{m} 0}}{\left(p_{\mathrm{m}}+P_{\mathrm{Lm}}\right)\left(p_{\mathrm{m} 0}+P_{\mathrm{Lm}}\right)}\right]\right\}。 \end{gathered} (17) 根据渗透率与孔隙率的立方关系
\frac{{k}_{\text{f}}}{{k}_{\text{f}0}}={\left(\frac{{\phi }_{\text{f}}}{{\phi }_{\text{f}0}}\right)}^{3}\text{ }。 (18) 将式(17)代入式(18),得裂隙渗透率比值为
\begin{aligned} \frac{k_{\mathrm{f}}}{k_{\mathrm{f} 0}}= & \exp \left\{\left(\frac{3}{K_{\mathrm{b}}}-\frac{3}{K_{\mathrm{f}}}\right)\left[\left(\bar{\sigma}-\bar{\sigma}_0\right)-\beta\left(p_{\mathrm{f}}-p_{\mathrm{f} 0}\right)\right]-\frac{3 \alpha}{K_{\mathrm{b}}}\left(p_{\mathrm{m}}-p_{\mathrm{m} 0}\right)+\right. \\ & \frac{3}{K_{\mathrm{f}}}\left[p(t)-p\left(t_0\right)\right]+3 \varepsilon_{\mathrm{Lm}} P_{\mathrm{Lm}} \cdot\left[\frac{p_{\mathrm{f}}-p_{\mathrm{f} 0}}{\left(p_{\mathrm{f}}+P_{\mathrm{Lm}}\right)\left(p_{\mathrm{f} 0}+P_{\mathrm{Lm}}\right)}-\right. \\ & \left.\left.\left(U_{\mathrm{e}} \cdot \frac{a}{b}+1\right) \cdot \frac{p_{\mathrm{m}}-p_{\mathrm{m} 0}}{\left(p_{\mathrm{m}}+P_{\mathrm{Lm}}\right)\left(p_{\mathrm{m} 0}+P_{\mathrm{Lm}}\right)}\right]\right\}。 \end{aligned} (19) 式(17),(19)即为本文裂隙孔隙率与渗透率模型的一般表达式。
2.3 特定边界条件下的渗透率模型
在实验室中,通常对煤岩试样应用恒围压边界条件[10, 18],即d\overline \sigma =0。根据式(19),并由Kb \gg Kf,可得恒围压条件下的渗透率模型为
\begin{gathered} \frac{k_{\mathrm{f}}}{k_{\mathrm{f} 0}}=\exp \left\{\frac{3 \beta}{K_{\mathrm{f}}}\left(p_{\mathrm{f}}-p_{\mathrm{f} 0}\right)-\frac{3 \alpha}{K_{\mathrm{b}}}\left(p_{\mathrm{m}}-p_{\mathrm{m} 0}\right)+\frac{3}{K_{\mathrm{f}}}[p(t)-\right. \\ \\ \left.p\left(t_0\right)\right]+3 \varepsilon_{\mathrm{Lm}} P_{\mathrm{Lm}} \cdot\left[\frac{p_{\mathrm{f}}-p_{\mathrm{f} 0}}{\left(p_{\mathrm{f}}+P_{\mathrm{Lm}}\right)\left(p_{\mathrm{f} 0}+P_{\mathrm{Lm}}\right)}-\right. \\ \left.\left.\quad\left(U_{\mathrm{e}} \cdot \frac{a}{b}+1\right) \cdot \frac{p_{\mathrm{m}}-p_{\mathrm{m} 0}}{\left(p_{\mathrm{m}}+P_{\mathrm{Lm}}\right)\left(p_{\mathrm{m} 0}+P_{\mathrm{Lm}}\right)}\right]\right\} 。 \end{gathered} 。 (20) 单轴应变条件通常被作为原位煤层的边界条件,在煤层气开采研究中被广泛应用[12, 26]。结合边界条件( {\varepsilon _{11}} = {\varepsilon _{22}} = 0 , \Delta {\sigma _{33}} = 0 ),根据式(1)可得
\begin{aligned} \Delta \sigma_{11}=\Delta \sigma_{22} & =\frac{v}{1-v} \Delta \sigma_{33}+ \\ & \frac{2 v-1}{1-v}\left(\alpha \Delta p_{\mathrm{m}}+\beta \Delta p_{\mathrm{f}}\right)-\frac{E_{\mathrm{b}}}{3(1-v)} \Delta \varepsilon_{\mathrm{bs}} \end{aligned} 。 (21) 根据平均应力的表达式,可得平均应力增量为
\overline{\sigma }-{\overline{\sigma }}_{0}\text{=}-\frac{1}{3}({\sigma }_{11}-{\sigma }_{110}+{\sigma }_{22}-{\sigma }_{220}+{\sigma }_{33}-{\sigma }_{330})\text{ }。 (22) 将式(16),(21)代入式(22)可得
\begin{array}{r} \bar{\sigma}-\bar{\sigma}_0=\frac{2-4 v}{3(1-v)}\left[\alpha\left(p_{\mathrm{m}}-p_{\mathrm{m} 0}\right)+\beta\left(p_{\mathrm{f}}-p_{\mathrm{f} 0}\right)\right]+ \\ \frac{2 E_{\mathrm{b}}}{9(1-v)}\left(1-U_{\mathrm{e}}\right)\left(\varepsilon_{\mathrm{ms}}-\varepsilon_{\mathrm{ms} 0}\right) \end{array}。 (23) 将式(23)代入式(19),并由Kb \gg Kf,可得单轴应变条件下的裂隙渗透率模型为
\begin{gathered} \frac{{{k_{\text{f}}}}}{{{k_{{\text{f}}0}}}} = \exp \left\{ {\frac{1}{{{K_{\text{f}}}}}\left[ {\frac{{1 + \nu }}{{1 - \nu }}\beta ({p_{\text{f}}} - {p_{{\text{f0}}}}) + 3(p(t) - p({t_0}))} \right] - \frac{\alpha }{{{K_{\text{b}}}}}\frac{{1 + \nu }}{{1 - \nu }} \cdot } \right. \hfill \\ {\text{ }}({p_{\text{m}}} - {p_{{\text{m}}0}}) + 3{\varepsilon _{{\text{Lm}}}}{P_{{\text{Lm}}}} \cdot \left[ {\frac{{{p_{\text{f}}} - {p_{{\text{f}}0}}}}{{({p_{\text{f}}} + {P_{{\text{Lm}}}})({p_{{\text{f}}0}} + {P_{{\text{Lm}}}})}}} \right. - \hfill \\ \left. {\left( {\frac{{2{E_{\text{b}}}(1 - {U_{\text{e}}})}}{{9{K_{\text{f}}}(1 - \nu )}} + \frac{a}{b}{U_{\text{e}}} + 1} \right) \cdot \frac{{{p_{\text{m}}} - {p_{{\text{m}}0}}}}{{({p_{\text{m}}} + {P_{{\text{Lm}}}})({p_{{\text{m0}}}} + {P_{{\text{Lm}}}})}}} \right\} \end{gathered}。 (24) 3. 煤岩多物理场耦合
3.1 气体运移方程
气体在基质中的运移遵循Fick定律,在裂隙中的运移遵循Darcy定律。结合气体运移质量守恒定律,气体在基质与裂隙中的运移方程为[6]
\frac{{\partial {c_{\text{m}}}}}{{\partial t}} + \nabla \cdot ( - D\nabla {c_{\text{m}}}) = - {Q_{\text{s}}}\text{,} (25) \frac{\partial {m}_{\text{f}}}{\partial t}+\nabla \cdot \left({\rho }_{\text{g}}\overrightarrow{{v}_{\text{g}}}\right)={Q}_{\text{s}}\text{ }。 (26) 式中:cm为基质中的气体体积浓度(kg/m3);mf为裂隙中的气体含量(kg/m3);D为气体扩散系数(m2/s);ρg为气体密度(kg/m3);\overrightarrow {{v_{\text{g}}}} 为达西流速矢量(m/s);Qs为基质与裂隙之间的气体交换速率(kg/(m3·s))。
基质中包含游离态和吸附态气体,气体吸附量可由Langmuir方程计算;而裂隙中只有游离态气体,所以在基质与裂隙中气体含量表达式分别为
{c}_{\text{m}}={\rho }_{\text{gm}}{\phi }_{\text{m}}+{\rho }_{\text{c}}{\rho }_{\text{ga}}\frac{{p}_{\text{m}}{V}_{\text{Lm}}}{{p}_{\text{m}}+{P}_{\text{Lm}}}\text{ }\text{,} (27) {m}_{\text{f}}={\rho }_{\text{gf}}{\phi }_{\text{f}}\text{ }。 (28) 式中:ρc为煤岩的密度(kg/m3);ρga为标准状况下的气体密度(kg/m3);VLm为Langmuir体积常数(m3/kg)。
根据理想气体状态方程可得气体密度表达式为
{\rho }_{\text{g}}=\frac{M}{RT}p\text{ }。 (29) 式中:M为气体的摩尔质量(kg/mol);R为摩尔气体常数(J/(mol·K));T为温度(K)。
基质与裂隙之间的气体交换速率表达式为
{Q_{\text{s}}} = D \cdot w({\rho _{{\text{gm}}}} - {\rho _{{\text{gf}}}}) = - \frac{{DwM}}{{RT}}p{\text{(}}t{\text{)}}。 (30) 式中:w为形状因子(m-2)。
裂隙中气体达西流速矢量表达式为
\overrightarrow{{v}_{\text{g}}}=-\frac{{k}_{\text{f}}}{\mu }\nabla {p}_{\text{f}}\text{ }。 (31) 式中:μ为气体的动力学黏度(Pa·s)。
联立式(25)~(31),可得煤基质与裂隙中的气体运移耦合控制方程
\begin{aligned} & {\left[\phi_{\mathrm{m}}+\frac{\rho_{\mathrm{c}} p_{\mathrm{a}} V_{\mathrm{Lm}} P_{\mathrm{Lm}}}{\left(p_{\mathrm{m}}+P_{\mathrm{Lm}}\right)^2}\right] \frac{\partial p_{\mathrm{m}}}{\partial t}+p_{\mathrm{m}} \frac{\partial \phi_{\mathrm{m}}}{\partial t}+} \\ & \nabla \cdot\left[-D\left[\phi_{\mathrm{m}}+\frac{\rho_{\mathrm{c}} p_{\mathrm{a}} V_{\mathrm{Lm}} P_{\mathrm{Lm}}}{\left(p_{\mathrm{m}}+P_{\mathrm{Lm}}\right)^2}\right] \nabla p_{\mathrm{m}}\right]=D w \Delta p(t), \end{aligned} (32) {\phi }_{\text{f}}\frac{\partial {p}_{\text{f}}}{\partial t}+{p}_{\text{f}}\frac{\partial {\phi }_{\text{f}}}{\partial t}+\nabla \cdot \left(-{p}_{\text{f}}\frac{{k}_{\text{f}}}{\mu }\nabla {p}_{\text{f}}\right)=-Dw\Delta p(t)\text{ }。 (33) 3.2 流固耦合关系
煤层气开采过程中气体流动与扩散会导致气体压力变化,进而使煤岩发生变形,而煤岩变形会引起其孔隙率与渗透率的改变,又反过来影响气体流动与扩散。这些物理场的交叉耦合关系如图 2所示。涉及多物理场耦合的煤岩渗透率演化过程可于COMSOL Multiphysics中进行数值模拟求解。
4. 非平衡模型验证
在过去关于基质与裂隙相互作用的研究中,建立的渗透率模型往往能够顺利匹配“渗透率-压力”试验数据[12-14],但无法成功匹配“渗透率-时间”试验数据。这是因为这些模型通过拟合出一个定值f来描述基质与裂隙相互作用,而实际上基质吸附变形对裂隙开度变化的贡献程度与时间有关,定值f无法描述非平衡演化过程。本部分使用COMSOL模拟煤试样注气过程,以得到式(20)控制下的渗透率演化预测结果。将预测曲线与Wei等[18]试验获取的长期渗透率数据进行对比,以验证本文渗透率模型的合理性。该试验在恒围压条件下进行,根据岩芯形状与其受力的几何对称性,在数值模拟中可将圆柱形煤试样简化为二维轴对称模型,如图 3所示。几何模型底部边界为辊支撑,其余边界设置6 MPa恒定压力,上边界设置3 MPa的注入压力,煤试样内初始气体压力为0.1 MPa。
煤体弹性模量E取3 GPa,煤基质弹性模型Em取10 GPa[20],泊松比 \nu 取0.34,则可分别计算得到煤体体积模量Kb与煤基质体积模量Km;煤裂隙体积模量Kf由式(10)计算;煤基质骨架(颗粒)体积模量Ks可使用下式计算[31]:
{K_{\text{s}}}{\text{ = }}\frac{{{K_{\text{m}}}}}{{1 - 3{\phi _{{\text{m0}}}}(1 - \nu )/\left[ {2(1 - 2\nu )} \right]}} 。 (34) 根据各类体积模量可计算得到Biot系数α与β的值。其他相关参数取自Wei等[20]的研究,将模拟所需的所有参数汇总于表 1。
表 1 用于理论模型验证的参数Table 1. Parameters used for verification of theoretical model变量 值 物理意义 E 3000 MPa 煤岩弹性模量 Em 10000 MPa[20] 煤基质弹性模量 \nu 0.34 泊松比 ρc 1400 kg/m³ 煤岩密度 ρga 1.98 kg/m³ 标况下CO2密度 M 0.044 kg/mol CO2的摩尔质量 μ 1.84×10-5 Pa·s CO2动力黏度 R 8.414 J/(mol·K) 摩尔气体常数 T 293 K 温度 pa 0.1 MPa 大气压 pinj 3 MPa[20] 注入压力 σ0 6 MPa[20] 围压 D 2×10-8 m²/s[20] CO2扩散系数 w 12 m-2[20] 形状因子 εLm 0.0518[20] Langmuir应变常数 VLm 0.017 m3/kg[20] Langmuir体积常数 PLm 4.48 MPa[20] Langmuir压力常数 φf0 1.2%[20] 初始裂隙孔隙率 φm0 2%[20] 初始基质孔隙率 kf0 1×10-17 m2[20] 初始裂隙渗透率 模型验证结果如下图 4所示,并与Lu等[12]建立的模型进行比较。Lu等[12]模型考虑了基质裂隙相互作用,f取值为1和0时可以表示渗透率演化的上下限。f=1时,基质吸附变形完全影响裂隙开度的变化,此时的预测曲线能够成功匹配前期的试验数据,这是因为在前期基质发生局部变形,主要影响裂隙开度变化。f=0时,基质吸附变形对裂隙开度变化没有影响,此时的预测曲线则在后期接近试验数据,这是因为后期逐渐接近最终平衡状态,煤岩发生整体变形,对裂隙渗透率影响不大。但Lu等[12]模型无法描述从初始平衡至最终平衡演化的过程。相比之下,本文模型可以很好地匹配整个试验过程的渗透率数据,有效解释了非平衡演化过程的渗透率趋势。
5. 非平衡模型应用
5.1 模型描述
煤层产气过程与气体注入煤样过程相反,但渗透率的演化机制相同[19]。因此,本文理论模型也可用于解释煤层气开采过程中渗透率的演化。本部分采用COMSOL研究式(24)控制下的煤层气开采过程的渗透率演化规律。如图 5所示,建立一个200 m×200 m的煤层几何模型,中心位置为坐标原点。模型下边界设置为固定约束,左右边界设置为辊支撑,上边界设置竖直向下的恒定边界荷载Fy,以此来模拟单轴应变条件。
煤层初始状态下的基质孔隙压力与裂隙压力相等,均设为4 MPa[24]。在中心位置处钻一半径为0.5 m的井口,井口边界压力为0.3 MPa[24]。取A点(50 m,50 m)以研究采气过程中气体压力与渗透率的演化规律。煤层气开采数值模拟中所用的参数列于表 2,其他参数与表 1所列相同。
表 2 用于现场开采煤层气数值模拟的参数Table 2. Parameters used for numerical simulation of CBM5.2 气体压力时空分布规律
图 6绘制了气体压力随时间的分布规律。开采前,A点的基质孔隙压力等于裂隙压力,均为初始值4 MPa。钻井后进行开采时,由于裂隙渗透率远大于基质渗透率,裂隙压力首先迅速降低,初始平衡状态被打破,煤储层开始非平衡态演化。而基质孔隙压力在103 h时才出现较为明显的下降趋势,且在后续的很长时间的采气过程中均大于裂隙压力。相应地,非平衡因子Ue随着气体压差的变化从0迅速增大后缓慢降低。在106 h之后才恢复至0,煤储层到达最终平衡状态。非平衡状态的演化持续了很长的时间。
图 7描述了与井点同一水平位置处(x轴处)的气体压力空间分布规律。对比图 7(a),(b)可知,裂隙压力与基质孔隙压力均从远离井点的位置向井点处不断降低,横坐标x=0处的气体压力为井底压力0.3 MPa。但裂隙压力随空间的降低幅度大于基质孔隙压力的降低幅度。这是因为裂隙中的气体不断从远离井点处向近井点处运输,近井点处的气体更快排出井外,因此近井点处的裂隙压力较低。而基质孔隙中的气体则主要是向裂隙中传质,传质过程发生在储层的各个位置,因此基质孔隙压力曲线趋于水平。图 6显示在非平衡演化过程中,储层A点处的裂隙压力在任意时间均小于基质孔隙压力。图 7则显示在非平衡演化过程中,裂隙压力在同一时间的任意位置均小于基质孔隙压力。
5.3 裂隙渗透率演化规律
为了分析采气过程中的裂隙渗透率演化规律,绘制了其随时间的演化曲线,由图 8可知,渗透率演化同样可分为5个阶段,经历了先减小后增大再减小的趋势。其演化原理与注气过程相同,但过程相反。阶段Ⅰ代表初始平衡状态,裂隙压力与基质孔隙压力均为初始值。阶段Ⅱ中,采气使裂隙压力迅速减小,裂隙渗透率因有效应力的增大而减小。同时基质桥产生解吸收缩会对裂隙两侧的基质产生一种类似“牵拉”的力,同样会使裂隙开度减小。在此阶段,初始平衡状态被打破,非平衡演化过程开始。阶段Ⅲ中,由于裂隙与基质之间的气体压差较大,气体解吸发生在靠近裂隙的基质区域,产生的局部收缩使裂隙开度增大,渗透率增大。阶段Ⅳ中,随着基质内部气体的排出,裂隙与基质的气体压差减小,基质局部收缩向整体收缩演化,基质吸附变形从主要增大裂隙开度变为主要减小煤体体积,渗透率逐渐降低。阶段Ⅴ中,基质孔隙压力与裂隙压力再次相等,煤岩进入最终平衡状态。因此,在开采前期,有效应力与.基质桥的收缩对渗透率的影响占主导作用;随着开采的进行,主导渗透率演化的因素变为基质局部收缩与整体收缩效应。
图 9显示了裂隙渗透率随孔隙压力的变化规律。采气过程为降压过程,孔隙压力从4 MPa降至0.3 MPa的整个过程中,渗透率比值先降低后升高再降低,其随裂隙压力和基质孔隙压力变化的趋势相同,但变化位置不同。渗透率比值从1降至最小值时,基质孔隙压力几乎不变,而裂隙压力降低至1.5 MPa左右。渗透率比值从最小值增大到最大值时,基质孔隙压力降低至0.75 MPa左右,裂隙压力降低至0.4 MPa左右。渗透率比值从最大值降至最终平衡状态下的比值时,基质孔隙压力与裂隙压力从不同的压力值分别降低至0.3 MPa。由图 9也可看出,基质孔隙压力的降低始终滞后于裂隙压力,且其二者共同决定渗透率的演化过程。值得注意的是,根据图 6的注气曲线可得,最终平衡状态下的渗透率因裂隙压力的升高而大于初始平衡状态下的渗透率,因此图 8,9采气曲线的最终平衡状态下的渗透率应该由于裂隙压力的降低而小于初始平衡状态下的渗透率,但曲线却呈现出了相反的结果,这可能是由适用于现场的单轴应变条件导致。图 10绘制了不同条件下的渗透率演化规律,针对这一现象做出了解释。
由图 10可得,恒围压条件下的和单轴应变条件下的渗透率演化趋势相同,但变化幅度不同,且最终平衡状态下的渗透率大小也不同。在恒围压条件下,最终平衡状态的渗透率小于初始平衡状态的渗透率。而单轴应变条件下,最终平衡状态的渗透率大于初始平衡状态的渗透率,这是因为单轴应变条件限制了x与y方向的应变,使得基质收缩更多的影响了裂隙开度变化而不是煤整体变形,因此渗透率变化幅度更大。此外,对比单轴应变条件下不同基质宽度与基质桥宽度之比(a/b)的3条曲线可得,a/b值不会影响最终平衡状态下的渗透率,但会影响渗透率的变化幅度。a/b越大,基质占比越大,基质桥占比越小,基质吸附变形带来的影响越大,基质桥对基质变形的阻碍作用越小,渗透率变化的幅度越大。
6. 结论
建立了一个考虑裂隙与基质气体压差和基质桥吸附变形的煤岩非平衡演化全过程渗透率模型;模型中定义了一个与气体压差相关的非平衡因子,以表征煤岩抽注过程中基质吸附变形在不同时刻对裂隙开度变化的影响程度;使用渗透率长期演化的试验数据对该模型进行了验证;并采用数值模拟进一步分析了原位开采煤层气时的气体压力时空分布与渗透率演化规律,具体得到以下3点结论。
(1)煤岩抽注过程中基质孔隙压力的变化总是滞后于裂隙压力的变化,裂隙与基质之间则会产生气体压差,进而造成煤结构的变形。这种变形既包括因有效应力改变而产生的力学变形,也包括非平衡状态下的基质非均匀吸附变形。
(2)煤试样注气和煤层采气过程的渗透率长期演化的原理相同,过程相反。其演化可分为5个阶段,首尾两阶段代表初始平衡与最终平衡状态,中间3个阶段代表非平衡状态,非平衡演化过程中煤岩渗透率受基质非均匀吸附变形的影响而发生非单调变化。
(3)用于煤储层的不同边界条件会影响渗透率演化曲线的变化幅度和最终数值,单轴应变条件下的渗透率变化幅度和最终数值均大于恒定围压条件下的渗透率变化幅度和最终数值。不同基质块与基质桥宽度之比会影响渗透率曲线的变化幅度,但不影响其最终数值。基质块与基质桥宽度之比越大,渗透率变化幅度增大。
-
表 1 用于理论模型验证的参数
Table 1 Parameters used for verification of theoretical model
变量 值 物理意义 E 3000 MPa 煤岩弹性模量 Em 10000 MPa[20] 煤基质弹性模量 \nu 0.34 泊松比 ρc 1400 kg/m³ 煤岩密度 ρga 1.98 kg/m³ 标况下CO2密度 M 0.044 kg/mol CO2的摩尔质量 μ 1.84×10-5 Pa·s CO2动力黏度 R 8.414 J/(mol·K) 摩尔气体常数 T 293 K 温度 pa 0.1 MPa 大气压 pinj 3 MPa[20] 注入压力 σ0 6 MPa[20] 围压 D 2×10-8 m²/s[20] CO2扩散系数 w 12 m-2[20] 形状因子 εLm 0.0518[20] Langmuir应变常数 VLm 0.017 m3/kg[20] Langmuir体积常数 PLm 4.48 MPa[20] Langmuir压力常数 φf0 1.2%[20] 初始裂隙孔隙率 φm0 2%[20] 初始基质孔隙率 kf0 1×10-17 m2[20] 初始裂隙渗透率 表 2 用于现场开采煤层气数值模拟的参数
Table 2 Parameters used for numerical simulation of CBM
-
[1] XIAO Z Y, WANG G, WANG C S, et al. Permeability evolution and gas flow in wet coal under non-equilibrium state: Considering both water swelling and process-based gas swelling[J]. International Journal of Mining Science and Technology, 2023, 33(5): 585-599. doi: 10.1016/j.ijmst.2022.11.012
[2] ROBERTSON E. Measurement and modeling of sorption- induced strain and permeability changes in coal[R]. Idaho Falls: Idaho National Laboratory, 2005.
[3] SEIDLE J P, JEANSONNE M W, ERICKSON D J. Application of matchstick geometry to stress dependent permeability in coals[C]// Proceedings of SPE Rocky Mountain Regional Meeting. Casper, 1992.
[4] PALMER I, MANSOORI J. How permeability depends on stress and pore pressure in coalbeds: a new model[C]// SPE Annual Technical Conference and Exhibition. Colorado, 1996: 539-544.
[5] SHI J Q, DURUCAN S. 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. doi: 10.1023/B:TIPM.0000018398.19928.5a
[6] WU Y, LIU J S, ELSWORTH D, et al. Dual poroelastic response of a coal seam to CO2 injection[J]. International Journal of Greenhouse Gas Control, 2010, 4(4): 668-678. doi: 10.1016/j.ijggc.2010.02.004
[7] WU Y, LIU J S, CHEN Z W, et al. A dual poroelastic model for CO2-enhanced coalbed methane recovery[J]. International Journal of Coal Geology, 2011, 86(2/3): 177-189.
[8] KUMAR H, ELSWORTH D, MATHEWS J P, et al. Effect of CO2 injection on heterogeneously permeable coalbed reservoirs[J]. Fuel, 2014, 135: 509-521. doi: 10.1016/j.fuel.2014.07.002
[9] 张宏学. 页岩储层渗流-应力耦合模型及应用[D]. 徐州: 中国矿业大学, 2015. ZHANG Hongxue. Seepage-Stress Coupling Model of Shale Reservoir and Its Application[D]. Xuzhou: China University of Mining and Technology, 2015. (in Chinese)
[10] LIU J S, CHEN Z W, ELSWORTH D, et al. Interactions of multiple processes during CBM extraction: a critical review[J]. International Journal of Coal Geology, 2011, 87(3/4): 175-189.
[11] LIU H H, RUTQVIST J. A new coal-permeability model: internal swelling stress and fracture–matrix interaction[J]. Transport in Porous Media, 2010, 82(1): 157-171. doi: 10.1007/s11242-009-9442-x
[12] LU S Q, CHENG Y P, LI W. Model development and analysis of the evolution of coal permeability under different boundary conditions[J]. Journal of Natural Gas Science and Engineering, 2016, 31: 129-138. doi: 10.1016/j.jngse.2016.02.049
[13] LIU T, LIN B Q, YANG W. Impact of matrix–fracture interactions on coal permeability: Model development and analysis[J]. Fuel, 2017, 207: 522-532. doi: 10.1016/j.fuel.2017.06.125
[14] 肖智勇, 王长盛, 王刚, 等. 基质-裂隙相互作用对渗透率演化的影响: 考虑基质变形和应力修正[J]. 岩土工程学报, 2021, 43(12): 2209-2219. doi: 10.11779/CJGE202112007 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): 2209-2219. (in Chinese) doi: 10.11779/CJGE202112007
[15] LIU J S, WANG J G, CHEN Z W, et al. Impact of transition from local swelling to macro swelling on the evolution of coal permeability[J]. International Journal of Coal Geology, 2011, 88(1): 31-40. doi: 10.1016/j.coal.2011.07.008
[16] 王刚, 肖智勇, 王长盛, 等. 基于非平衡状态的煤层中气体运移规律研究[J]. 岩土工程学报, 2022, 44(8): 1512-1520. doi: 10.11779/CJGE202208016 WANG Gang, XIAO Zhiyong, WANG Changsheng, et al. Gas transport in coal seams based on non-equilibrium state[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(8): 1512-1520. (in Chinese) doi: 10.11779/CJGE202208016
[17] PENG Y, LIU J S, WEI M Y, et al. Why coal permeability changes under free swellings: New insights[J]. International Journal of Coal Geology, 2014, 133: 35-46. doi: 10.1016/j.coal.2014.08.011
[18] WEI M Y, LIU J S, SHI R, et al. Long-term evolution of coal permeability under effective stresses gap between matrix and fracture during CO2 injection[J]. Transport in Porous Media, 2019, 130(3): 969-983. doi: 10.1007/s11242-019-01350-7
[19] WEI M Y, LIU J S, ELSWORTH D, 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. doi: 10.1016/j.ijrmms.2018.12.021
[20] WEI M Y, LIU J S, ELSWORTH D, et al. Impact of equilibration time lag between matrix and fractures on the evolution of coal permeability[J]. Fuel, 2021, 290: 120029. doi: 10.1016/j.fuel.2020.120029
[21] ZENG J, LIU J S, LI W, et al. Evolution of shale permeability under the influence of gas diffusion from the fracture wall into the matrix[J]. Energy & Fuels, 2020, 34(4): 4393-4406.
[22] ZENG J, LIU J S, LI W, et al. A process-based coal swelling model: Bridging the gaps between localized swelling and bulk swelling[J]. Fuel, 2021, 293: 120360. doi: 10.1016/j.fuel.2021.120360
[23] HUANG Y F, LIU J S, ELSWORTH D, et al. A transient dual porosity/permeability model for coal multiphysics[J]. Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 2022, 8(2): 40. doi: 10.1007/s40948-022-00348-8
[24] LIU X X, CHEN L, SHENG J C, et al. A Non-Equilibrium multiphysics model for coal seam gas extraction[J]. Fuel, 2023, 331: 125942. doi: 10.1016/j.fuel.2022.125942
[25] 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. doi: 10.1016/j.coal.2018.01.005
[26] 张宏学, 刘卫群. 非平衡解吸状态下页岩气储层渗透率演化机制[J]. 岩土力学, 2021, 42(10): 2696-2704. ZHANG Hongxue, LIU Weiqun. Permeability evolution mechanism of shale gas reservoir in non-equilibrium desorption state[J]. Rock and Soil Mechanics, 2021, 42(10): 2696-2704. (in Chinese)
[27] GAO Q, LIU J S, HUANG Y F, et al. A critical review of coal permeability models[J]. Fuel, 2022, 326: 125124. doi: 10.1016/j.fuel.2022.125124
[28] DETOURNAY E, CHENG A H D. Fundamentals of poroelasticity[M]//Analysis and Design Methods. Amsterdam: Elsevier, 1993: 113-171.
[29] JIANG C Z, ZHAO Z F, ZHANG X W, et al. Controlling effects of differential swelling index on evolution of coal permeability[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2020, 12(3): 461-472. doi: 10.1016/j.jrmge.2020.02.001
[30] 吴宇. 煤层中封存二氧化碳的双重孔隙力学效应研究[D]. 徐州: 中国矿业大学, 2010. WU Yu. Study on Double Pore Mechanical Effect of Carbon Dioxide Storage in Coal Seam[D]. Xuzhou: China University of Mining and Technology, 2010. (in Chinese)
[31] CHEN M, HOSKING L J, SANDFORD R J, et al. Dual porosity modelling of the coupled mechanical response of coal to gas flow and adsorption[J]. International Journal of Coal Geology, 2019, 205: 115-125. doi: 10.1016/j.coal.2019.01.009