马万里，陈世江

(内蒙古科技大学矿业研究院，内蒙古 包头 014010)

岩石的非均质性向来都是国内外学者研究的重点，不同学者采用不同数值计算软件对非均质岩石的破裂特征进行研究[1-2]。岩石的破裂特征对冲击地压、顶板冒落、巷道支护具有重要意义。

为了描述岩石的非均质性，不同学者采用不同的方式进行刻画，唐春安等[3]通过假设岩石参数服从Weibul方式对岩石弹性模量、抗压抗拉强度进行赋值；田志等[4]通过数字岩石物理的形式对地层沉积过程中的非均质性进行表征；罗荣等[5]通过蒙特卡洛法对岩石非均质性进行描述；陈沙等[6]以图像处理为媒介，利用Matlab灰度识别获取岩石的截面特征，实现了岩石的非均质表征。随着计算机分析手段的进化，CT扫描三维重构技术得到了极大的发展，不同研究方向学者以CT扫描技术为基础，研究了岩石的细观破坏、渗流、结构变化等，获得了较多有益成果，极大发展了岩石力学[7-9]。岩石破裂方面：谢和平等[10]以能量耗散角度对岩石破裂过程不同阶段进行了精细刻画；程昊等[11]采用C++对FLAC进行二次开发，实现了岩石破裂过程中声发射参量导出；朱泽奇等[12]分析了花岗岩破裂及声发射特性，明确了非均质对岩石局部应变产生重要影响；霍丙杰等[13]以FLAC数值计算软件分析了房式采空区煤柱塑性区分布规律，明确了煤柱拉剪复合破坏机制。为了更好地描述岩石每个单元的破碎程度，基于单元损伤、单元安全度、破坏接近度等描述岩石应力状态及三维重构裂隙越来越来受到人们的重视[14-15]，朱万成等[16]以弹性损伤建立了非均质岩石应力-损伤本构模型；王峰[17]以MC准则为基础，定义了以单元应力状态点距包络线距离为基础的破坏接近度概念。

以上研究丰富和发展了岩石非均性表征方法和岩石破裂过程中裂隙演化、能量释放特征，岩石的非均质性同样对水力压裂、爆破等有着重大影响，建立了基于MC准则的应变软化本构模型，通过FLAC内嵌的FISH语言及Matlab编程，实现了岩石的非均质赋参及弹性模量、内聚力、内摩擦角随塑性应变弱化，探讨了不同均质度岩石破裂过程及不同裂隙处理方式的优缺点，为今后非均质爆破试验及数值计算提供基础性工作。

1 非均质岩石应变软化本构模型建立及验证

1.1 岩石非均质性的描述

煤岩体物理力学性质和颗粒分布、缺陷、节理等分布表现出较强的差异性，局部应力、应变导致不同位置应力分布不同，假设其弹性模量、内聚力、内摩擦角、抗拉强度等参数服从weibul分布，且无相关性，可表示为式(1)[18]。

(1)

式中：ε*为微元强度；m为均质度系数，表征岩石的非均质程度，其值越大，表明岩石均质性较高；M为某一力学参数的均值，表征了岩石某一力学性能强度。

1.2 基于MC准则的应变软化本构模型

FLAC3D内置MC本构模型为理想弹塑性模型，即进入塑性阶段后应力随着应变的增加保持不变，很难描述岩石的应变软化过程。其内置的应变软化模型是基于单元的塑性剪切应变对单元的内聚力、内摩擦角、膨胀角的参数进行弱化，对拉破坏单元的适用性较弱，因此，本文在FLAC内置应变软化模型基础上进行修正，建立以塑性应变为参数弱化依据参量，则其内聚力、内摩擦角、弹性模量与塑性应变关系可表示为式(2)。

(2)

式中：c0、φ0、E0和cr、φr、Er分别为初始和残余内聚力、内摩擦角、弹性模量；εp*为初始和残余弹性模量数值之和一半对应的塑性应变。

因此，修正后的MC准则为式(3)[19]。

(3)

式中：Nφ=(1+sinφ)/(1-sinφ)；φ为摩擦角；c为内聚力；σt为抗拉强度；c(εp)为随塑性应变变化的内聚力，εp为0时表示岩石材料处于弹性阶段，εp大于0时表示岩石进入塑性阶段，此时其内聚力、内摩擦角、弹性模量等参数将发生弱化。

若以内聚力的弱化程度对损伤变量进行描述，损伤变量定义为式(4)。

(4)

根据连续介质损伤力学理论，其本构关系可以表示为式(5)[20]。

σ=Eε(1-D)

(5)

结合式(1)～式(5)得到基于MC准则的非均质岩石应变软化(应力-损伤)本构模型。

1.3 模型验证

本文建立的应变软化(应力-损伤)模型为多参数弱化本构模型，和FLAC内置应变软化模型具有一定的相似性，为了验证模型的合理性及正确性，通常需要实验室试验反演获取模型各参数数值，然后将实验室和FLAC数值计算结果进行对比。由于未进行实验室试验，为了验证模型的可行性，本文以单参量轴向应变(F)作为微元强度值衡量岩石损伤程度，并和杨圣奇等[21]试验数据进行对比，对比结果及参数取值见图1。由图1可知，理论模型无法体现岩石压缩过程中的初始压密阶段，弹性阶段匹配性较高，同时屈服破坏阶段应力跌落贴合度较高。岩石单轴压缩表现出较强的弹脆性，均值度m取值较大，为了理论曲线和实验室曲线较为接近，同样可能导致模型参数的失真，但总体上验证了模型的可行性和正确性，表明本文建立模型具有较强的适用性，可以反映岩石全应力-应变整个过程。

图1 实测曲线-理论曲线对比Fig.1 Comparison of measured curve and theoretical curve

2 非均质岩石破裂特征数值计算

2.1 数值计算模型建立及边界条件

为了验证将上文建立的非均质岩石应变软化本构模型引入FLAC，分析岩石非均质性对岩石破裂特征影响及破裂特征描述，本小节建立数值计算模型尺寸为100 mm×50 mm×50 mm立方体模型(图2)，模型一共25 000个网格，上下表明施加位移速度边界5e-7m/step。初始参数见表1，同时计算过程中对弹性单元数量、能量进行监测。

图2 数值计算模型Fig.2 Numerical calculation model

表1 单轴抗压有限差分法模拟数据Table 1 Simulation data of uniaxial compressionfinite difference method

2.2 非均质性及应变软化求解过程

为了表征岩石的非均质性，本文通过Matlab编程生成weibul函数txt文件，数据点和单元网格数量保持一致，通过FLAC读取文件的形式实现对不同参数的对应赋参，同时由于本文不考虑数据点之间的关联性，因此采用多个weibul数据文档进行赋值。具体FISH编程思路为：将模型每个单元力学参数值和weibul函数txt文件数值进行一一对应，同时利用parse功能实现数据的分隔，并采用额外变量显示，均质度为7的内聚力和内摩擦角分布特征见图3。

图3 岩石材料非均质描述Fig.3 Description of heterogeneity of rock material

建立FISH函数对单轴压缩过程中单元内聚力、内摩擦角、弹性模量等参数与塑性应变关系进行描述，计算一步求解每个单元的塑性应变，并根据式(2)对模型各参数进行修改，并监测单元的弹塑性状态、数量及弹性能，再计算一步并不断循环，最终实现非均质岩石单轴压缩破裂全过程数值计算，具体求解过程见图4。

图4 数值计算求解流程Fig.4 Numerical calculation flow

2.3 单轴压缩全应力-应变各阶段

以均质度m=7为例分析单轴压缩过程中能量及裂隙演化全过程(图5)，单轴压缩初期，岩石内部缺陷较大的点即力学参数较弱的点发生破裂，但整体数量较小，呈随机分布，随着加载步数增大，破裂点数量不断增大，岩石仍处于弹性阶段，此时，破裂点开始出现汇集、融合等现象，这些为宏观裂隙的出现提供了必要条件。岩石进入屈服阶段后，开始出现大量拉伸破裂点，一旦进入破坏阶段，大量单元状态由正在剪切破坏状态向过去剪切状态转变，破裂点从随机分布的无序状态向局部剪切带形成的有序状态演变，这和混沌理论较为一致，最终形成“X”型剪切带。岩石弹性能演化和整个应力应变过程贴合度较高(图6)，且在应力峰值时到达最大值，能量演化曲线呈现凹型演化和实验室应力应变曲线较为类似。若以岩石破裂点作为声发射参量，可以看出随着加载步数的增加岩石的总破裂点基本恒定，且声发射参量在岩石峰值时稍靠右侧达到最大值，之后由于岩石大面积卸载，导致声发射参量不断减小，这个过程为岩石内部颗粒重新装配过程，本文称之为平静期(图7)。

图5 岩石破裂过程破裂点分布Fig.5 Distribution of fracture points in the process of rock fracture

图6 应力能量曲线Fig.6 Stress energy curve

图7 声发射曲线Fig.7 Acoustic emission curve

3 讨 论

3.1 不同均质度能量演化及声发射参量

不同均质度条件下岩石应力应变曲线如图8所示。岩石的弹性模量、峰值强度和均质度呈正相关，即岩石的材料参数分布越统一，高强度材料参数单元数量越多，岩石的宏观强度越高，这是均质材料无法体现的。同时随着均质度的增加岩石的脆性程度逐渐增加，表现为峰后阶段应力跌落速度越快，这和内部颗粒的胶结程度、颗粒强度有较大关系。均质

图8 不同均质度应力应变及弹模分布Fig.8 Stress strain and elastic modulus distribution of different homogeneity

度为3、5、7对应的抗压强度分别为9.43 MPa、11.27 MPa、12.21 MPa，和均质度呈对数函数关系。

从不同均质度单元破坏数量及声发射参量可以看出，岩石试样压缩初始阶段，均质度较小的试件已有部分破裂单元，随着均质度的增大，破裂点出现位置对应应变逐渐增加，表明了高均质材料具有较强的稳定性。但从单元破坏总数量上来看，单元总破坏数量和均质度并无直接关联。不同均质度对应的声发射特征差异性较大，均质度较小时，声发射事件数基本呈n型分布，仅在屈服阶段出现小幅度提升；

随着均质度增大，声发射事件数分布形态向倒V型转变(图9)，即峰前声发射事件数逐渐增大，到达峰值时突然大幅度提高；随着均质度增大声发射频率及强度表现出不同的形式，这和前人总结的声发射特征具有高度的一致性，即群震型、前震-主震-余震型和主震型3种典型模式[22]。

图9 不同均质度单元破坏及声发射参量Fig.9 Failure and AE parameters of different homogeneity units

3.2 岩石破裂形态描述

为了更好地表征岩石破裂形态，不同学者以不同的思路进行处理，通常我们采用FLAC自带的塑性区显示方式，认为岩石进入塑性区即进入破坏状态，由之带来的问题则是无法显示岩石的破裂程度，如果仅从塑性区对某关键层位岩层稳定性判定将会带来误差，并且相对保守。

本文则通过对比以塑性区、破坏接近度、损伤值(本文定义)等方式对破裂点的描述，分析各种形式的优缺点，以期获得更合理的表征形式。其中，关于塑性区的定义本文不再赘述，破坏接近度是空间应力的一点沿最不利应力路径到屈服面的距离与相应的最稳定参考点在相同罗德角方向上沿最不利应力路径到屈服面的距离之比，见式(6)。

(6)

式中：I1、J2为应力张量的第一不变量、第二不变量；θσ为罗德角。

本文同样定义了损伤值表达式，已在前文进行说明，不同形式描述岩石破裂特征见图10。由图10可知，三者都可以反映出岩石的破裂特征，其中，塑性区表征程度较为一般，它仅反映了岩石破裂点分布，而损伤值则将剪切带刻画的更为细致，损伤值为1时表示岩石已完全破坏，不同损伤值表明了岩石所处的应力状态；损伤值为0时表示岩石处于弹性状态。破坏接近度则更为细致，它将单元的应力状态统一起来进行体现，从弹性、塑性、破坏状态对每个单元的应力状态进行刻画，破坏接近度小于1时表明岩石处于弹性状态，但其弥补了损伤值的不足，即使单元处于弹性状态，它仍可对较危险单元进行划分；破坏接近度大于2时表示单元处于破裂状态，可以根据其值判定岩石较为危险位置。综上所述，基于破坏接近度的单元接近度评判法则有着较强的优越性。

图10 不同破裂表征形式Fig.10 Different fracture characterization forms

4 结 论

1) 本文建立了非均质岩石应力-损伤本构模型，并采用前人数据验证了模型的正确性。通过Matlab软件和FISH语言实现了FLAC非均质岩石单轴压缩破裂数值计算。

2) 不同均质度单轴压缩数值计算表明：单轴抗压强度和弹性模量和均质度呈对数函数关系，随着均质度系数的增大，岩石由弹塑性向弹脆性演变。

3) 岩石的破裂形态由杂乱无章向局部有序演变，对不同描述破裂特征方式进行了讨论，确定了单元破坏接近度的优越性。