渗透压作用下黄山石灰石矿区炭质泥页岩剪切蠕变特性分析

严 珊，胡 斌，魏二剑，马利遥，田 磊

（1.武汉科技大学资源与环境工程学院，湖北 武汉 430081；2.冶金矿产资源高效利用与造块湖北省重点实验室，湖北 武汉 430081）

我国四川省黄山石灰石矿区大范围存在含缓倾软弱夹层的二叠系石灰岩地层[1]，经过调查发现，这类软弱夹层主要是含炭质的薄层泥页岩，由于我国西南地区雨水丰沛，软弱夹层在雨水浸泡、疏干的长期循环作用下，其蠕变特性将发生显著变化，软弱夹层极易沿着破坏面发生压剪破坏，从而导致矿山边坡滑坡。

目前，国内外学者在岩石蠕变力学相关领域已经进行了大量研究：秦哲等[2]对经历不同饱水-失水循环次数的岩样进行下三轴蠕变试验，从细观角度分析了水岩作用对岩石的损伤机理；王玮玮等[3]对炭质泥岩开展了三轴压缩蠕变试验，探讨了不同干湿循环次数下炭质泥岩孔隙的演化特征和蠕变规律；FUJII 等[4]开展了砂岩在饱水状态与干燥状态下的蠕变试验，发现饱水状态下，砂岩的环向应变更大；李鹏等[5]为探究含水率对软弱夹层剪切蠕变特性的影响，开展了软弱夹层在多种含水率下的剪切蠕变试验；于永江等[6]通过开展不同含水率软岩的直剪蠕变试验，建立了能够描述含水软岩蠕变特性的非线性剪切蠕变模型。上述研究表明，目前对于水岩作用下岩石蠕变特性方面的试验研究主要集中在试验前对岩样进行饱水-干燥循环作用以及岩样的含水率方面，这与实际情况中岩体主要处于压缩状态下持续产生变形的同时，仍不断经历雨水浸泡、疏干的反复循环作用过程存在一定差异。

因此，以黄山石灰石矿区炭质泥页岩为研究对象，开展渗透压作用下的剪切蠕变试验，为了研究其蠕变特性，基于传统的西原模型，提出了一种能同时描述岩石黏弹塑性特性的非线性蠕变模型，并利用该模型进行剪切蠕变参数辨识。以此为基础，运用大型有限元软件FLAC3D对炭质泥页岩的室内剪切蠕变试验进行数值模拟，通过数值模拟方法获得了室内剪切蠕变试验中炭质泥页岩的变形特征，这些结果拓展了试验成果和实测资料，使之可以更合理地评价炭质泥页岩在渗透压作用下的剪切蠕变特性，也为众多含二叠系软弱夹层的露天矿山边坡长期稳定性分析提供一定的理论支撑。

1 矿区工程地质条件

黄山石灰石矿山位于四川盆地边缘低山～低中山地段。矿区地势南高北低，最高标高为1 229 m，最低标高为640 m，相对标高589 m，属低中山区。矿区气候以四川盆地亚热带湿润性季风气候为主，全年阴湿多雨，雨季多集中在6 月到9 月。在大量降雨和周边地震的影响下，雨水入渗流入边坡内部软弱夹层，造成软弱夹层力学参数不断劣化，内部节理裂隙贯通，加之软弱夹层明显的流变效应，极易引起滑坡，如图1 所示的东采区滑坡、老鹰嘴滑坡等，危及矿山生产安全及矿山山脚下村民的生命安全。

图1 黄山石灰石矿区数字化全貌图Fig.1 Digital panorama of Huangshan limestone mining area

2 渗透压作用下炭质泥页岩剪切蠕变试验

2.1 试验概况

室内试验的炭质泥页岩试样采自黄山石灰石矿区，室内试验使用课题组研发的降雨渗流-爆破振动耦合模拟软岩剪切流变试验系统，如图2 所示。将所取岩块送至加工厂加工成150 mm×75 mm×75 mm的长方体试样，试样包含三个等孔距的渗流孔，孔距为37.5 mm，孔深度为37.5 mm，制备完成的试样如图3 所示。

图2 软岩剪切流变试验系统Fig.2 Shear rheological test system of soft rock

图3 软岩试样Fig.3 Soft rock specimen

2.2 试验方案及结果分析

本次室内试验是对四组炭质泥页岩（J-1、J-2、J-3、J-4）进行直剪试验，试验结果见表1。再对另外四组炭质泥页岩（J-5、J-6、J-7、J-8）进行不同渗透水/气压循环次数下的剪切蠕变试验，分别为循环0 次、4 次、8 次、12 次。对试样施加1 MPa 法向压力和五级剪切应力的过程中交替施加0.5 MPa 渗透水/气压。以施加0.5 MPa 渗透水压4 h 后再施加0.5 MPa 渗透气压4 h 为一次循环。需要说明的是，本次室内试验设计的渗透水/气压的循环在施加第四级剪切应力前结束，在施加法向压力和剪切应力的同时进行渗透水/气压的循环是本次试验的一个创新点，更好地还原了实际情况中软弱夹层主要处于压缩状态下持续产生变形的同时，仍不断经历雨水浸泡、疏干的反复循环作用过程。剪切蠕变试验方案设计见表2 和图4。剪切蠕变试验法向压力和渗透水/气压都为恒定值。

表1 软岩瞬时强度参数Table 1 Instantaneous intensity parameters of soft rock单位：MPa

表2 剪切蠕变试验方案设计Table 2 Shear creep test protocol design

图4 剪切蠕变试验方案设计Fig.4 Protocol design of shear creep test

图5 为不同循环次数下的剪切蠕变曲线。由图5 可知，在前四级剪切应力作用下，炭质泥页岩经历了两个阶段，即减速蠕变阶段和稳定蠕变阶段，在第五级剪切应力作用下，炭质泥页岩则经历了完整的蠕变三阶段。在减速蠕变阶段中，炭质泥页岩的剪切应变刚开始随时间增大较快，但其应变率随时间逐渐减小。在稳定蠕变阶段中，曲线近似成直线，炭质泥页岩的剪切应变随时间匀速增加。在加速蠕变阶段中，炭质泥页岩的剪切应变随时间的增长而加速增加，其应变率也随时间迅速增加，直至炭质泥页岩发生破坏。此外，炭质泥页岩岩样在每一级剪切应力加载的瞬间都会产生瞬时应变，且随着渗透水/气压循环次数的增加，瞬时应变也随之增加。在前四级剪切应力作用下，其剪切蠕变变形都非常小。以渗透水/气压循环0 次为例，在施加第一级剪切应力1.78 MPa 时，炭质泥页岩岩样达到稳定蠕变阶段时发生的剪切位移为0.021 mm，在施加第二级剪切应力2.16 MPa 时，炭质泥页岩岩样达到稳定蠕变阶段时发生的剪切位移为0.025 mm，在施加第三级剪切应力2.54 MPa 时，炭质泥页岩岩样达到稳定蠕变阶段时发生的剪切位移为0.029 mm，在施加第四级剪切应力2.91 MPa 时，炭质泥页岩岩样达到稳定蠕变阶段时发生的剪切位移为0.032 mm。

图5 不同循环次数下的剪切蠕变曲线Fig.5 Shear creep curves at different cycle times

2.3 炭质泥页岩的长期强度

岩体的长期强度是评价工程长期稳定性的重要参量[7]，目前应用最广泛的是稳态流变速率法，该方法具有计算简单、直观方便的优点。根据剪切蠕变试验曲线计算不同循环次数下、各级剪切应力下炭质泥页岩岩样的稳态蠕变速率，并采用指数函数拟合得到岩样的稳态蠕变速率-剪切应力关系曲线，如图6 所示。根据该方法确定了在渗透压作用下炭质泥页岩剪切蠕变的长期强度，炭质泥页岩的长期强度即为由稳定蠕变过渡到加速蠕变的极限应力。炭质泥页岩的蠕变破坏强度即炭质泥页岩发生剪切破坏时的强度，由试验即可确定，结果见表3。

表3 软岩长期强度与蠕变破坏强度参数Table 3 Long-term strength and creep failure strength parameters of soft rock

图6 蠕变速率与剪切应力关系曲线Fig.6 Relationship curve between creep rate and shear stress

3 改进的西原模型及参数辨识

3.1 改进的西原模型

传统的西原模型由虎克体、黏弹性体以及黏塑性体串并联而成[8]，如图7 所示。由于传统的西原模型难以描述炭质泥页岩岩样非线性加速蠕变阶段的蠕变特性，在传统西原模型的基础上，提出了一种能同时描述岩石黏弹塑性的非线性蠕变模型，该模型在西原模型的基础上串联了一个非线性黏塑性体，并引入了两个改进的非线性牛顿体黏壶，如图8 所示。E1为瞬时弹性模量，E2、E3均为黏弹性模量，η1、η2、η3均为黏滞系数，ε1、ε2、ε3、ε4为图8 中各部分对应的应变，当剪切应力小于长期强度τs时，该非线性黏塑性体不发挥作用，当剪切应力大于长期强度τs时，岩样进入加速蠕变阶段，该非线性黏塑性体发挥作用[9]。

图7 西原模型Fig.7 Nishihara model

图8 改进的西原模型Fig.8 Improved Nishihara model

夏才初等[10]研究指出直接将蠕变参数中的黏滞系数η替换为η（t）得到的结果将是错误的蠕变方程。因此，本文参照相关文献[11-15]对黏滞系数η2、η3进行等效的非线性化处理，针对图8 中ε3部分和ε4部分非线性黏塑性体构造相应的本构方程，见式（1）和式（2）。

式中：τ0、τs分别为剪切应力和长期强度；ts为岩样由稳定蠕变阶段转为加速蠕变阶段对应的时刻；n为蠕变指数；Δt为单位时间，主要用来保持量纲一致，取1 h。故改进的西原模型一维蠕变方程见式（3）和式（4）。

3.2 剪切蠕变模型参数辨识

根据炭质泥页岩在渗透压作用下的剪切蠕变试验曲线，识别改进的西原剪切蠕变模型相应的参数，其参数识别结果见表4。由表4 可知，模型表达式与试验数据拟合效果很好，相关性系数高，R2均在0.96以上，如图9 所示。由此可知，改进的西原剪切蠕变模型能够很好地模拟炭质泥页岩完整的蠕变三阶段。因此，所建立的改进的西原剪切蠕变模型是合理的。

表4 改进的西原剪切蠕变模型拟合参数Table 4 Fitting parameters of improved Nishihara shear creep model

图9 试验曲线与拟合曲线对比图Fig.9 Comparison of test curve and fitting curve

采用指数函数拟合得到炭质泥页岩的弹性模量-渗透水/气压循环次数关系曲线与黏滞系数-渗透水/气压循环次数关系曲线如图10 所示。受渗透水/气压循环次数的影响，炭质泥页岩的平均瞬时弹性模量、平均弹性模量、平均黏滞系数均不断减小，进一步说明了降雨反复作用影响下炭质泥页岩力学强度降低，水理性质变差，从而影响矿山边坡稳定性。

图10 模型参数变化规律Fig.10 Variation pattern of model parameters

4 数值模拟分析

4.1 数值计算模型

将自定义的蠕变本构模型进行二次开发，数值计算模型的形状采用长方体，尺寸与炭质泥页岩试样的尺寸相同，长宽高分别为150 mm×75 mm×75 mm，共8 803 个单元，3 385 个节点。采用的边界条件为试样底部约束法向位移，试样上半部分的左面约束横向位移。采用的加载条件与室内试验中渗透水/气压循环4 次下的加载条件相同，试样上表面施加1 MPa 的法向应力，试样下半部分的右面施加剪切应力，剪切应力为1.69～3.20 MPa，应力梯度为0.38 MPa，分五级加载，每级剪切应力加载时间为24 h，共120 h。采用自定义的蠕变本构模型计算，计算参数见表4。

4.2 数值模拟结果分析

图11 为数值模拟计算结果图。由图11 可知，最大位移值位于试样的右下角部分，由于试样上半部分的左面施加了横向约束，所以试样的左上角部分位移值为0，且随着剪切应力的不断增加，X轴方向的位移也在不断增大。在第一级剪切应力作用下，最大位移值为0.18 mm；在第二级剪切应力作用下，最大位移值为0.32 mm；在第三级剪切应力作用下，最大位移值为0.47 mm；在第四级剪切应力作用下，最大位移值为0.62 mm；在第五级剪切应力作用下，最大位移值为1.17 mm。X轴方向的位移云图的变化趋势与实际情况也比较相符，图11（b）可以看到明显的剪切变形。

图11 数值模拟计算结果Fig.11 Numerical simulation calculation results

对试样上（0.15，0，0.037 5）点进行位移监测，得到如图12（a）所示的位移监测图。由于在试验过程中采用分级加载方式，剪切蠕变呈台阶型，需利用Boltzmann 迭加原理对试验数据进行二次处理，处理后的剪切蠕变曲线如图12（b）所示。由此可见，使用FLAC3D数值模拟软件得到的剪切蠕变曲线与室内试验剪切蠕变曲线总体趋势相同，从而验证了室内试验和所建立的非线性黏弹塑性本构模型的正确性。

图12 室内试验结果与数值模拟结果对比Fig.12 Comparison of indoor test results and numerical simulation results

5 结论

1）在渗透水/气压的循环作用下，随着剪切应力的增加，炭质泥页岩试样的剪切蠕变曲线呈现了完整的蠕变三阶段。不同的渗透水/气压循环次数下，试样的含水率不同，剪切蠕变曲线也会有所区别。

2）基于室内渗透压作用下炭质泥页岩剪切蠕变试验得到的剪切蠕变曲线，在西原模型的基础上，串联了一个非线性黏塑性体，并引入了两个改进的非线性牛顿体黏壶，建立了改进的西原剪切蠕变本构模型，能够很好地拟合炭质泥页岩试样蠕变的三个阶段，具有较好的实用性。

3）运用FLAC3D软件对炭质泥页岩试样以相同的室内试验条件进行数值模拟计算，结果表明：最大位移值位于炭质泥页岩试样的右下角部分，可以看到明显的剪切变形。使用FLAC3D数值模拟软件得到的X轴方向的位移监测曲线与室内试验剪切蠕变曲线趋势相同，都呈台阶状且包含完整的蠕变三阶段，与实际结果比较吻合。