煤层甲烷俗称煤矿瓦斯，它是引起煤矿瓦斯灾害的根源。根治煤矿瓦斯灾害和实现瓦斯资源高效利用的根本措施就是对其进行高效抽采。目前，高瓦斯单一煤层采取密集钻孔、水压致裂、水力冲孔、开采解放层强化抽采、水力割缝和注水驱气等方法进行煤层瓦斯的抽采^[1–2]，上述强化抽采瓦斯技术对提高中国煤层瓦斯采收率有重要作用，除存在一定安全隐患的深孔松动爆破外，这些增透方法的主要特点是在煤层中形成相对单一的通道来提高瓦斯抽采率，而煤层是由吸附甲烷约90%的微孔基质系统和割理裂隙系统组成的孔隙裂隙双重介质，煤岩储层的孔裂隙系统影响煤层甲烷在储层中富集和渗流的机理，因此，充分改造煤体割理裂隙系统的渗透性，促进基质微孔系统的瓦斯解吸扩散速度，是提高煤层瓦斯抽采率的关键。

超临界CO₂是一种优良的有机溶剂，具有较高的溶解能力和良好的流动性，能迅速渗透到多孔介质的微观孔裂隙中，除了具有非超临界CO₂的驱替置换作用，还具有明显优于非超临界CO₂的增透作用，溶解煤体中的部分有机物，充分改善煤微观孔隙和裂隙的形态、大小、充填情况、开裂程度等微观结构特征，加快煤层气在孔裂隙中的渗流扩散，提高低渗透煤层的渗透性^[3–4]。目前已见相关超临界CO₂增透方面的实验研究。孙可明等^[5]在低渗透煤层中注入超临界CO₂，得到了超临界CO₂作用后煤体渗透率和渗透系数随有效体积应力、温度以及孔隙压力的变化规律。杨涛等^[6]进行了超临界CO₂抽提低渗透煤试验，验证了超临界CO₂对提高煤层低渗透性的可行性。岳立新等^[7]在宏观超临界CO₂增透实验基础上进行了增透煤的微观成像实验，宏微观结合研究了超临界CO₂作用后煤微观孔隙率和渗透率的关系。上述研究成果为利用超临界CO₂改造低渗透煤层的渗透性及煤层气的高效抽采提供了重要的理论依据。目前，利用超临界CO₂改善低渗透煤层的渗透性大多还仅限于实验室的试件研究，在数值模拟方面多见于非超临界CO₂在煤体中渗流规律的数值模拟研究。叶建平等^[8]基于渗透率变化规律建立数值模型，分析了CO₂注入量、频率和注入方式对井组或单井的产量、采收率的影响。唐书恒等^[9]开展了注CO₂驱替煤层甲烷的实验模拟研究，分析了驱替过程中CH₄组分的产出特征。冯启言等^[10]建立了CO₂–CH₄气固耦合的有限元数值模型，模拟分析了CO₂–CH₄二元气固耦合机理。可见CO₂在非超临界状态下对煤体渗透性影响方面的数值模拟已经有了一定的研究基础，但现有模拟研究均基于非超临界CO₂的驱替置换作用建立相应的数值模型，而当CO₂达到超临界状态后，对煤体具有较强的溶解增透作用，除了具有非超临界CO₂的驱替置换作用外，还可以充分改造低渗透煤层割理裂隙系统的渗透性来提高裂隙系统瓦斯渗流速度和促进基质微孔系统的瓦斯解吸扩散速度，以有效提高煤层瓦斯的抽采率，这是超临界CO₂在增透低渗透煤层方面优于非超临界CO₂的显著特点。鉴于此，作者基于超临界CO₂增透后煤微观孔裂隙演化特征，考虑孔隙压力、温度、体积应力及超临界CO₂溶解增透作用的影响，建立超临界CO₂增透煤热流固耦合力学模型，实现低渗透煤层注超临界CO₂增透规律数值模拟，研究成果将对利用超临界CO₂改造低渗透煤层的渗透性具有非常重要的战略意义。

1 超临界CO₂增透煤热流固耦合力学模型

在力学耦合模型建立之前需进行一些假设：1）煤体变形符合小变形；2）煤体中仅流动单相的CO₂流体，并呈饱和状态；3）CO₂在煤体孔裂隙中流动时满足达西定律；4）煤体中热量的传递方式为流体与煤体骨架之间的传导和流体间的对流；5）CO₂流体和煤体骨架之间可瞬时达到局部的热平衡；6）惯性力和体积力忽略。

1.1 孔隙率和渗透率动态演化模型

依据前期超临界CO₂增透煤微观孔裂隙随温度及孔隙压力的演化规律^[11]，以及相关文献关于体积应力对孔隙率的影响^[12]，得到了煤微观孔隙率随孔隙压力P、温度T及体积应力 ${\sigma _{\rm{V}}}$ 变化的孔隙率方程。

$\phi {}_1 = {\phi _0}{{\rm e}}^{aP - bT - c{\sigma _{{\rm V}}}}$

(1)

当超临界CO₂作用于煤体后，除了温压条件及体积用力对煤微观孔裂隙的影响外，由于超临界CO₂较强的溶解增透作用，能溶解煤中的部分有机化合物，从而改变煤微观孔裂隙的密度、形态以及连通性等特征，因此煤微观孔隙率的变化还应该考虑超临界CO₂的溶解增透作用。超临界CO₂的溶解度满足Chrastil方程^[13–14]：

$\ln y = {\alpha _1}\ln \rho + {\alpha _2}/T + {\alpha _3}$

(2)

即：

$y = {{\rm e}^{{\alpha _1}\ln \rho + {\alpha _2}/T + {\alpha _3}}}$

(3)

式中： $y$ 为物质在超临界CO₂中的溶解度，kg/m³； ${\alpha _1}$ 、 ${\alpha _2}$ 、 ${\alpha _3}$ 为参数，其中 ${\alpha _1} = 4.5$ ， ${\alpha _2} = - 11041$ ， ${\alpha _3} = - 1.2$ ； $T$ 为温度，K； $\rho $ 为超临界CO₂的密度，kg/m³，按照式（4）计算^[15]：

$\rho = \frac{{\rho {}_{\rm{c}}P{T_0}}}{{T{P_0}{\textit{Z}}}}$

(4)

式中： ${\textit{Z}}$ 为CO₂的压缩系数； ${P_0}$ 为标况压力， ${P_0}$ =0.1 MPa； $P$ 为孔隙压力，MPa； $T$ 为温度，K； ${T_0} = 273\;{\rm K}$ ，令 $A = \displaystyle\frac{{{\rho _{\rm{a}}}T{}_0}}{{{P_0}{\textit{Z}}}}$ ， $A = \displaystyle\frac{{{\rho _{\rm{a}}}T{}_0}}{{{P_0}{\textit{Z}}}}$ ，则 $\rho = \displaystyle\frac{{AP}}{T}$ 。

则由于超临界CO₂溶解增透作用引起的煤微观孔隙率为：

$\phi {}_2 = {\phi _0}{{\rm e}^{{\alpha _1}\ln \rho + {\alpha _2}/T + {\alpha _3}}}$

(5)

由式（1）和（5），得出考虑孔隙压力、温度、体积应力及超临界CO₂溶解增透作用影响的孔隙率方程为：

$\phi = {\phi _0}{{\rm e}^{aP - bT - c{\sigma _{\text{V}}} + {\alpha _1}\ln AP/T + {\alpha _2}/T + {\alpha _3}}}$

(6)

式中： ${\phi _0}$ 为初始孔隙率； $a$ 、 $b$ 、 $c$ 为拟合系数，其中 $a$ 、 $b$ 由前期微观实验得到， $a = 0.1189$ ， $b = 0.1913$ ， $c$ 参照文献^[15]取值0.0001； ${\sigma _{\rm{V}}}$ 为体积应力， ${\sigma _{\rm{V}}} = {\sigma _1} + 2{\sigma _2}$ ，其中 ${\sigma _1}$ 和 ${\sigma _2}$ 分别为煤单元体中的第一主应力和第二主应力，MPa。

根据Palmer等提出的煤岩渗透系数与孔隙率变化规律呈三次指数关系^[16–17]，得到渗透率方程为：

$k = {k_0}{({{\rm e}^{\scriptstyle aP - bT - c{\sigma _{\rm{V}}} + {\alpha _1}\ln AP/T + {\alpha _2}/T + {\alpha _3}}})^3}$

(7)

式中， $k$ 和 ${k_0}$ 分别为渗透率和初始渗透率。

1.2 煤体变形场控制方程

1）煤体变形几何方程

${\rm d}{\varepsilon _{ij}} = \frac{1}{2}\left( {\frac{{\partial \left( {{\rm d}{u_i}} \right)}}{{\partial {x_j}}} + \frac{{\partial \left( {{\rm d}{u_j}} \right)}}{{\partial x{}_i}}} \right)$

(8)

式中： ${u_i}$ 为位移分量， ${\rm d}{\varepsilon _{ij}}$ 为总的应变增量。

2）平衡方程

$\frac{{\partial \left( {{\rm d}{\sigma _{ij}}} \right)}}{{\partial {x_j}}} + {\rm d}{F_i} = 0$

(9)

式中， ${\rm d}{F_i}$ 为体力分量增量。

3）弹塑性损伤本构方程

煤岩内部存在大量的微裂纹、微孔洞等缺陷，超临界CO₂作用于煤体后，煤体中微观孔裂隙不断演化，引起煤体的损伤，根据增透后煤微观孔隙率的变化定义损伤变量：

$d = d(\phi )$

(10)

由渗透率与孔隙率之间的3次方关系及渗透率和损伤之间的关系方程，可以推出孔隙率与损伤之间的函数关系：

$k = {k_0}{\left( {\frac{\phi }{{{\phi _0}}}} \right)^3} = {k_0}{{\rm e}^{{{\left[ {\left( {{A_1}d} \right)} \right]}^{{B_1}}}}}$

(11)

则：

$\begin{array}{l}d = \displaystyle\frac{1}{{{A_1}}}{\left[ {3\ln (\displaystyle\frac{\phi }{{{\phi _0}}})} \right]^{{B_1}}}= \displaystyle\frac{1}{{{A_1}}}{[3(aP - bT - c{\sigma _e} + {\alpha _1}\ln AP/T + {\alpha _2}/T + {\alpha _3})]^{{B_1}}}\end{array}$

(12)

式中： $d$ 为损伤变量； ${A_1}$ 和 ${B_1}$ 是常数，依据文献^[18] ${A_1}$ 取10， ${B_1}$ 取1.7。

则增量形式的弹塑性本构关系为：

$\left\{ {\begin{aligned}& {{\rm d}{\sigma _{ij}\rm } = \left( {1 - {d}} \right){ D}_{ijkl}^e\left( {{\rm d}{\varepsilon _{ij}} - {\rm d}\varepsilon _{ij}^p} \right) - {\delta _{ij}}\alpha {\rm d}P - {\beta _s}{\delta _{ij}}{\rm d}T},\\& {{\rm d}\varepsilon _{ij}^p = {\rm d}\lambda \displaystyle\frac{{\partial G\left( {{\sigma _{ij}}} \right)}}{{\partial {\sigma _{ij}}}}}\end{aligned}} \right.$

(13)

式中： ${ D}\mathop {_{ijkl}^e}\limits^{} $ 为初始弹性刚度矩阵； ${\rm d}\varepsilon _{ij}^e$ 为弹性应变增量； ${\rm d}\varepsilon _{ij}^p$ 为塑性应变增量； ${\beta _s}$ 为热膨胀系数； $\alpha $ 为Biot系数； ${\delta _{ij}}$ 为Kronnecker数； ${\rm d}\lambda $ 为硬化参数； $G\left( {{\sigma _{ij}}} \right)$ 为塑性势函数。塑性势函数 $G$ 选取Drucker-Prager双曲线函数。

1.3 渗流场方程

依据达西定律，忽略重力的影响，超临界CO₂流动速度矢量 $ { q} $ 的表达式为：

${ q}= - \frac{k}{\mu }\nabla P$

(14)

式中： $k$ 为渗透率， $\mu $ 为超临界CO₂的黏度。

煤体中包含的吸附态CO₂和游离CO₂的总质量为：

$m = \rho \phi + \left( {1 - \phi } \right){\rho _{\rm c}}{\rho _{\rm s}}\frac{{{V_{\rm L}}P}}{{P + {P_{\rm L}}}}$

(15)

式中， $\phi $ 为煤微观孔隙率； ${\rho _{\rm c}}$ 为标况下CO₂的密度，kg·m^–3； ${\rho _{\rm s}}$ 为煤体密度，kg·m^–3； ${V_{\rm L}}$ 为Langmuir体积常数； ${P_{\rm L}}$ 为Langmuir压力常数。

由式（15）得到超临界CO₂含量的变化率为：

$\begin{aligned}[b]\displaystyle\frac{{\partial m}}{{\partial t}}= & \left( {\displaystyle\frac{{A\phi }}{T} + \left( {1 - \phi } \right)\displaystyle\frac{{{\rho _{\rm a}}{\rho _{\rm c}}{V_{\rm L}}{P_{\rm L}}}}{{{{\left( {P + {P_{\rm L}}} \right)}^2}}}} \right)\displaystyle\frac{{\partial P}}{{\partial t}}+\\ & \left( {\displaystyle\frac{{AP}}{T} - \displaystyle\frac{{{\rho _{\rm a}}{\rho _{\rm c}}{V_{\rm L}}P}}{{P + {P_{\rm L}}}}} \right)\displaystyle\frac{{\partial \phi }}{{\partial t}} - \displaystyle\frac{{A\phi P}}{{{T^2}}}\displaystyle\frac{{\partial T}}{{\partial t}}\end{aligned}$

(16)

式中：

$\frac{{\partial \phi }}{{\partial t}} = \phi \left( {(a + \frac{{{\alpha _1}}}{P} + c\bar \alpha )\frac{{\partial P}}{{\partial t}} - (b + \frac{{{\alpha _1}}}{{{T^2}}} + \frac{{{\alpha _2}}}{{{T^2}}})\frac{{\partial T}}{{\partial t}}} \right)$

(17)

由式（14）和（16），依据质量守恒定律得渗流场方程为：

$\begin{array}{l}\left( {\displaystyle\frac{{A\phi }}{T} + \left( {1 - \phi } \right)\displaystyle\frac{{{\rho _{\rm a}}{\rho _{\rm c}}{V_{\rm L}}{P_{\rm L}}}}{{{{\left( {P + {P_{\rm L}}} \right)}^2}}} + \phi (\displaystyle\frac{{AP}}{T} - \displaystyle\frac{{{\rho _{\rm a}}{\rho _{\rm c}}{V_{\rm L}}P}}{{P + {P_{\rm L}}}})(a + \displaystyle\frac{{{\alpha _1}}}{P} + c\bar \alpha )} \right)\displaystyle\frac{{\partial P}}{{\partial t}} - \quad\quad\left( {\phi (\displaystyle\frac{{AP}}{T} - \displaystyle\frac{{{\rho _{\rm a}}{\rho _{\rm c}}{V_{\rm L}}P}}{{P + {P_{\rm L}}}})(b + \displaystyle\frac{{{\alpha _1}}}{{{T^2}}} + \displaystyle\frac{{{\alpha _2}}}{{{T^2}}}) + \displaystyle\frac{{A\phi P}}{{{T^2}}}} \right)\displaystyle\frac{{\partial T}}{{\partial t}} + \nabla \cdot (\displaystyle\frac{{APk}}{{T\mu }} \cdot \nabla P) = 0\end{array}$

(18)

式中： $\phi = {\phi _0}{e^{aP - bT - c{\sigma _{\rm V}} + {\alpha _1}\ln \scriptstyle\frac{{AP}}{T} + {\alpha _2}/T + {\alpha _3}}}$ ， $k = {k_0}{\left( {{{\rm{e}}^{aP - bT - c{\sigma _{\rm{v}}} + {\alpha _1}\ln AP/T + {\alpha _2}/T + {\alpha _3}}}} \right)^3}$ ，

$\begin{array}{l}\displaystyle\frac{{\partial \phi }}{{\partial t}} ={\phi _0}{{\rm{e}}^{aP - bT - c{\sigma _{\rm{v}}} + {\alpha _1}\ln AP/T + {\alpha _2}/T + {\alpha _3}}}\left( {(a + \displaystyle\frac{{{\alpha _1}}}{P} + c\alpha )\displaystyle\frac{{\partial P}}{{\partial t}} - (b + \displaystyle\frac{{{\alpha _1}}}{{{T^2}}} + \displaystyle\frac{{{\alpha _2}}}{{{T^2}}})\displaystyle\frac{{\partial T}}{{\partial t}}} \right)。\end{array}$

1.4 温度场控制方程

假定热传导方程满足傅立叶定律，因为CO₂流体和固体之间瞬时达到局部热平衡，则CO₂在煤体中流动的总的能量守恒方程为：

$\begin{array}{l}\displaystyle\frac{\partial }{{\partial t}}\left[ {(1 - \phi ){\rho _{\rm s}}{c_{\rm s}}\nabla T + \phi \;{\rho _{\rm c}}{c_{\rm c}}\nabla T} \right] + (1 - \phi )T{\beta _{\rm s}}K\displaystyle\frac{{\partial {\varepsilon _{\rm V}}}}{{\partial t}} + \nabla \cdot \left( {{\rho _{\rm c}}c{}_{\rm c}\; { q}\Delta T\;} \right) = \;\;\nabla \cdot \left( {\lambda \nabla T} \right)\end{array}$

(19)

式中： ${\rho _{\rm s}}$ 为煤体密度； ${c_{\rm S}}$ 为煤体比热； ${\beta _{\rm S}}$ 为热膨胀系数； ${\varepsilon _{\rm V}}$ 为体积应变； ${\lambda _{\rm s}}$ 为煤体导热系数； $K$ 为体积模量， $K = E/3(1 - 2\mu )$ ， $E$ 为弹性模量， $\mu $ 为煤体的泊松比； ${\rho _{\rm c}}$ 为CO₂的密度； ${c_{\rm c}}$ 为CO₂比热； ${\lambda _{\rm c}}$ 为CO₂导热系数； $\lambda = \left( {1 - \phi } \right){\lambda _{\rm s}} + \phi {\lambda _{\rm c}}$ 。

1.5 定解条件

1）温度场边界条件及初始条件

初始条件： $T\left( {{t_0}} \right) = T{}_0$ ；

边界条件： $T = \tilde T$ ，（在 ${\varGamma _1}$ 内）； $k\nabla Tn = q$ ，（在 ${\varGamma _2}$ 内）。

2）渗流场边界条件及初始条件

初始条件： $P\left( {{t_0}} \right) = {P_0}$ ；

边界条件： $\sigma \cdot n = s$ （在 ${\varGamma _\sigma }$ 内）。

3）变形场边界条件及初始条件

初始条件： $u\left( {{t_0}} \right) = {u_0}$ （在 $\varOmega $ 内）；

边界条件： $\sigma \cdot n = s$ （在 ${\varGamma _\sigma }$ 内）； $u = \tilde u$ （在 ${\varGamma _u}$ 内）。

2 超临界CO₂增透煤数值模拟 2.1 马脊梁矿概况

山西大同马脊梁矿瓦斯含量丰富，随着采煤深度的逐年增加，该矿井的瓦斯涌出量不断增多，但由于煤层渗透性差致使瓦斯抽采率较低。根据马脊梁矿13 250工作面的煤层地质情况及瓦斯抽采现状，概化一个算例进行数值模拟。该工作面走向长度为1 860 m，工作面长185 m，平均倾角2.8°，煤层厚度2 m，处于约400 m的埋深。煤层相关参数为：瓦斯压力0.65 MPa，瓦斯含量为11.802 5 m³/t，煤孔隙率为0.92%，透气性系数为1.59 m²/(MPa·d)。在胶带巷和轨道巷内布置本煤层钻孔，采用垂直巷道走向钻孔，孔深100 m，孔距2 m，开孔孔径150 mm，终孔孔径113 mm，扩孔深度15 m，在工作面切眼口以外6 m处施工钻孔，钻孔开孔高度为底板以上4 m的位置。

2.2 数值模型建立

1）数值模拟几何模型

将煤层注CO₂转化为平面应变问题，模型的几何长度4 m，宽度2 m，钻孔半径为0.075 m。模型网格划分采用2阶Quad型Structured网格，单元类型为标准线性CPE4P，将模型划分成3 456个几何单元，几何模型及网格划分如图1所示。

2）边界条件

依据土的容重及埋深，计算初始地应力为10 MPa。煤层初始孔隙压力为0.65 MPa，初始温度为293 K，注入CO₂的孔隙压力为10 MPa，温度为338 K，模型的左右边界、底边和钻孔边缘施加相应的对称约束。模拟中的主要参数见表1。

3）模拟方案

在挖孔过程中，孔周必然形成应力集中，从而使得地应力平衡受到影响，为了取消挖孔对地应力平衡的影响，在模拟时进行三步分析：第一步，地应力平衡步，建立未开挖的模型，设定边界条件，施加载荷，进行初始地应力平衡计算；第二步，钻孔开挖卸载，将模型中钻孔单元移除，此时由于开挖卸载将会引起失衡；第三步，对开挖后的模型注入超临界CO₂，进行计算。

2.3 模拟结果分析

1）应力场变化规律

体积应力分布如图2所示，一定温度和压力的超临界CO₂注入后，由于钻孔引起孔周煤体的应力集中，以及超临界CO₂与煤体的流固耦合作用，注气孔周围体积应力分布呈现出明显的方向性，且各向传播的强度不等。随着注气时间的增加，煤体内超临界CO₂的流量增多，煤体内体积应力呈降低的变化规律。

2）温度场变化规律

温度随时间的变化规律如图3所示，煤层初始温度为298 K，当温度为338 K的超临界CO₂注入煤体后，逐渐向周围煤体孔裂隙内渗透流动，并将高温流体具有的热量经热传导方式向周围煤体传递，使煤体温度逐渐升高，并不断向煤体纵深方向传递。当注气时间t=20 h时，超临界CO₂进入孔周煤体孔裂隙中的流量较少，对孔周围煤体温度影响的范围较小，随着注气时间的延长，当注气时间延至90 h时，超临界CO₂流体通过渗流扩散进入更大范围煤体的孔裂隙中，温度传播的范围进一步扩大，孔周径向约1.2 m范围内的煤体温度都达到了CO₂的超临界温度。

图4为注气孔周围煤体温度沿图示路径传播的规律曲线，距注气孔较近范围煤体内，温度曲线变化较陡，温度增加较快，这是因为近注气孔煤体内CO₂流量较多，与煤体发生热传导的速率较快，煤体温度变化较为明显。随着距注气孔距离的增加，煤体孔裂隙中超临界CO₂流体的流量逐渐减少，通过热传导传递给煤体的热量逐渐减少，煤体温度增加的速率和幅度逐渐降低。对比可见，注气时间越长，煤体温度升高的幅度越大。

3）渗流场变化规律

图5为单孔注气过程中，孔隙压力随时间传播的变化云图，超临界CO₂注入后，注气孔周围煤体内孔隙压力有了明显增加，在较短的注气时间t=20 h时，注气孔周围明显形成了由10 MPa孔隙压力向周围煤体逐渐减小的压力梯度，并呈现出近似椭圆形的孔隙压力传递区。随着注气时间的持续，超临界CO₂逐渐向孔周煤体孔裂隙中渗流扩散，煤体内孔隙压力进一步提高，孔隙压力传递的范围不断扩大。当注气时间达到t=72 h时，孔隙压力传播范围达到了注气孔周围径向约1 m的区域。可见，为取得较好的增透效果，需要保证一定的注气时间。

图6为孔隙压力变化的路径图，在一定范围内，煤体孔隙压力较初始孔隙压力0.67 MPa增加明显，不同注气时间下，孔隙压力的传播呈现出明显的非线性变化规律，均表现为随着距注气孔距离的增加，孔隙压力增加的速率逐渐减小，直至趋于稳定。对比不同注气时间下孔隙压力的变化可以看到，注气时间越长，孔隙压力增加的幅度越大，说明在其它注气条件相同的条件下，增加注气时间，有助于超临界CO₂在煤体中的渗流扩散，扩大孔隙压力传递的范围，提高增透效果。

4）煤微观孔隙率变化规律

图7为煤微观孔隙率随时间的变化云图，超临界CO₂的注入引起了注气孔周围煤微观孔裂隙的演化，孔隙率逐渐增加。随着注气时间的延长，孔隙率增加的范围逐渐增扩大，当t=80 h时，在注气孔周围1.67 m范围内的孔隙率都有了一定程度的增加。这是由于煤是一种多孔隙、裂隙的介质，一定压力梯度的超临界CO₂注入后，对煤体孔隙、裂隙产生向外扩张的作用，使孔裂隙尺寸增大，渗流通道数量增加，此外由于超临界CO₂对煤体的增透作用，使煤体内新的孔裂隙不断萌生、演化，并与原有的孔裂隙相互贯通、扩展，也促进了煤微观孔裂隙的有效发育。

从孔隙率随注气时间变化的等值线图（图8）可以看到，注气孔周围孔隙率等值线变化的疏密程度不同，距注气孔越近，等值线越密集，随着距注气孔距离的增加，孔隙率等值线逐渐稀疏，说明距注气孔越近，孔隙率增加的速率越快；距注气孔越远，孔隙率增加的速率逐渐减慢。这是因为距离注气孔越近，煤体孔隙内超临界CO₂的流量越多，对煤体的增透致裂作用越强，煤体内演生出的孔裂隙越多。

超临界CO₂注入后煤微观孔隙率的演化说明了超临界CO₂在改造低渗透煤层渗透性方面的特有优势，除了具有非超临界CO₂对煤层瓦斯的驱替置换解吸作用外，还具有明显的溶解作用，显著改善了煤微观孔裂隙的连通性，提高了煤的微观孔隙率，提供了更多的供煤层气运移的渗流通道，为有效提高煤层气的采收率开辟了新的方向。

5）煤体损伤变化规律

图9为煤体损伤随注气时间的变化云图，注气孔边缘煤体损伤最为明显，当t=15 h时，注气孔边缘的损伤量达到了0.36，当注气时间为96 h时，孔边缘处损伤增大到了0.38，这是由于注气孔边缘煤体内孔隙压力较高，超临界CO₂流体在较高压力下极易进入煤体的微孔隙和微裂隙中，促使煤微观孔裂隙演化，引起煤体损伤，并逐步向周围煤体中扩展，损伤增加。随着注气时间的延长，煤体损伤范围不断扩大，并向注气孔周围煤体内呈不规则波形传播，具有明显的方向性。

图10为沿图示路径方向煤体损伤的变化曲线，在不同的注气时间条件下，由注气孔向周围煤体，损伤的变化规律是一致的：即距注气孔越近，损伤程度越大，损伤增加的越快。注气时间越长，超临界CO₂对煤体孔裂隙的致裂增透作用越强，损伤增加的幅度越大。

损伤的模拟规律说明了经超临界CO₂作用后，煤体内部原有的各级孔裂隙不断演化，局部产生开裂，萌生新的孔隙，微观孔裂隙的演化发展，最终引起煤体损伤。

6）渗透系数变化规律

图11为渗透系数随注气时间的变化云图，超临界CO₂注入后煤体渗透系数相对于初始态提高了3个数量级，进一步证实了超临界CO₂在改造低渗透煤层渗透性方面的显著优点。模拟结果与前期实验结果是一致的^[14]。

随着注超临界CO₂时间的延长，煤体渗透系数由注气孔向周围煤体纵深方向进一步扩展，范围逐渐增大。当注气时间达到t=85 h时，注气孔周围较大区域内煤体的渗透系数较注气前有了较为明显的提高，可见注气时间越长，注气孔周围煤微观孔隙率增加的范围越大，渗透系数向周围煤体传播的范围随之扩大。

由注气孔向周围煤体渗透系数传播的强度看，距注气孔越近，渗透系数增加的程度最大，随着距离的增加，渗透系数传递的强度减弱，这是由于距离注气孔越远，超临界CO₂浓度逐渐减小，对煤体的增透致裂作用逐渐减弱，孔裂隙的演化速率逐渐减小，渗透系数增加的程度也呈现减小的趋势。

图12为渗透系数变化的路径图，在距注气孔较近范围内，渗透系数增加的速率较快，随着距离的增加，渗透系数增加的速率逐渐趋于平稳。当注气时间t=80 h时，渗透系数的增幅明显变大。在近注气孔区域，由于应力集中的影响，煤体渗透系数波动较大。不同注气时间下渗透系数的对比说明，增加注气时间，有助于煤微观孔裂隙的演化，使渗透系数的增幅加大，有效提高煤体的渗透性。

3 结　论

依据增透煤微观孔裂隙演化规律，建立超临界CO₂增透煤热流固耦合力学模型，针对山西大同马脊梁矿13250工作面瓦斯抽采现状，建立平面几何应变模型，实现低渗透煤层注超临界CO₂增透规律数值模拟，得到以下一些结论：

1）考虑孔隙压力、温度、体积应力及超临界CO₂溶解增透作用影响建立的增透煤热流固耦合力学模型和依据增透煤微观孔隙率的变化确定损伤变量的方法，对利用超临界CO₂改造低渗透煤层渗透性的方法是有效的。

2）煤层注超临界CO₂过程中，注气孔周围煤体的体积应力、温度及孔隙压力发生了明显的变化，且随距注气孔距离的增加，影响程度逐渐减弱，并趋于稳定。

3）超临界CO₂的注入促进了煤微观孔裂隙的演化，注气孔周围煤微观孔隙率较注气前提高两个数量级，损伤程度逐渐增大，煤体渗透系数较注气前提高3个数量级。

4）随着距注气孔距离的增加，煤微观孔隙率和损伤增加的速率逐渐减慢，渗透系数增加的程度呈降低的变化规律；随着注气时间的延长，煤微观孔裂隙演化程度提高，孔隙率和损伤程度增加，为流体提供更多的渗流通道，有效提高了低渗透煤层的渗透性。