煤层气是一种非常规天然气， 与煤相伴而生赋存于煤层之中。作为一种新型能源， 煤层气资源具有很大的发展潜力。由于煤层储气层渗透率往往较低， 且煤层具有低孔、 低压的特点， 使得煤层气开采难度增大、 开采效率降低， 必须采取工程措施以提高煤层气的开采量。目前， 最常用的方法是水力压裂技术， 但由于不同地域的煤层具有不同的地质特征， 水力压裂技术不能适用于所有的煤层。如果将人工冻结技术引入煤层改造之中， 通过冰水相变拓展、 贯通裂纹， 形成裂纹网络， 进而可达到增加煤层气产出的效果。

煤岩作为一种复杂的地质材料， 原始存在着大量不规则跨尺度的裂隙， 在外载作用下， 煤岩内部原有的裂隙会发生变化， 局部还可能形成新的裂隙。随着裂隙结构的发展， 细小裂隙贯通形成宏观裂缝， 进而导致煤岩体发生失稳破坏， 因此， 煤岩的变形破坏过程实质上是一个十分复杂的损伤演化过程^［1］。此外， 随着人工冻结技术的发展以及冻结技术在煤炭、 煤层气开采中的应用， 冻结状态下煤岩的损伤破坏行为成为岩体力学研究的难点问题^［2-5］。

冻结煤岩与冻土有着紧密的联系， 但又存在着明显的差异， 其根本区别在于煤岩体中存在裂隙， 裂隙将整个岩体分割成不连续的块体， 因而不能用连续介质力学方法（如有限元法）来模拟岩体的力学变形行为。同时， 由于裂隙中冰的胶结或润滑作用， 不断地改变了岩体的力学性质及变形特性^［6-8］， 因此， 冻结煤岩的损伤破坏过程更为复杂。早期的研究大多通过宏观试验来分析岩体的力学变形行为， 如进行岩体的单轴压缩、 单轴拉伸、 常规三轴、 剪切、 劈裂试验等条件力学性质方面的研究^［9-12］。此外， 代高飞等^［13］利用CT成像技术研究了脆性煤岩单轴压缩破坏全过程的细观损伤规律； 杨永杰等^［14］基于MTS试验机和扫描电镜， 从强度和变形特性的角度研究了煤岩的微细观演化过程。这些研究在一定程度上揭示了煤岩在荷载作用下损伤演化特性。

随着计算机技术的发展， 人们更多地依赖计算机模拟手段来分析加卸载条件下裂隙岩体的变形、 强度与破坏特性。唐春安等^［15-16］研发了RFPA软件， 将岩体材料性质引入到数值模拟中， 实现了对岩体破坏过程的全程模拟。钟卫平等^［17-18］运用ANSYS软件对裂隙岩体在压剪荷载作用下的应力场及位移场进行了数值模拟研究， 并利用断裂力学理论分析了纯剪切应力状态和双向压应力状态作用下岩体的断裂机制。还有一些能量判别准则被提出并引入到数值模拟中， 从而实现了煤岩体损伤及裂纹扩展的模拟^［19-21］。以上对裂隙岩体的计算模拟都采用了基于连续介质理论的有限元方法， 对裂隙岩体的损伤破坏过程及力学性质有了一定的认识。但是有限元方法受限于连续介质的假设， 对于裂隙岩体这种非完全连续又非完全散体的结构来说， 存在一定的限制条件。

本文从非连续介质理论角度出发， 采用离散元方法， 对冻结裂隙煤岩在加载条件下的变形、 强度及破坏特性进行数值分析， 并标定煤岩的细观力学参数， 以期通过本文的研究可以从非连续介质理论对冻结裂隙煤岩的力学性质及破坏机理进行分析， 为进一步拓宽离散元的应用领域， 并为今后裂隙煤岩或煤层气开采方面的研究提供参考。

煤岩是一种介于非连续介质又非完全散体介质之间的特殊材料， 由于其特殊的形成过程及后期的地质运动， 使得煤岩体中包含着大量的裂隙、 节理、 破碎带甚至是断层等不连续的结构， 这使得煤岩既不完全连续， 也不完全散体。若仅用传统的连续介质力学去解释岩体裂隙、 节理、 破碎断裂等非连续力学行为， 则会导致分析结果与实际状态之间存在较大差异。本文使用离散元模拟的方法来描述这种不完全属于连续介质的材料， 把裂隙煤岩看成是由离散的煤块和煤块间的节理面组成的材料， 煤块既能移动、 转动和变形， 节理面也可以被压缩、 分离或滑动， 从而可以较为真实的模拟煤岩中的不连续面。

由于离散单元法将煤岩看成是由单独的煤块通过节理面连接而成的材料， 所以需要描述煤块与煤块间的力学接触行为。在离散元理论中， 相互接触的分散块体间的相对运动会在连接处产生力和力矩， 并作用在相互连接的块体上， 因此， 与该模型相对应的接触力可以分解为在节理面上的法向力F_n 和切向力F_s， 如图1（a）。

在位移ΔU内弹性力增量为：

式中： k_n 表示接触处的法向刚度； k_s 表示接触处的切向刚度； ΔU_n 表示接触处的法向位移； ΔU_s 表示接触处的切向位移。

单就一个块体来说， 其受力状况如图1（b）， 相应的运动方程如下：

块体上的合力为：

合力矩为：

加速度为：

角加速度：

式中： F表示合力； F_i 表示分力； M表示总力矩； e_ij 表示基矢量； x_j 表示矢径； u¨表示加速度； θ¨表示角加速度； m表示块体的质量； I表示块体的转动惯量。

离散元的求解方法可分为静态松弛法和动态松弛法。本文采用动态松弛法， 动态松弛法是为解决隐式问题引入的方法， 采用向后时步迭代的计算思想， 把非线性静力学问题转化为动力学问题。该方法实质就是对临界阻尼振动方程进行逐步积分。动态松弛法的基本运动方程如下：

式中： u表示块体形心的位移； μ表示黏性阻尼系数； k表示刚度系数； f表示块体所受的外荷载； t表示时间。

块体离散单元法的求解过程不涉及到大型矩阵运算， 比较简单、 稳定， 因此可以节省计算时间， 同时允许单元发生大位移， 克服了传统有限元方法和边界单元法的小变形假设， 用来求解非线性问题十分方便。

岩体破坏遵循Mohr-Coulomb基本准则， 如图2所示， 图中曲线A到B点的屈服包络线f（σ₁，σ₃）=0对应的Mohr-Coulomb剪切屈服准则fs=0可以表示为：

曲线B点到C点对应的拉破坏准则ft=0可以表示为：

式中： φ表示内摩擦角； c表示黏聚力； σ^t 表示岩体的抗拉强度； 其中Nϕ=1-sinϕ1+sinϕ。

式（8）中， 当fs<0时， 发生剪切屈服； 式（9）中， 当fs>0时， 发生拉张屈服.

如果岩体的抗拉强度比σ₃小， 根据图2中的fs=0及σ₁=σ₃两条直线确定的抗拉强度最大值为：

由于煤岩属于节理结构的岩体， 因此选用能描述岩体内部节理的节理面接触-库伦滑移理论， 来模拟岩体中的裂隙、 断层及层面。对于完整结构面所受的法向应力满足下式：

式中： T_max表示破坏面上的拉力； σ表示节理面的抗拉强度； A_c 表示节理面面积。

结构面所受的最大剪力应满足：

式中： Fmaxs表示破坏面上的剪力； C表示破坏面上的黏聚力； φ表示破坏面上的内摩擦角。 当结构面张开或滑动后， 公式（11）和（12）中的相关指标用节理残余参数表示。

以山西沁水盆地煤系地层岩样的室内单轴试验为基础进行离散元模拟。将取得的整体岩块经过钻取岩心、 切割、 打磨几道程序， 制成直径为50 mm， 高为100 mm的圆柱形试样， 高径比符合国际岩体力学学会规定2.0~2.5的要求^［22］。将试样饱水后分别置于-15 ℃、 -10 ℃、 -4 ℃、 20 ℃的恒温环境下进行不同温度下的单轴压缩试验［图3（a）和3（b）］。轴向采用应力控制模式， 以0.5 MPa·s^-1的速度加载直至试件破坏。

根据地质勘查报告， 可知岩样的初始物理参数如表1。

三维离散元数值模型的建立分为三步： ①在试样制作完成后， 进行CT扫描（图4为20 ℃煤样的CT图像）， 获得岩样中上、 中间、 中下三个断面的裂隙度； ②建立离散体的结构空间， 煤岩试样的模拟尺寸采用试验的实际尺寸［图3（d）］； ③参考每个煤样CT扫描获得的裂隙度， 采用随机的方法来生成煤柱内的节理， 以保证节理的裂隙度与真实情况较为接近。

离散元数值试验是通过模型参数得到的结果来逼近室内试验的结果， 模型细观参数的选取不能直接采用宏观试验参数， 需要反复调试细观参数直至模拟曲线接近试验数据。表2中列出了不同温度下煤岩单轴压缩试验的模拟细观参数， 其中K为体积模量， G为剪切模量。

在-15 ℃、 -10 ℃、 -4 ℃、 20 ℃的恒温环境下， 进行煤岩的三维离散元数值模拟。图5（a）~5（d）分别为在四个温度下煤岩单轴压缩数值试验与室内试验应力-应变关系曲线的比较结果。从图中可以看出， 在单轴压缩情况下， 数值试验曲线与室内试验的应力-变关系存在一定的差异， 这主要是由于室内试验的煤柱存在天然缺陷， 以及煤柱两端面的不平整， 在加载的条件下会发生接触压密的过程， 随着煤柱缺陷的完全压密， 其应力-变曲线迅速增加。而理论模拟曲线不存在压密阶段， 因此模型煤柱加载后的应力应变关系呈线性增加， 出现短暂的塑性变形后产生破坏。由图5（a）~5（d）以看出， 单轴抗压强度的试验值与模拟理论值的相差不大， 总体来看， 单轴模拟曲线能够反映煤岩的压缩破坏过程。

从图5可以看出， 当煤柱的温度为-15 ℃时， 试验曲线在应变0.3%之前压缩变形较为明显， 之后应变曲线迅速增加， 并在应变为0.80%时发生脆性破坏， 而模拟曲线的屈服应变则为0.84%， 二者差别很小， 单轴抗压强度的试验值和理论值分别为5.23 MPa和4.96 MPa， 二者相差仅为5.2%； 当煤柱温度为-10 ℃时， 试验曲线在应变0.2%之前压缩变形较为明显， 随着煤柱压密， 在应变0.86%处发生脆性破坏， 而模拟曲线的屈服应变也为0.86%， 二者相同， 单轴抗压强度的试验值和理论值分别为5.18 MPa和4.66 MPa， 二者相差10.0%； 当煤柱温度为-4 ℃时， 试验曲线在应变0.4%之前的压缩变形比较明显， 之后应变曲线迅速增加， 并在应变为0.88%时发生脆性破坏， 而模拟曲线的屈服应变同样是0.88%， 二者相同， 单轴抗压强度的试验值和理论值分别为4.90 MPa和4.35 MPa， 二者相差11.2%； 当煤柱的温度为20 ℃时， 试验曲线在应变为4%之前的压缩变形较为明显， 随着煤柱压缩， 在应变为0.97%时发生脆性破坏， 而模拟曲线的屈服应变为0.94%， 二者差别很小， 单轴抗压强度的试验值和理论值分别为4.75 MPa和4.61 MPa， 相差12.4%。以上模拟结果和试验结果对比表明， 随着温度的降低， 模拟值更加接近真实值。

通过以上的分析发现， 随着煤柱温度的降低， 曲线的峰值应力增大。出现这样的现象是因为冻结煤岩中的未冻水会随着温度的降低而逐渐转变成冰， 而未冻水存在于煤岩的微裂隙中， 未冻水的状态变化使得煤岩的微裂隙得到冰的填充。由于煤岩的破坏主要发生在这些微裂隙、 节理面等上， 冰的黏结作用使得这些薄弱连接面得到了加固， 使其不易发生错动， 从而使得煤岩的强度得到提高。同时， 从四条曲线的走势来看， 煤岩的破坏类型均为脆性破坏， 这说明在此温度范围内， 煤岩的破坏类型与温度无关。但是温度的不同使得发生破坏的应变不同。当温度为-15 ℃、 -10 ℃、 -4 ℃、 20 ℃时， 煤岩离散元模拟曲线的应力-应变曲线峰值分别发生在应变为0.84%、 0.86%、 0.88%、 0.94%处。这说明随着温度的升高， 出现了破坏滞后现象， 这与煤岩的脆性变化相关。在低温环境下， 煤岩的脆性较为明显， 随着温度的升高， 其脆性有所降低， 但并不会改变其破坏形式。这是因为随着温度的降低， 煤岩中冰的含量增加， 冰对煤岩体的黏结作用增强， 从而使得煤岩体的整体性增强， 表现出更多的脆性特性。温度升高， 冰的含量减小， 冰对煤岩的黏结作用减弱， 表现出脆性减弱的性质。但是由于煤岩在变形发展过程中， 煤岩体内部缺陷不断发展， 微小裂隙逐渐贯穿形成大的裂缝， 使其变形迅速发展， 最终发生脆性破坏。

图6给出了温度为-15 ℃时， 模拟曲线在A点、 B点和C点的变形矢量图。A、 B、 C三点的选取分别为应变0.50%（弹性阶段应变）、 应变0.84%（破坏应变）、 应变1.20%（加载结束的应变）。取以上三个阶段的位移矢量变化可以描述煤岩体变形破坏过程中的变化规律。在弹性阶段， 煤岩体内的位移分布比较均匀， 从煤柱顶面到底面逐渐减小， 主要为竖向变形， 表面没有明显的凸出或凹陷， 煤样变形量较小， 呈现出完整柱状； 随着外荷载的增加， 煤样内的应力逐步达到其屈服状态， 达到屈服点的煤柱仍然呈现完整柱状， 没有明显破坏； 屈服后的煤岩从节理面开始发生滑移， 从而变形急剧增加， 煤样出现一条变形错断带， 并发生由表及里的剥离破坏， 此时， 水平位移大于竖向位移。发生这种现象的原因是， 由于荷载较小时， 煤柱内的节理面并未发生滑移， 各煤岩块仅发生弹性变形； 当荷载增大时， 破坏从裂隙处开始发展， 裂隙处相接触的块体发生相互错动， 这使得变形急剧增加； 水平方向由于没有围压的约束， 相对于竖向位移值更大； 由于煤柱内的变形主要为水平位移， 使得煤柱表面的部分被顶出而发生剥离破坏， 这与煤柱试验现象完全一致［图3（c）］。

图7、 图8、 图9分别为温度为-10 ℃、 -4 ℃、 20 ℃时煤岩的数值模拟曲线在A、 B、 C三点处的水平位移和竖向位移， 其中A、 B、 C三点的选取与之前提到的选取方法相同。由云图可知， 其变形规律及破坏形式与温度为-15 ℃时煤岩的数值模拟结果基本相同。这说明温度不会影响煤岩在荷载作用下的变形规律和破坏形式。

通过对四种温度下煤岩的单轴压缩数值试验， 可以得到煤岩在不同温度下的弹性模量.由数值试验结果可知， 煤岩的弹性模量随温度的降低而增加。出现这样的现象是因为随着温度的降低， 煤岩中未冻水转化成冰， 冰的胶结作用使得煤岩体克服了裂隙结构的缺陷， 整体性增强， 从而抵抗弹性变形的能力增强； 温度较高的煤岩体中含冰量少， 冰的胶结作用减弱， 从而使的煤岩的整体性下降， 抵抗弹性变形的能力有所减弱。对煤岩的弹性模量进行线性拟合， 可以得到如下结果， 如图10所示。

线性关系可以用如下关系式表示：

式中： E是煤岩弹性模量； E₀是参考弹性模量， 取1 MPa； T是温度； T₀是参考温度取-1 ℃。

当温度为正温时， 煤岩体的弹性模量不再受温度的影响。

通过对四种温度下煤岩的单轴压缩数值试验， 可以得到煤岩在不同温度下的强度值。由数值试验结果可知， 煤岩的强度值随温度的降低而增大。出现这样的现象是因为随着温度的降低， 煤岩中的未冻水冻结成冰， 冰填充了煤岩的裂隙结构， 同时冰的胶结作用也使得在裂隙面相互接触的块体黏结增强， 从而煤岩的强度增加； 随着温度的升高， 煤岩中的冰含量减小， 冰对节理面上相互接触的块体黏结减弱， 在外载作用下， 煤岩中的软弱裂隙面出现缝隙或缝隙增大， 变形迅速发展从而破坏。对煤岩的强度进行拟合， 可以得到如下结果（图11）。

线性关系可以用如下关系式表示：

式中： σ_c 煤岩抗压强度； σ₀是参考强度值， 取1 MPa； T是温度； T₀是参考温度值， 取-1 ℃。

当温度为正温时， 煤岩体的强度不再受温度的影响。

采用本文的细观参数对冻结裂隙煤岩的单轴抗压强度试验进行三维离散元数值模拟， 通过对不同温度下煤岩的应力-应变曲线、 位移以及其力学参数进行分析， 可以得到以下结论：

（1） 在单轴压缩情况下， 由于室内试验煤柱存在天然缺陷， 加载后会发生接触压密过程， 因而试验初始阶段的数值模拟曲线与试验应力-应变关系有一定的偏差， 但随着煤柱的完全压密， 偏差越来越小， 离散元模拟曲线总体上能反映冻结裂隙煤岩的压裂破坏过程。

（2） 在弹性压缩阶段， 煤岩内的位移呈层状分布， 且变形量较小； 屈服后的煤岩将从节理面发生滑移， 煤岩的变形急剧增加， 并且水平位移稍大于竖向位移， 煤柱表面将由表及里逐渐发生剥离破坏。

（3） 随着冻结温度的降低， 煤岩的弹性模量、 强度呈线性规律增加； 但温度变化不会改变煤岩的破坏形式。

本文应用力学试验和数值模拟相结合的研究手段， 揭示了冻结裂隙煤岩细观结构损伤演化与宏观力学变形过程之间的定性关系， 为从微观损伤角度研究岩体工程的安全稳定奠定了理论基础和技术储备。但鉴于自然裂隙岩体的多变性和复杂性， 本文仍然存在以下问题：

（1） 采用了传统的Mohr-Coulomb强度理论来分析， 而Mohr-Coulomb强度理论将剪切滑移面理想为平直面， 只是涉及了煤岩介质在平面情况下的摩擦系数， 忽略了实际情况下岩体断裂面极为粗糙的现象和本质。

（2） 基于煤样CT扫描获得了裂隙/节理几何尺寸和分布状态， 以此采用节理随机生成方法建立了裂隙煤岩的数值模型。但受CT扫描技术和图像处理技术所限， 不能获取煤样真实的裂隙/节理分布状况， 因此， 急需研发高精度、 大尺度岩体三维扫描技术， 以便获取准确的岩体细观结构状态。

（3） 由于冻结裂隙煤岩中冰的存在， 使得煤岩内部的接触、 冰对煤岩的胶结作用等更加复杂， 这都导致了模拟值与真实值之间存在一定的偏差。因此， 需要优化接触模型， 并且还需考虑冻结煤岩中冰的胶结作用， 建立煤岩细观参数与宏观响应之间的关系， 从而精确解释冻结裂隙煤岩的损伤破坏机理及力学性质。