花岗岩颗粒流模型循环压缩作用下能量特征分析*

龙恩林，陈俊智

(昆明理工大学 国土资源工程学院，云南 昆明 650093)

0 引言

金属矿产资源作为世界各国重要的工业原料，对国民经济发展和科技进步有着举足轻重的作用，就我国而言，矿产资源经过多年的开采，位于浅部的资源逐年减少或已枯竭，因此我国金属矿山资源开采正处于向深部全面推进的阶段，进入深部开采后，岩体的结构和力学性质会发生重大变化，浅部的硬岩可能会变为软岩，弹性体可能变为塑性体，尤其是在深部高地应力作用下的硐室开挖、采场爆破、巷道支护等施工性周期荷载扰动作用会进一步扩大，从而使得岩体力学性质更加难以把控，而能量转化作为岩体各种物理变化的本质，研究能量变化对解决工程实际问题有着重要意义[1]。

因此，有许多学者对岩石在受力过程中的能量变化进行了大量研究，如邓华锋等[2]研究了砂岩在单轴循环加卸载作用下应力-应变滞后现象对能量参数计算的影响并提出了能量参数的修正计算方法；何明明等[3]建立了单轴循环荷载下砂岩的耗散能随循环次数变化的能量演化方程；张楚旋等[4]探讨了不同应力路径下花岗岩卸荷破坏过程的能量耗散特征，认为能量耗散与时间呈非线性关系，且随着施加围压的增大，相同变形情况下能量耗散也越大；张志镇等[5]分析了红砂岩在加载过程中的能量分配规律，并认为岩石的弹性能密度-应力曲线更能够反映材料的力学性质，表现出比应力-应变曲线的优越性；苏国韶等[6]利用自主开发的真三轴岩爆试验系统实现了岩爆碎块耗能的定量分析，得出不同加载速率下，岩爆碎块耗能均以剪切耗能为主的结论；ZHANG等[7]使用SEM对不同加载速率下大理岩的断口形貌进行了分析，结果表明岩石动态断裂时的能量明显大于静态时的；JIANG等[8]对煤样进行了三轴加卸载试验，认为累积耗散能量呈指数函数形式增长；SUN等[9]研究了砂岩三轴循环加卸荷下的损伤演化，提出了1种基于等增幅多重循环加卸荷试验的岩石损伤演化模型。以上能量研究主要基于宏观力学参数(应力-应变曲线)积分计算求得，而在积分求解时曲线的滞后效应等必然导致结果的误差，因此有学者从细观力学角度出发，采用数值模拟软件定量计算出试件在加卸载过程中的能量大小，如王云飞等[10]结合室内单轴试验和PFC颗粒流数值模拟软件获得了花岗岩的细观力学参数，研究了不同围压下的花岗岩能量演化机制；温韬[11]通过PFC数值模拟对锦屏砂岩进行了循环加卸载数值试验，分析了其在破坏过程中摩擦能、动能以及应变能等细观能量的演化规律；陈子全[12]用PFC数值模拟软件的破裂源机理对岩石的破裂及其能量机理进行了分析；姜耀东等[13]运用FLAC数值模拟软件结合岩石CT扫描结果，进行了多载速率下单轴数值模拟，得出了加载速率增加将会使能量耗散和积聚更早发生的结论。

就以上研究成果而言，对常规受力下岩石破坏时的能量转化研究较多，但对花岗岩在周期荷载作用下应力-应变曲线不同变形阶段的能量变化机理研究较少，而花岗岩作为地下开采所面临的普遍岩体，因此，本文在室内单轴试验的基础上，结合PFC颗粒流数值模拟软件对各阶段的能量变化(总能量、弹性应变能和耗散能)进行详细研究，其结果对深入了解花岗岩在周期荷载下的力学特性和能量破坏机制有一定参考意义，同时也为矿山对施工扰动下矿柱的安全性分析给予一定帮助。

1 花岗岩室内试验与数值模拟

1.1 花岗岩室内试验

本次试验试件取自云南某矿山，将其按照工程岩体试验方法标准制备成100 mm×50 mm×50 mm的长方体试件，选取其中完整性较好地试样进行单轴压缩试验，从而得到相关力学参数以便进行后续数值试验。室内试验系统采用TAW-2000岩石三轴试验机，加载时采用应力控制，加载速率为0.5 kN/s。室内试样宏观力学参数：平均强度为75.2 MPa，弹性模量为32.31 GPa，泊松比为0.16。

1.2 数值模型与试验

颗粒间的相互作用以颗粒流接触本构模型来表征，在PFC颗粒流数值模拟软件内置的众多接触本构模型当中，平行黏结模型因其在拉伸或剪切断裂时可以有效模拟岩石类材料而被广泛应用[14-16]。因此本文采用平行黏结模型来表征花岗岩颗粒间的胶结物，数值试验计算模型与室内试验试件尺寸保持一致，在结合JENSEN等[17]学者关于试件尺寸D和内部颗粒平均粒径d50的比值大于30～40时，颗粒数量的多少对模型结果不会产生显著影响的前提下，通过试错法多次调节细观参数，以数值模型求解得到的应力-应变曲线与室内试验曲线相吻合为准则完成了模型的标定，其标定结果对比如图1所示，初始数值模型与构架如图2所示。

图1 室内试验和数值试验应力-应变对比Fig.1 Comparison of stress-strain diagrams for laboratory tests and numerical tests

图2 初始数值模型与构架Fig.2 Initial numerical model and architecture

数值试验方案为常规单轴压缩和循环加卸载，循环应力水平约为峰值应力的20%～40%，循环次数为1，2，5，10，25，50，100，200以及300，尽可能模拟实际工程中的循环扰动，从而得到逼近于实际的结论。

2 能量特征分析

2.1 能量原理

能量转化作为材料各种参数变化的本质特征，贯穿于岩石变形破坏的始终，谢和平等[18]根据热力学第一定律，假设岩体在外力作用下产生变形且该过程中没有与外界进行热交换，则有：

U=Ud+Ue

(1)

式中：U为外力对岩石所做的总功；Ud为受力过程中岩石所耗散的能量；Ue为受力过程中岩石所存储的弹性应变能。

图3为可释放弹性应变能与耗散能对应关系，其中阴影部分为岩石在受力过程中存储于自身的弹性应变能，Ud为卸荷弹性模量Eu与应力-应变曲线所围成的面积。

图3 可释放应变能与耗散能对应关系Fig.3 Corresponding relationship between releasable strain energy and dissipated energy

2.2 能量转化分析

岩石在受力致破坏过程中涉及到众多能量转化，如弹性变形对应的弹性势能、塑性变形对应的塑性势能、裂纹生成与扩展所需的表面能、岩石破坏后产生动能、辐射能以及其他各种尚未发现的能量形式等等。现采用PFC数值模拟软件对岩石破坏过程中主要的能量转化(总能量、弹性应变能与耗散能)进行分析。

图4为2次循环下各能量与裂纹随应力-应变的变化曲线，将曲线分为4个阶段：

图4 能量与裂纹随应力-应变的变化曲线Fig.4 Curves of energy and crack change with stress-strain

1)初始压密阶段(OA)：PFC数值模拟软件在初始建模时就对模型内部颗粒进行了重生压密，颗粒间接触相对均匀，因此应力-应变曲线偏向于线性化增长且压密不明显，就该阶段能量而言，外界对试件所做的总功基本都转化为弹性应变能，没有裂纹产生，能量耗散基本为零。

2)弹性阶段(AB)：A点处总功与弹性应变能开始分离，产生的少量耗散能主要用于循环时颗粒间的摩擦，粒间黏结键的少量断裂，微裂纹的萌生等，因此耗能较小，试件宏观上表现为线弹性；循环结束后随着加载的进行，弹性应变能存储速率和能量的耗散速率逐渐加快，且二者皆呈现出线性稳定增长的趋势。

3)屈服阶段(BC)：试件进入屈服阶段后，虽然其能量转化仍以弹性应变能的存储为主，但耗散能不断增加，内部损伤的累积开始，并于B点处产生少量可见裂纹，随着裂纹数量的不断产生，耗散能增长速率也不断增大，呈现出“下凸”型增长；相反，弹性应变能的存储开始减缓，且其增长速率于峰值点C附近减小为0。总体而言该阶段弹性应变能增速开始变缓，耗散能增速不断增加，可见裂纹开始显现并迅速增加。

4)峰后阶段(CD)：弹性应变能在峰值点C附近达到储能极限，前期存储的弹性应变能开始急剧释放，主要用于内部新裂纹的产生以及宏观破裂面的贯通形成，裂纹数量也于该阶段激增直至试件全面破裂而停止；而相应的耗散能也因试件损伤的加剧开始迅速增长，并处于该阶段能量占比的主导地位(46.8%～60.4%)。

从以上可以看出，在岩石受力致失稳破坏的整个过程中，能量急剧释放是导致岩石失稳破坏的真正原因，峰值点C之前主要是弹性应变能的存储和能量的缓慢耗散；峰值点C之后则是弹性应变能的急剧释放，耗散能的快速增加，裂纹的不断增长。

2.3 耗散能特征分析

岩石能量的耗散是岩石失稳破坏的本质属性，反映的是岩石内部微观缺陷(裂纹、孔洞)的产生、扩展和贯通，材料强度不断弱化并丧失的过程，有学者以循环过程中形成的滞回环面积大小来描述耗散能[2]。图5为循环荷载作用下单次循环耗散能数值与循环次数的关系曲线，由图5可知，颗粒流数值模型在循环荷载过程中耗散能的发展随循环次数大致呈现出线性变化，该变化过程可划分为2个阶段，即：初期等速发展阶段(AB)以及后期加速发展阶段(BC)。其中初始阶段，耗散能随着循环次数等速稳定增长，这主要是由于颗粒流数值模型密实度较高，内部微观缺陷较少，压密阶段消耗的能量较少，每循环1次就形成1个滞回环，耗散能稳定缓慢增长；而后期加速发展阶段(BC)，在循环次数达到一定量级时，多次循环下岩石内部累积的损伤开始显现，能量主要消耗在颗粒间的摩擦、内部微小裂纹的萌生和生成新的塑性区上，每多循环1次，内部损伤便得到累积，因此在高循环次数下，循环耗散能数值于后期更趋向于加速增长的状态。

图5 耗散能随循环次数变化的关系曲线Fig.5 Relationship curves between dissipated energy and number of cycles

图6为峰值强度之前耗散能随轴向应力变化的关系图。从图6可知，不同循环次数下耗散能随轴向应力变化的散点图几乎重合，这主要是数值模型的同一性所致；从速率上看，散点图的拟合曲线呈现为非线性增长，先慢后快，初期主要是由于颗粒流模型内部接触良好，初始缺陷耗能较少所致，对应于图4中OA段；随后试件进入弹性阶段，对应于图4中AB段，虽然此时试件内部被进一步压实，但是由于应力的加大和集中，内部依旧会有微小裂纹的发育扩展等行为，导致耗能速率慢慢增大；当应力达到60 MPa附近时，增速达到最大，此后其能量的消耗主要用于裂纹的发展，对应于图4中BC段。

图6 耗散能随应力变化的关系曲线Fig.6 Relation curve of dissipated energy with stress

2.4 耗散能的变化方程

材料的非弹性性质在任何荷载形式下都存在，即使是在低应力的条件下依旧不是完全按照弹性性质进行力学响应，在循环荷载作用下加卸载曲线路径不重合，每次加卸载都形成滞回环，这主要由能量的损耗所致，有学者将其定义为材料的阻尼，并进行大量的试验得出材料的阻尼耗能与应力呈现对数关系[3]。

从图6拟合结果中可以发现，材料的耗散能与应力呈现出幂函数关系，即符合上述阻尼理论：

Ud=Kσn

(2)

式中：K，n分别为岩石相关力学性质和循环次数与应力有关的参数。

通过图6中的拟合可得K，n的值，见表1。

表1 不同循环次数下K与n的值Table 1 Values of K and n under different numbers of cycles

假设岩石在循环荷载过程中K，n的值与循环次数N有以下关系：

K=a+bNC

(3)

n=d+eNf

(4)

式中：a,b,c,d,e,f分别为与岩石力学性质相关的参数。

则式(2)可以改写为：

Ud=(a+bNc)σ(d+eNf)

(5)

则式(5)为一定循环次数下耗散能随应力的变化方程，其中当a=0.000 285 7,b= 0.000 043 92,c=0.640 7,d=2.039 7,e=-0.004 55,f=0.768 22时可以得到较好的结果，如图7所示。

图7 耗散能随应力变化的关系曲线Fig.7 Relation curve of dissipated energy with stress

3 结论

1)试件峰值强度前吸收的总能量、积聚弹性应变能和耗散能都随着轴向应力呈现出非线性增长，总能量速率增长最快，弹性应变能次之，耗散能最慢。

2)弹性应变能在峰值强度附近达到储能极限，增长速率逐渐降为0；峰后则是前期积累的弹性应变能急剧释放，耗散能随裂纹的发展而快速增长。

3)单次循环耗散能数值随着循环次数的增长而加大，说明周期荷载作用下岩石内部的损伤在高循环次数下有着扩大的趋势。

4)基于材料阻尼理论，建立了一定循环次数下耗散能随应力变化的演化方程，其计算结果与试验数据吻合性较好。