张学朋 ，王 刚 ，蒋宇静，吴学震 ，王者超，黄 娜

（1.山东科技大学 矿山灾害预防控制省部共建国家重点实验室培育基地，山东 青岛 266590；2.山东科技大学山东省土木工程防灾减灾重点实验室，山东 青岛 266590；3.山东大学 岩土与结构工程研究中心，济南 250061）

1 引 言

在外部荷载的作用下岩石内部的原生裂隙不断演化，新生裂隙会持续萌生和发展，岩石不断产生损伤，最终形成宏观裂缝，导致岩体失稳破坏[1]，这一破裂演化规律是揭示岩石破坏机制、探究岩体本构关系、解决工程稳定性的基础[2]。室内试验和现场原位测试只能观察到岩石的宏观破裂，无法观察到细观结构上微裂纹的张开、扩展及贯通等现象。基于颗粒流理论的PFC 程序可以将细观物理力学性质参数与宏观力学参数联系起来模拟岩石的基本力学特性，并能实时监测试验中岩石内微裂纹的萌生、扩展及贯通过程，实现岩石在受压致裂过程中宏细观过程和性质的统一。

目前应用颗粒流方法进行岩石力学特性细观模拟的研究主要集中在参数的校核及微裂纹的破裂机制研究[3-7]，而岩石变形破坏过程是能量的复杂转化过程，实质上是能量耗散与能量释放的综合结果。众多学者[8-11]对岩石破坏过程的能量变化规律进行了大量试验研究，并取得了很多有价值的研究成果。但室内试验难以将能量的变化过程与岩石微观破裂机制有效地联系在一起，有关两方面相联系的离散元数值模拟研究还鲜见报道。

本文基于颗粒流理论，借助PFC2D的黏结颗粒模型（BPM）建立双轴压缩模型，通过“试错法”得出与室内试验宏观力学参数相适应的细观物理力学性质参数，以此为依据模拟岩石压缩试验中微裂纹破裂演化规律以及能量的变化规律，深入研究岩石的宏细观破裂机制。

2 花岗岩细观物理力学性质参数

2.1 双轴压缩试验的实现

进行颗粒流模拟试验时，首先生成试样模型墙体（见图1）将颗粒包围，然后通过移动墙体模拟施加围压和加载过程，给定模型顶、底部墙体的速度来模拟应变控制加载方式，侧向墙体的运动由伺服控制程序自动控制，在整个试验过程中围压保持不变。记录整个试验过程中墙体的位移、不平衡力，通过后处理得到试件的宏观变形过程数据。

图1 花岗岩试件及PFC 模型Fig.1 PFC model and granite specimen

颗粒流模型尺寸与室内试验试件尺寸一致，选用φ 50 mm×100 mm，花岗岩试件及PFC 模型，见图1。计算模型中，最小颗粒半径取0.3 mm，粒径比选为1.66[4]。

2.2 细观物理力学性质参数初始值确定方法

模型颗粒的细观参数主要有颗粒接触模量 Ec、颗粒法向刚度与切向刚度比值 kn/ks、摩擦系数f，平行黏结模型的细观参数有平行黏结半径乘子λ、黏结模量、黏结法向刚度与切向刚度比值 kn/ks和法向和切向黏结强度、。

首先，通过室内试验确定出材料的宏观力学参数，即弹性模量E、泊松比ν 和抗压强度值σc及抗剪强度参数c、φ。本文中选用的是三轴压缩试验结果，选取抗压强度值为岩石三轴抗压强度。通过对宏观力学参数进行分析，初步确定颗粒接触模量Ec和平行黏结模量，这里的接触模量不同于宏观的弹性模量，往往比宏观弹性模量大[7,12]。

颗粒刚度初始值为

式中：R[A]、R[B]为两接触颗粒的半径。

平行黏结初始法向及切向刚度值为

2.3 花岗岩细观物理力学性质参数确定

对于模拟岩石材料，需要设置黏结模型来表征颗粒之间胶结物的存在，黏结模型分为接触黏结模型和平行黏结模型。使用接触黏结模型时，由于黏结断裂后只要颗粒仍保持接触，接触刚度仍然有效，因此接触黏结的断裂不会对宏观刚度产生太大影响，这与岩石的破裂机制不符。使用平行黏结模型时，宏观刚度由接触刚度和黏结刚度组成，黏结破裂会立即导致宏观刚度的下降，因此，平行黏结通过在拉伸，或剪切断裂时刚度相应的降低可更逼真地模拟岩石类材料。本文选用平行黏结模型。

由于PFC 模拟采用细观物理力学性质参数表征颗粒及黏结的力学性质，且这些细观参数无法从室内试验直接获取，因此，在数值模型进行计算之前需要对模型的细观物理力学性质参数进行标定。在此过程中，需要进行一系列与室内试验或现场条件类似的模型试验，并将模拟结果与室内试验或原位测试结果进行对比，采用“试错法”反复改变细观参数[12]，直到模型的宏观力学响应满足要求为止。其基本过程见图2。

图2 细观物理力学性质参数校准过程Fig.2 Calibration process of mesoscopic physico-mechanical parameters

本文选用黄岛国家石油储备库地下水封石油洞库花岗岩的力学参数，室内三轴压缩试验采用长春朝阳试验仪器有限公司研制的TAW-2000微机控制电液伺服岩石三轴试验机。颗粒流程序中通过“试错法”反复调整对比，使得颗粒流模型的宏观力学参数接近花岗岩室内试验力学参数，最终确定细观物理力学性质参数，见表1。校核后，花岗岩室内试验与BMP 模型模拟宏观弹性力学参数吻合较好，弹性模量和泊松比分别为28.7、28.4 MPa 和0.230 0、0.228 5。

表1 花岗岩细观物理力学性质参数Table 1 Mesoscopic physico-mechanical parameters of granite

3 BPM 模型模拟与室内试验宏观参数校核对比

将标准校正的细观物理力学性质参数用于颗粒流模型中，在围压3、6、10 MPa 作用下分别进行了花岗岩压缩试验模拟。对比图3 中BPM 模型及花岗岩室内试验应力-应变曲线可以看出，随着围压的增加，试件的屈服应力及峰值强度都有所增加，且BPM 模型与花岗岩试验结果吻合较好。

图3 不同围压模型应力-应变曲线Fig.3 Stress-strain curves at different confining stresses

图4为花岗岩试件与BPM 模型在不同围压作用下的破坏形态。从模拟结果与试验结果可以看出，由于围压较小，试件主要表现为单斜面剪切破坏，且两者吻合较好。

峰值强度包络线反映出峰值强度与围压的变化关系[13]，可以用摩尔-库仑准则来表征岩石的峰值强度包络线，采用以下形式：

式中：σ1为最大主应力；σ3为最小主应力；c为黏聚力；φ为内摩擦角。

图4 花岗岩试件及BPM 模型破坏形态Fig.4 Failure patterns of specimens of granite and BPM model

表2为在不同围压作用下花岗岩及BPM 模型的峰值强度值。应用式（6）处理得到花岗岩及BPM模型的黏聚力和内摩擦角如图5 所示。对比分析可知，BPM 模型的黏聚力比花岗岩的大，而内摩擦角比花岗岩的小。这一差异主要是BPM 模型采用圆形颗粒造成的，Potyondy 等[4]、Cho 等[14]曾指出BPM 模型存在这一问题。

表2 不同围压作用下花岗岩与BPM 模型的峰值强度Table 2 Peak strength of granite and BPM model at different confining stresses

图5 摩尔-库仑强度包络线Fig.5 Envelope of Mohr-Coulomb strength

从宏观力学参数、应力-应变曲线以及试件破坏形态对比来看，BPM 模型模拟结果与室内试验结果具有良好的一致性。

4 花岗岩微裂纹破裂演化及能量变化规律

以围压为10 MPa 的BPM 模型双轴压缩为例，通过分析微裂纹破裂演化及能量变化规律对花岗岩的破坏机制进行更加深入的研究。

4.1 花岗岩压缩过程中微裂纹破裂演化规律

BPM 模型在受力过程中，颗粒间的黏结强度小于颗粒间所传递的强度时，颗粒黏结就会发生断裂，即对应于岩石内部的微裂纹[15]。裂纹扩展过程中应变能将以弹性波的方式快速释放，产生次声波、声波或超声波，即声发射现象(AE)。一般声发射过程线如图6 所示[16]。图7为BPM 模型的轴向应力-应变曲线、裂纹数量变化曲线及轴向应力与体积应变的关系曲线。图8为声发射频度曲线。为研究模型内部裂纹分布演化规律，在应力-应变曲线中设置了5个监测点A～E，分别对应于轴向应变0.0%，0.169%，0.388%，0.507%，0.939%。图9为各监测点对应试件的裂纹分布。

图6 一般声发射过程曲线Fig.6 Curves of general acoustic emission process

图7 BPM 模型轴向应力-应变曲线及力学特性Fig.7 Axial stress-strain curves and mechanical behavior of BPM model

图8 BPM 模型轴向应力-应变曲线及声发射频度曲线Fig.8 Curves of axial stress-strain and acoustic emission frequency of BPM model

图9 BPM 模型在变形破坏过程中微裂纹发育状态Fig.9 Meso-crack development state of BPM model in process of deformation and failure

由图7～9 可以看出，受到外部荷载作用后试件在B 测点（34%σc，宏观上的启裂强度，σc为峰值强度）开始产生微裂纹，发出第一次声发射，在声发射过程线上定义为初始发射点a（见图6 及图8中a 点）。从初始发射到应力屈服点C（78.6%σc），声发射事件一直比较稀疏，振幅随着裂纹的张开、扩展有所增大，主要在0～8 次变化范围内，表明试件由弹性阶段向稳定破裂传播阶段过渡。此阶段体积应变表现出线性变化且体积随荷载增大而逐渐减小，试件呈现受压缩状态。

当荷载超过屈服应力点C 后，试件从稳定破裂传播阶段转向非稳定破裂传播阶段，其应力-应变曲线呈下凹型，变形增大，体积应变曲线经过短暂的不变阶段后出现反弯，开始发生体积膨胀，产生不可逆的塑性变形，裂纹进一步增加、扩展，声发射事件呈现快速增长趋势，其值增加到57 次，说明试件已经达到临界失稳点b（图6 及图8 上b 点）。超过b 点以后，声发射事件反而会有所衰减，主要是由于此时裂纹扩展产生的新自由面急剧增加而消耗部分能量，在声发射过程线上则呈现出能量的吸收谷，其底点定义为聚能点（图6 及图8 上c 点）。待裂纹扩展暂时平衡之后，随着荷载增加，裂纹重新开始压密，能量开始再次积聚。这一阶段体积膨胀的主要原因是：在不断加载的过程中，由于颗粒间破裂的发生，微裂隙的增加、扩展、贯通等现象出现，使试件内部孔隙不断地增加，促使其在宏观上表现为体积增大。

轴向应力超过峰值D 点（σc）以后，试件表现出脆性发展特征，轴向应力大幅度下降，裂纹数量急剧增加，微裂纹进一步扩展为宏观裂缝，声发射事件出现最大值92 次（图6 及图8 中d 点）。随后轴向应力曲线平稳发展，裂纹数量保持恒定的速率增长，声发射事件降至较低水平，变化范围主要是在0～10 次，表明试件又恢复了稳定破裂的发展态势。此阶段试件破坏主要是沿着宏观裂缝产生摩擦滑动，最终形成宏观断裂面（如图9 中E 点）。

4.2 花岗岩压缩过程中能量耗散规律

图10为BPM 模型在双轴压缩过程中能量耗散曲线。由图可见，屈服应力（图8 中C 点）之前，摩擦能、动能、应变能和黏结能的和基本与边界能相等。黏结能和应变能所占比例较大，这部分能量与裂纹的产生和驱动有关，其消长与模型材料的劣化有关。代表裂纹作用的摩擦能则与之相反，它们之间是此消彼长的关系。因为裂纹产生要克服黏结能，然后在应变能的驱动下扩展，裂纹产生以后摩擦部分才开始起作用。在试件达到峰值强度后，黏结能和应变能急剧减小，摩擦能急剧增加，摩擦能所占比例随裂纹进一步扩展逐步提高，由此可以看出，摩擦作用是残余强度的主要提供者。试件整个变形过程中动能所占比例不大，与加载过程及试件内部动态平衡有关，说明试件变形不是很剧烈，裂纹稳定扩展贯通。

图10 BPM 模型破裂过程中能量耗散曲线Fig.10 Energy dissipation laws during BPM model fracture process

通过对花岗岩试件微裂纹演化规律及能量耗散规律的研究，从能量的角度揭示了其受压变形破坏的细观力学机制：在轴压和围压的共同作用下，花岗岩内部能量开始发生积聚，颗粒之间首先克服黏结能产生微裂纹，然后在应变能驱动下微裂纹进一步扩展，磨擦作用渐渐开始起作用且其比例随着裂纹的进一步扩展逐步提高。达到峰值强度以后应变能和黏结能急剧释放，摩擦能随着裂纹的扩展、贯通急剧增加且居主导地位。残余阶段黏结能在边界能和应变能的驱动下多转化为其他能量形式，颗粒间很难再克服黏结强度产生大量新的裂纹，岩石最终的失效形式主要是沿宏观裂缝滑动破坏。

5 结 论

（1）微裂纹破裂演化规律：岩石压缩初始阶段（图8 中AB 阶段）试件内部不存在裂纹，达到宏观启裂强度B 点以后，开始产生裂纹，并发出第一次声发射。此后，裂纹的扩展主要经历平稳发展、急剧增加、平稳发展3个阶段，表明试件主要经历弹性阶段，稳定破裂传播阶段，非稳定破裂传播阶段以及残余阶段的稳定破裂传播阶段。

（2）花岗岩在受力变形破坏的过程中积聚的能量以声发射的形式得到释放：初始发射以后，声发射稀疏发展，在0～8 次的变化范围内；屈服应力C点过后，出现小幅度增加，增加到57 次，说明试件达到临界失稳点；随后由于裂纹扩展产生的新自由面急剧增加而消耗部分能量，从而试件进入能量的吸收谷；能量积聚后声发射急剧增加，峰值应力后达到最大值92 次，试件开始失效破坏。

（3）能量耗散规律：屈服应力之前，摩擦能、动能、应变能和黏结能的和基本与边界能相等。黏结能和应变能所占比例较大，这部分能量与裂纹的产生和驱动有关，其消长与模型材料的劣化有关。代表裂纹作用的摩擦能则与之相反，它们之间是此消彼长的关系。在试件达到峰值强度后，黏结能和应变能急剧减小，摩擦能急剧增加，摩擦能所占比例随裂纹进一步扩展、贯通逐步提高且居主导地位。试件整个变形过程中动能所占比例不大，与加载过程及试件内部动态平衡有关。

（4）从能量的角度揭示了其受压变形破坏的细观力学机制：在轴压和围压的共同作用下，花岗岩内部能量开始发生积聚，颗粒之间首先克服黏结能产生微裂纹，然后在应变能驱动下微裂纹进一步扩展，磨擦作用渐渐开始起作用且其比例随着裂纹的进一步扩展逐步提高。达到峰值强度以后，应变能和黏结能急剧释放，摩擦能随着裂纹的扩展、贯通急剧增加且居主导地位。残余阶段黏结能在边界能和应变能的驱动下多转化为其他能量形式，颗粒间很难再克服黏结强度产生大量新的裂纹，岩石最终的失效形式主要是沿宏观裂缝滑动破坏。

[1]苏承东,张振华.大理岩三轴压缩的塑性变形与能量特征分析[J].岩石力学与工程学报,2008,27(2):273-280.SU Cheng-dong,ZHANG Zhen-hua.Analysis of plastic deformation and energy property of marble under pseudo-triaxial compression[J].Chinese Journal of Rock Mechanics and Engineering,2008,27(2):273-280.

[2]蔡美峰,何满潮,刘东燕.岩石力学与工程[M].北京:科学出版社,2002:157-166.

[3]POTYONDY D O.A bonded-disk model for rock:relating micro properties and macro properties[C]//Proceedings of the Third International Conference.Santa Fe:[s.n.],2002:340-345.

[4]POTYONDY D O,CUNDALL P A.A bonded-particle model for rock[J].International Journal of Rock Mechanics and Mining Sciences,2004,41:1329-1364.

[5]FAKHIMI A,CARVALHO F,ISHIDA T,et al.Simulation of failure around a circular opening in rock[J].International Journal of Rock Mechanics and Mining Sciences,2002,39(4):507-515.

[6]STAVROPOULOU M.Modeling of small-diameter rotary drilling tests on marbles[J].International Journal of Rock Mechanics and Mining Sciences,2006,43(7):1034-1051.

[7]徐小敏,凌道盛,陈云敏,等.基于线性接触模型的颗粒材料细-宏观弹性常数相关关系研究[J].岩土工程学报,2010,32(7):991-998.XU Xiao-min,LING Dao-sheng,CHEN Yun-min,et al.Correlation of microscopic and macroscopic elastic constants of granular materials based on linear contact model[J].Chinese Journal of Geotechnical Engineering,2010,32(7):991-998.

[8]谢和平,彭瑞东,鞠杨,等.岩石破坏的能量分析初探[J].岩石力学与工程学报,2005,24(15):2604-2608.XIE He-ping,PENG Rui-dong,JU Yang,et a1.On energy analysis of rock failure[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(15):2604-2608.

[9]刘新荣,郭建强,王军保,等.基于能量原理盐岩的强度与破坏准则[J].岩土力学,2013,34(2):305-310.LIU Xin-rong,GUO Jian-qiang,WANG Bao-jun,et al.Investigation on mechanical properties and failure criterion of salt rock based on energy principles[J].Rock and Soil Mechanics,2013,34(2):305-310.

[10]杨圣奇,徐卫亚,苏承东.大理岩三轴压缩变形破坏与能量特征研究[J].工程力学,2007,24(1):136-141.YANG Sheng-qi,XU Wei-ya,SU Cheng-dong.Study on deformation failure and energy properties of marble spe-cimen under triaxial compression[J].Engineering Mechanics,2007,24(1):136-141.

[11]许国安,牛双建,靖洪文,等.砂岩加卸载条件下能耗特征试验研究[J].岩土力学,2011,32(12):3611-3617.XU Guo-an,NIU Shuang-jian,JING Hong-wen,et al.Experimental study of energy features of sandstone under loading and unloading[J].Rock and Soil Mechanics,2011,32(12):3611-3617.

[12]Itasca Consulting Group Inc.Manual of particle flow code in 2-dimension(Version 3.10)[M].Minneapolis:Itasca Consulting Group Inc,2004.

[13]余华中,阮怀宁,褚卫江.大理岩脆-延-塑性转换特性的细观模拟研究[J].岩石力学与工程学报,2013,32(1):55-64.YU Hua-zhong,RUAN Huai-ning,CHU Wei-jiang.Mesoscopic simulation study of brittle-ductile-plastic transition character of marble[J].Chinese Journal of Rock Mechanics and Engineering,2013,32(1):55-64.

[14]CHO N,MARTIN C D,SEGO D C.A clumped particle model for rock[J].International Journal of Rock Mechanics &Mining Sciences,2007,44(7):997-1010.

[15]HAZZARD J F,YOUNG R P,MAXWELL S C.Micromechanical modeling of cracking and failure in brittle rocks[J].Journal of Geophysical Research(Solid Earth),2000,105(7):1978-2012

[16]龙湘桂,汪瑞强.岩石力学试验中的声发射测试[J].岩土工程学报,1981,3(2):69-74.LONG Xiang-gui,WANG Rui-qiang.Acoustic emission experiments of rock mchanics[J].Chinese Journal of Geotechnical Engineering,1981,3(2):69-74.