图1 PSC模型
Fig.1 PSC model
胡光辉，徐 涛，陈崇枫，杨学凯
(东北大学 岩石破裂与失稳研究中心, 沈阳 110819)
摘 要： 为从细观角度探究脆性岩石的蠕变失稳过程及失稳机理，该文基于三维颗粒流程序(PFC 3D )考虑岩石的时效变形损伤过程，引入岩石细观单元时效损伤的应力腐蚀模型，建立了基于离散元方法的岩石时效变形损伤破裂模型，并通过单轴压缩及单轴蠕变的室内实验和数值模拟对比验证了所建立的时效变形损伤破裂模型的合理性。数值模拟再现了岩石的初始蠕变、稳态蠕变和加速蠕变三个蠕变阶段，同时模拟结果表明，在单级加载条件下，随着应力水平提高，稳态蠕变应变率显著增大，岩石蠕变失效时间逐渐缩短，初始轴向应变、初始侧向应变和初始体应变不断增大，且细观裂纹扩展形式与单轴压缩破坏形式基本相同，都是以拉伸裂纹为主，裂纹的增长速率随着时间增加而不断增大，尤其在第三蠕变阶段裂纹增长速率迅速增大；在分级加载试验过程中，模型的轴向应变、侧向应变和体应变以及裂纹最终扩展形态与单级加载基本相同；此外将三维蠕变模拟结果与二维模拟结果进行对比，结果显示三维模型拟合程度更高。
关键词： 岩石蠕变；三维颗粒流程序；应力腐蚀模型；蠕变应变率；裂纹扩展
岩石蠕变是指岩石在低于其瞬时抗压强度的恒定外载荷作用下，变形不断增加而导致岩石发生失稳破坏的现象 [1] 。岩石蠕变性质是岩石的重要力学特性之一，很多岩石工程，如地下核废料储存洞室、岩石边坡的长期稳定性都与岩石蠕变有着密切关系。因此，开展岩石蠕变的相关研究有着重要的理论指导和工程实践意义 [2―4] 。目前人们对于岩石蠕变的研究多采用理论模型分析和室内试验研究。如杨圣奇和徐鹏 [5] 基于有效应力概念建立了岩石非线性损伤流变模型，该模型能够较好地描述岩石的衰减、稳态和加速流变阶段；Heap等 [6] 开展了相同孔压、不同围压作用下的岩石蠕变试验，基于微结构分析和声发射技术研究了岩石蠕变过程中声发射及孔隙体积变化规律；蒲成志等 [7] 以流变学理论为基础，将损伤因子引入弹性元件中，给出了考虑损伤效应的五元件变参数非线性蠕变损伤模型，并通过试验验证了该模型能够很好地描述蠕变的时效非线性特征。
近年来，数值分析方法已成为岩石蠕变特性分析的重要手段 [8―10] ，如 Lockner和 Madden [8] 基于重整化群理论运用有限差分方法开展了脆性岩石的蠕变变形模拟；Xu等 [9] 考虑材料的非均匀性基于材料力学性质的劣化特性，运用有限元方法开展了脆性岩石的细观蠕变损伤失稳机理研究；Chen和Konietzky [10] 基于亚临界裂纹扩展理论运用块体离散元软件开展了脆性岩石的蠕变失稳数值模拟研究，Tran等 [11] 在研究大坝的稳定性时提出了rock-ageing模型模拟坝体的蠕变问题，随后 Wang等 [12] 基于 rock-ageing理论提出了能够在 Particle Flow Code (PFC)中模拟蠕变的bond-ageing模型；Potyondy [13] 在应力腐蚀理论基础上基于二维颗粒流程序PFC 2D 建立了平行粘结应力腐蚀模型，在此基础上一些学者 [14―15] 也应用二维颗粒流程序开展了岩石蠕变失稳的离散元模拟研究。上述研究基本都是关于岩石蠕变的二维模拟，然而岩石都是三维结构，二维模型无论是在微观结构还是近似程度上都与三维模型存在差别，若能开展岩石蠕变损伤的三维数值模拟研究，对于进一步认识岩石蠕变的细观失稳过程及其失稳机理具有重要的意义。杨振伟等 [16] 基于三维颗粒流程序采用伯格斯(Burgers)模型模拟了岩石蠕变过程，分析了伯格斯模型中的弹性系数、粘性系数以及摩擦因数对瞬时强度特性和流变特性的影响；Liu等 [17] 针对盐岩的蠕变特点提出了一种非线性的粘弹塑性蠕变模型，并将此蠕变模型在FLAC 3D 软件中实现，通过模拟结果与室内蠕变试验结果的对比发现该模型能够有效模拟盐岩在低应力水平下长期蠕变特征以及在高应力水平下出现加速蠕变特征，但是该蠕变模型需要标定的参数较多；其他学者也建立了不同的三维蠕变模型，如 Li等 [18] 建立了非定常的西原正夫模型研究胶结类矸石充填材料的三维蠕变特性，赵延林等 [19] 通过研究岩石在蠕变各阶段的非线性蠕变特征提出了新的非线性弹粘塑性流变模型，并提供了非线性元件参数的确定方法。虽然以上方法都较好地模拟了岩石在各蠕变阶段的变形特征，尤其是加速蠕变阶段的变形特征，但这些模型包含较多需要确定的参数，求解这些模型参数增加了研究的工作量且无法监测模型内部的损伤状况。基于此，本文在三维颗粒流程序PFC 3D 中引入应力腐蚀模型，建立基于离散元方法的岩石时效变形损伤破裂模型，并利用此模型对岩石在不同应力水平下的蠕变失稳破坏过程进行数值模拟研究，从细观角度揭示其蠕变失稳机理，同时构建了二维蠕变模型，在相同加载条件下对二维模拟结果与三维模拟结果进行对比，以期对岩石工程长期稳定性等问题提供一些参考和理论指导。
颗粒流的研究始于20世纪50年代，1971年Cundall率先提出离散单元法 [20] ，随后 Cundall和Strack [21] 开发出适用于岩土力学的颗粒流程序PFC，其本质是利用显式差分算法和离散元理论开发的细观力学程序。颗粒流模型将离散介质视为颗粒集合体，颗粒与颗粒之间采用接触模型粘结起来，并从细观力学角度研究介质系统宏观力学特征和力学响应。在进行颗粒流数值模拟时，通过颗粒总的运动轨迹反映材料的宏观变化，而颗粒的运动轨迹受制于颗粒的细观参数，所以有效地确定颗粒间的细观参数是进行颗粒流数值模拟的关键，由于颗粒间细观参数的不确定性和复杂性，目前还没有统一的公式能够建立颗粒间的细观参数和材料宏观参数间对应定量关系，故需要不断地进行参数调试，这个过程称为参数标定。基于这些原因，本文主要是通过两步骤参数标定建立能反映岩石蠕变特性的细观模型：第一步，标定数值模型瞬态细观参数，使数值模型的弹性模量和峰值应力等力学参数与物理模型相匹配；第二步，标定岩石蠕变细观参数，使岩石数值模拟试验曲线与室内蠕变试验曲线相匹配。
岩石在恒定载荷作用下其力学参数随时间不断劣化，犹如腐蚀液不断对岩体的内部结构和力学性质进行腐蚀，故称为应力腐蚀。Potyondy [13] 针对岩石蠕变过程中发生的应力腐蚀，提出了采用颗粒流离散元方法模拟岩石静态疲劳破坏过程的一种应力腐蚀模型(parallel-bonded stress corrosion model，PSC)。PSC模型是将岩石离散为球形颗粒，颗粒之间的接触采用平行粘结模型如图1(a)表示。在PSC模型中，首先假设应力腐蚀只作用于颗粒之间的平行粘结接触，根据应力和裂纹扩展速度间的关系使颗粒间的胶结物逐渐收缩直至完全消失，从而模拟裂纹的形成或扩展。在应力腐蚀过程中，裂纹的产生或者扩展的驱动力是两颗粒间胶结处表面的局部应力，而且只有当颗粒间的接触为受拉状态且拉应力超过起裂应力阈值时，应力腐蚀才开始发生。
本文为了研究岩石在蠕变过程中的微裂纹扩展情况，通过引入离散裂隙网络 discrete fracture network (DFN) [22―23] ，将每一次的接触断裂视为一次微裂纹的产生。通过统计产生 DFN的数量，可以对岩石蠕变过程中的微裂纹进行监测；同时记录每次产生 DFN的位置、方向等信息，可以在程序中实时显示微裂纹的位置和贯通情况。
图1 PSC模型
Fig.1 PSC model
参考图2(b)，PSC模型的微裂纹扩展速度模型可以表示为：
式中： 
为颗粒间平行粘结接触表面处最大拉应力； 
为应力阈值； 
为颗粒间平行粘结抗拉强度； 
为平行粘结接触直径；β 1 和β 2 为材料常数。
根据PFC 3D 基本理论，颗粒间接触处的应力计算准则为：
式中： 
、 
分别为接触处的法向应力和切向应力； 
、 
分别为接触处的法向力和切向力； 
、 
分别为接触处惯性矩和极惯性矩；| 
|、| 
|分别为接触处的法向力矩和切向力矩； 
为平行粘结接触半径；A为平行粘结接触面积(PFC 3D 中接触为圆形，A=π 
)； 
为转动惯量系数(在PSC模型中设为1.0)。
当接触处的 
＞ 
或者 
＞ 
时，程序判定该平行粘结接触破坏，同时在数值模型中移除相应的力、力矩和刚度等力学参数，并将该处记录为一次微裂纹产生。
模型瞬态细观参数标定通常是根据室内单轴压缩试验中岩石试件所表现的宏观力学参数，如峰值应力、弹性模量以及泊松比作为标定基准，通过建立PFC数值模型，调节颗粒细观参数和接触细观参数使数值模型表现的宏观力学性质与室内岩石试件的力学参数相匹配。
室内试验采用φ50 mm×100 mm 岩石试样进行单轴压缩试验，其各宏观力学参数为：峰值应力47.27 MPa，弹性模量12.4 GPa，泊松比0.29。本文根据室内单轴压缩试验，对数值模型的细观力学参数进行标定。颗粒与颗粒之间的接触采用平行粘结接触(parallel-bonded contact)，建立φ50 mm×100 mm数值颗粒模型，颗粒细观参数和平行粘结接触细观参数见表1。
表1 数值模型细观参数
Table1 Micro-parameters of numerical model
采用表1中的细观参数，建立PFC静态数值模型，并进行单轴压缩试验，其应力-应变曲线如图2。由于岩石试件内部存在各种微孔洞和微裂纹，故在进行压缩试验时，随加载应力的增大，这些微孔洞和微裂纹首先被压密，应力应变曲线下凹，之后岩石进入弹性阶段。由于岩石弹性模量是以岩石在弹性阶段的斜率计算为准，故将试验曲线向左平移，以消去初始裂隙压密阶段 [24] ，可以看出模拟曲线与试验曲线总体上较好地相吻合。
图2 单轴压缩试验
Fig.2 Uniaxial compression test
由式(1)可知，能够控制 PSC模型的裂纹扩展速度的关键参数是 
、 
、β 1 和β 2 ，但 
属于颗粒间粘结材料的强度参数，是控制模型的静态强度，而 
、β 1 和β 2 是控制模型的长期强度。所以为分析PSC中各参数对模型长期强度的影响，分别对 
、β 1 、β 2 进行参数敏感性分析。
在进行参数敏感性分析之前，首先根据室内单轴蠕变试验结果，进行初期的试算从而确定各参数的大致范围 
：1 MPa～25 MPa，β：1×10 -15 m/s～ 1 25×10 -15 m/s，β 2 ：15～30。
根据上述各参数的大致范围，设计了 18组模拟试验，设每组载荷为70%应力水平，在试验过程中同时记录蠕变失效时间 t f (以加速蠕变阶段中轴向应变发生突变的时间为准)，并令T f =log t f 。具体参数设置如表2。
表2 PSC模型参数敏感性分析
Table2 Parameter sensitivity analysis of PSC model
1)根据试验1～试验6可知：β 2 与T f 近似呈指数关系，随着β 2 的增大，T f 急剧减小，且当β 2 取0时，蠕变失效时间为一常数，其分析结果如图3所示。
图3 β 2 与T f 的关系曲线
Fig.3 Relationship between β 2 and T f
2)根据试验7～试验12可知：β 1 与T f 大致呈幂函数关系，随着β 1 的增大，T f 有急剧减小的趋势，且当β 1 取0时，模型将长期处于稳定蠕变阶段，不会发生蠕变断裂；当β 1 大于某个数值时，试样瞬时断裂，其分析结果如图4所示。
图4 β 1 与T f 的关系曲线
Fig.4 Relationship between β 1 and T f
3)根据试验13～试验18： 
与T f 呈线性关系，随着 
的增大，T f 逐渐增加。原因在于增大 也就提高了颗粒之间激发PSC模型的门槛值，在相同加载条件下减少了发生损伤的接触的数量，模型更难被破坏。
图5 
与T f 的关系曲线
Fig.5 Relationship between 
and T f
通过对 PSC模型参数的敏感性分析可知，β 1 和β 2 对蠕变失稳时间的影响较大， 
与模型的失稳时间T f 呈正相关。所以在确定PSC模型的各参数时，可以先将 
设一个较大的值，并将β 1 和β 2 从较大值逐渐向较小值过渡，将模型蠕变失稳时间与室内蠕变失稳时间调整到同一数量级，之后在确定 
。
根据已经标定的瞬态细观参数，通过调试PSC模型的主要参数β 1 、β 2 和 
，将单轴蠕变数值模拟结果与室内单轴蠕变试验结果(70%应力水平)进行对比，确定一组合理的 PSC参数： 
=10 MPa，β 1 =8×10 -15 m/s，β 2 =18.5。
图6 单轴蠕变试验
Fig.6 Uniaxial creep test
图6为数值模拟结果与室内蠕变试验结果的对比。由图6可以看出数值模拟结果同物理实验结果较好地相吻合，模拟结果不仅能够有效模拟出室内蠕变试验的初始蠕变、稳态蠕变和加速蠕变三个阶段及其相关特征(初始瞬时弹性应变ε e 、稳态蠕变速率 2 ε˙和蠕变失效时间t f )，而且还很好地展示了微裂纹萌生、发展及其贯通形成宏观剪切带的过程。
基于上述蠕变模型开展不同应力水平(50%、60%、70%、80%、90%的单轴抗压强度)下的岩石单轴压缩蠕变试验，试验结果如图7所示。
图7 不同应力水平下岩石轴向蠕变曲线
Fig.7 Axial strain creep curves of rock samples at different stress levels
由图7可见，在50%和60%应力水平下，岩石轴向应变基本不发生变化，在此应力水平下岩样蠕变曲线经历初始蠕变和稳态蠕变两个阶段并保持长期稳定。在70%、80%和90%应力水平下，岩石出现典型的初始蠕变、稳态蠕变和加速蠕变三个阶段，而且稳态蠕变应变率随着应力水平的提高而增大，如图8所示，其稳态蠕变阶段的蠕变应变率分别为 5.79×10 -9 s -1 、3.32×10 -8 s -1 和 1.10×10 -7 s -1 ，可见应力水平的一个微小的变化，都会导致岩石在蠕变过程中的应变速率的巨大变化。
图9为不同应力水平下的侧向应变随时间变化关系曲线。由图9可见，岩石在单轴压缩蠕变过程中侧向应变与轴向应变有相同的变化趋势。
图 10为不同应力水平与对应的稳态蠕变应变率( 2 ε˙)关系曲线，由图可见，无论是轴向稳态蠕变率还是侧向稳态蠕变率均与应力水平近似成线性关系，数值模拟结果同已有实验结果相一致 [6, 25] 。从细观角度分析，随着应力水平的不断提高，模型内部接触达到应力腐蚀阈值而断裂的数量增多，在单位时间内模型的损伤程度显著增大，在宏观上表现为应变增加，应变速率增大。在单级加载稳定蠕变阶段，模型侧向蠕变率大于轴向蠕变率，表明在此阶段，裂纹萌生导致的侧向变形大于压缩引起的轴向变形，这也是岩石蠕变实验过程中岩样经历体积缩小、扩容并膨胀的主要原因。
图8 不同应力水平下岩石轴向应变率随时间变化曲线
Fig.8 Axial creep rate vs time curves of rock samples at different stress levels
图9 不同应力水平下岩石侧向应变随时间变化曲线
Fig.9 Lateral strain creep curves of rock samples at different stress levels
图10 稳态蠕变率( 2 ε˙)与应力关系
Fig.10 Relationship curves between steady creep rate ( 2 ε˙)and stress levels
不同应力水平下岩石蠕变过程中其体积变化曲线如图11所示。由图11可见，在较低的应力水平下(如图中50%和60%应力水平)，岩样的体积经历初始蠕变阶段的压缩和稳态蠕变阶段的膨胀，但由于轴向压缩变形程度大于侧向膨胀变形程度，模型试样的体积整体表现为减小，而在高应力水平下(如图7和图9中70%、80%和90%应力水平)，在初始蠕变阶段，由于轴向压缩变形程度大于侧向膨胀变形程度，模型的体积整体表现为减小；在稳定蠕变阶段，模型侧向蠕变率大于轴向蠕变率(如图10所示)，模型发生侧向膨胀变形，体积收缩趋势逐渐减小；在加速蠕变阶段，因颗粒接触单元断裂的数量级的增长，模型岩样的侧向膨胀变形程度显著大于轴向压缩变形，模型的体积迅速从压缩状态过渡到膨胀状态，岩样体积急剧增加，宏观表现为体积应变曲线的急剧增大及岩样的失稳破坏，如图 11所示。该现象已被文献[6]中的试验所证实。
图11 不同应力水平下岩石体应变随时间变化曲线
Fig.11 Volume strain vs time curves of rock samples at different stress levels
由于在50%和60%应力水平下，岩石内部颗粒间很少有接触应力达到或超过应力阈值，几乎没有微裂纹产生，为了深入研究岩石在发生蠕变失稳破坏时的微裂纹扩展情况，故在此增加75%和85%应力水平的单轴蠕变试验。
由于不同应力水平下，岩石蠕变破坏时间差异较大，为了研究蠕变过程中微裂纹数量的变化趋势，将蠕变时间按照t/t f 进行归一化。t f 为岩石蠕变破坏时间，t为蠕变损伤时间。
图12 不同应力水平下微裂纹数量随t/t f 变化趋势
Fig.12 Variation of number of microcracks with relative time t/t f at different stress levels
图12为不同应力水平下微裂纹数量随t/t f 变化曲线，从图12可以看出，微裂纹的数量变化大致有以下三个阶段。1)初始阶段：在此阶段微裂纹的数量基本为0，而且随着应力水平的提高，模型处于此阶段的时间逐渐减少。从细观角度分析，随应力水平的提高，模型内部达到或超过应力阈值的接触数量逐渐增多，模型的损伤速率逐渐增大。2)稳定增长阶段：在此阶段，微裂纹数量近似成线性增加，而且随应力水平的提高，该阶段的斜率逐渐增大。从细观角度分析，随应力水平的提高，模型内部的损伤逐渐增加，破坏加剧，微裂纹增长速率增大。3)急剧增加破坏阶段：在此阶段，微裂纹数量呈指数型增长。从细观角度分析，由于在该阶段之前，模型内部的已经产生了一定数量的微裂纹，这些微裂纹的产生改变了模型内部的微观结构，使得模型的承载能力降低，微裂纹数量迅速增加，模型发生失稳破坏。在各应力水平下，拉伸微裂纹数量和剪切微裂纹数量有同样的变化趋势，如图13所示。
图13 不同应力水平下微裂数量变化曲线
Fig.13 Curves of microcracks number at different stress levels
为了进一步揭示蠕变过程中损伤的时空演化过程，记录了不同应力水平作用下蠕变过程中微裂纹数量和微裂纹的空间位置信息，如图 14所示。图 14中深色代表拉伸破坏，浅色代表剪切破坏。从微裂纹分布图 14可以看出，单轴蠕变情况下模型的破坏是以拉伸破坏为主，剪切破坏为辅，最终形成宏观剪切破坏带。这与单轴压缩条件下的岩石破坏模式十分相似。随应力水平的提高，模型内部产生微裂纹的空间位置更加离散，除在宏观破坏带处产生大量的微裂纹外，在宏观破坏带之外，也产生了较多数量的微裂纹。这些微裂纹的产生，加剧了模型的失稳破坏程度，缩短了蠕变失稳时间。
图14 不同应力水平下微裂纹分布图
Fig.14 Microcracks at different stress levels
在室内试验中，由于各岩石试样的力学性质很难完全相同，所以有时试验数据会有很大的离散性，故通常采用分级加载试验，以降低试验数据的离散性和缩短试验时间。在PFC 3D 中，通过控制墙体的位移，实现对岩石模型的分级加载。从60%应力水平到90%应力水平，共分为4个加载等级(60%、70%、80%和 90%)，每一级加载模拟 4 h，在载荷达到90%应力水平时，继续保持应力恒定，直到模型破坏。
模型在分级加载条件下，轴向应变与横向应变都具有蠕变特征。在每一级应力加载初始阶段模型产生弹性变形，且弹性变形较大，之后轴向应变与横向应变进入稳定发展阶段。在较低应力水平下，模型进入蠕变阶段后，轴向应变增长缓慢，最后趋于稳定，随应力水平的提高，模型应变值迅速增大，最终模型在高应力条件下发生失稳破坏。从细观角度分析，在较低应力水平下模型内部接触达到应力腐蚀阈值而断裂的数量稳定，轴向应变与侧向应变基本保持不变(如图15中的70%、80%和90%应力水平)，在蠕变破坏后期，模型内部微观承载结构的力学性质进一步劣化，模型内大量的接触无法继续承载，从而导致模型在宏观上应变值迅速增加，最终发生失稳破坏。
图15 分级加载试验
Fig.15 Stress-stepping test
为了能够更好地比较三维和二维离散元 PSC模型在模拟岩石蠕变试验结果上的准确性，本文还进行了在相同加载条件下的二维单轴蠕变模拟试验。同建立三维离散元蠕变模型过程一样，在建立二维岩石蠕变模型之前，首先根据单轴压缩试验标定出二维数值模型所需要的各细观参数，最终标定的二维单轴压缩试验曲线如图16所示。
图16 二维与三维模型单轴压缩模拟计算结果
Fig.16 Uniaxial compression strength test results of 2D and 3D models
由图 16可知，二维离散元模型同样能够有效模拟出岩石在单轴压缩过程中所表现的各宏观力学性质，并能够很好地拟合试验曲线。在此二维数值模型基础上，参照建立三维蠕变模型的过程，建立能够反映岩石蠕变特征的二维模型，并根据PSC参数敏感性分析所得的结果，确定一组二维PSC模型参数，标定结果如图17所示。
图17 二维与三维模型单轴蠕变模拟计算结果
Fig.17 Uniaxial compression creep test results of 2D and 3D models
由图 17可见，采用二维模型同样能大致模拟出岩石蠕变各阶段的特征，但与三维模型的模拟结果不同之处在于初始减速蠕变阶段和加速蠕变阶段。在初始阶段，三维模型模拟结果能够有效地模拟出岩石的瞬时弹性应变，而后进入到减速蠕变阶段，二维模型同样模拟出了岩石在蠕变初期阶段蠕变速率递减的规律，但其轴向应变值低于室内试验值和三维模拟结果，与室内试验结果差距较大；在加速蠕变阶段，三维模拟结果很好地模拟出了岩石在加速蠕变阶段中轴向应变迅速增大直到破坏的过程，而二维模型虽然模拟出了岩石出现加速蠕变的大致时间，但其轴向应变值增加缓慢，蠕变速率远小于室内试验和三维模拟结果。分析其原因，可能是由于三维颗粒集合体存在堆积效应，在相同颗粒数量条件下，三维颗粒之间的接触数量多于二维模型的概率大为增加，使得三维模型在发生破坏前相较与二维模型更加稳定，二维颗粒结构受接触的形成和断裂的影响更大，故在初期减速蠕变阶段和后期加速蠕变阶段二维模型的轴向应变是随着接触的形成和破坏而逐步增大，相较于三维模拟结果并非是突然变化的，即在初期轴向应变小于三维模拟，而后期轴向蠕变速率远小于三维蠕变结果。
(1)基于三维颗粒流程序 PFC 3D 引入应力腐蚀PSC模型，建立了基于离散元方法的岩石时效变形损伤破裂过程模型，并应用此模型开展了岩石单轴压缩和单轴压缩蠕变试验，并与室内试验进行了对比验证，发现该模型能够反映岩石的蠕变特性并有效模拟岩石的初始蠕变、稳态蠕变和加速蠕变三个阶段。
(2)单级加载条件下，岩石是否发生加速蠕变与所受的载荷大小、力学性质等有密切关系。当岩石所受载荷低于岩石的长期强度时，岩石内部未发生损伤或损伤程度很小，岩石的轴向应变未出现加速阶段，但在所受载荷高于长期强度时岩石的轴向应变出现明显的加速阶段，模型失稳破坏。侧向蠕变应变和体积蠕变应变曲线有着相似的变化趋势。
(3)引入离散裂隙网络DFN模型监测岩石蠕变过程中微裂纹的时空演化规律，发现在岩石单级蠕变加载试验过程中，微裂纹的扩展同样可分为三阶段：初始阶段、稳定增长阶段和急剧增加破坏阶段。岩石蠕变损伤失稳过程中裂纹以拉伸裂纹为主，剪切裂纹为辅，最终形成宏观剪切带。
(4)分级加载条件下，在低应力水平时岩石内部颗粒孔隙逐渐被压实，体应变收缩的趋势最终趋于稳定。随应力水平的提高，微裂纹在蠕变稳态阶段后期(即加速蠕变阶段)急剧增长，导致岩样体积膨胀，发生失稳破坏。
在相同加载条件下，二维模型同样能大致模拟岩石蠕变过程，但与室内试验和三维模拟结果相比，误差主要出现在减速蠕变阶段和加速蠕变阶段。
参考文献：
[1] 许宏发. 软岩强度和弹模的时间效应研究[J]. 岩石力学与工程学报, 1997, 16(3): 246―251.Xu Hongfa.. Time dependent behaviours of strength and elasticity modulus of weak rock [J]. Chinese Journal of Rock Mechanics and Engineering, 1997, 16(3): 246―251. (in Chinese)
[2] 李良权, 徐卫亚, 王伟, 等. 基于流变试验的向家坝砂岩长期强度评价[J]. 工程力学, 2010, 14(11): 127―136.Li Liangquan, Xu Weiya, Wang Wei, et al. Estimation of long-term strength for xiangjiaba sandstone based on creep tests [J]. Engineering Mechanics, 2010, 14(11):127―136. (in Chinese)
[3] 王贵君, 张磊, 张昱文, 等. 盐岩流变-损伤-破坏特性的实验研究[J]. 工程力学, 2013, 30(4): 288―293.Wang Guijun, Zhang Lei, Zhang Yiwen, et al.Experimental investigation on the creep-damage-rupture characteristics of salt rock [J]. Engineering Mechanics,2013, 30(4): 288―293. (in Chinese)
[4] 周广磊, 徐涛, 朱万成, 等. 基于温度-应力耦合作用的岩石时效蠕变模型[J]. 工程力学, 2017，35(10): 1―9.Zhou Guanglei, Xu Tao, Zhu Wancheng, et al. A time-dependent thermo-mechanical creep model of rock[J]. Engineering Mechanics, 2017, 35(10): 1―9. (in Chinese)
[5] 杨圣奇, 徐鹏. 一种新的岩石非线性流变损伤模型研究[J]. 岩土工程学报, 2014, 36(10): 1846―1854.Yang Shengqi, Xu Peng. A new nonlinear rheological damage model for rock [J]. Chinese Journal of Geotechnical Engineering, 2014, 36(10): 1846―1854.(in Chinese)
[6] Heap M J, Baud P, Meredith P G, et al. Time-dependent brittle creep in Darley Dale sandstone [J]. Journal of Geophysical Research: Solid Earth, 2009, 114(B7): 1―22.
[7] 蒲成志, 曹平, 张春阳, 等. 考虑时效损伤劣化的变参数非线性蠕变损伤模型[J]. 工程力学, 2017, 34(6):17―27.Pu Chengzhi, Cao Ping, Zhang Chunyang, et al. Variable parameters nonlinear creep damage model of rock with consideration of aging, damage and deterioration [J].Engineering Mechanics, 2017, 34(6): 17―27. (in Chinese)
[8] Lockner D A, Madden T R. A multiple-crack model of brittle-fracture [J]. Journal of Geophysical Research-Solid Earth, 1991, 96(B12): 19643―19654.
[9] Xu T, Tang C A, Zhao J, et al. Modelling the time-dependent rheological behaviour of heterogeneous brittle rocks [J]. Geophysical Journal International, 2012,189(3): 1781―1796.
[10] Chen W, Konietzky H. Simulation of heterogeneity,creep, damage and lifetime for loaded brittle rocks [J].Tectonophysics, 2014, 633: 164―175.
[11] Tran T H, Vénier R, Cambou B. Discrete modelling of rock-ageing in rockfill dams [J]. Computers and Geotechnics, 2009, 36(1): 264―275.
[12] Wang Y, Zhao Z, Song E. Discrete element modeling of the effect of particle shape on creep behavior of rockfills[J]. International Journal of Environmental, Engineering,2017, 11(9): 803―807.
[13] Potyondy D O. Simulating stress corrosion with a bonded-particle model for rock [J]. International Journal of Rock Mechanics and Mining Sciences, 2007, 44(5):677―691.
[14] 刘宁, 张春生, 褚卫江. 深埋大理岩破裂扩展时间效应的颗粒流模拟[J]. 岩石力学与工程学报, 2011,30(10): 1989―1996.Liu Ning, Zhang Chunsheng, Chu Weijiang. Simulating time-dependent failure of deep marble with particle flow code [J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(10): 1989―1996. (in Chinese)
[15] 孙金山, 陈明, 姜清辉, 等. 锦屏大理岩蠕变损伤演化细观力学特征的数值模拟研究[J]. 岩土力学, 2013 34(12): 3601―3608.Sun Jinshan, Chen Ming, Jiang Qinghui, et al. Numerical simulation of mesomechanical characteristics of creep demage evolution for Jingping marble [J]. Rock and Soil Mechanics, 2013, 34(12): 3601―3608. (in Chinese)
[16] 杨振伟, 金爱兵, 周喻, 等. 伯格斯模型参数调试与岩石蠕变特性颗粒流分析[J]. 岩土力学, 2015, 36(1):240―248.Yang Zhenwei, Jin Aibing, Zhou Yu, et al. Parametric analysis of Burgers model and creep properties of rock with particle flow code [J]. Rock and Soil Mechanics,2015, 36(1): 240―248. (in Chinese)
[17] Liu X, Yang X, Wang J. A nonlinear creep model of rock salt and its numerical implement in FLAC 3D [J].Advances in Materials Science & Engineering, 2015,2015(5): 1―8.
[18] Li C, Ma Q, Chen G, et al. Study on Creep characteristics of cemented waste rock backfills [J]. Open Civil Engineering Journal, 2015, 9(1): 957―961.
[19] 赵延林, 曹平, 文有道, 等. 岩石弹黏塑性流变试验和非线性流变模型研究[J]. 岩石力学与工程学报, 2008,27(3): 477―486.Zhao Yanlin, Cao Ping, Wen Youdao, et al. Elastoviscoplastic rheological experiment and Nonlinear rheological model of rocks [J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(3): 477―486. (in Chinese)
[20] Cundall P A. A computer model for simulating progressive, large-scale movements in blocky rock system [C]// Proceedings of the Symposium of the International Society for Rock Mechanics, Society for Rock Mechanics (ISRM), France, Nancy. 1971:10018723276.
[21] Cundall P A, Strack O D L. A discrete numerical model for granular assemblies [J]. Geotechnique, 1979, 29(1):47―65.
[22] 王培涛, 杨天鸿, 于庆磊, 等. 基于离散裂隙网络模型的节理岩体渗透张量及特性分析[J]. 岩土力学, 2013,24(2): 448―455.Wang Peitao, Yang Tianhong, Yu Qinglei, et al.Permeability tensor and seepage properties for jointed rock masses based on discrete fracture network model[J]. Rock and Soil Mechanics, 2013, 24(2): 448―455.(in Chinese)
[23] Bisdom K, Nick H M, Bertotti G. An integrated workflow for stress and flow modelling using outcrop-derived discrete fracture networks [J].Computers & Geosciences, 2017, 103: 21―35.
[24] Hao S W, Wang H Y, Xia M F, et al. Relationship between strain localization and catastrophic rupture [J].Theoretical and Applied Fracture Mechanics, 2007,48(1): 41―49.
[25] Brantut N, Heap M J, Meredith P G, et al.Time-dependent cracking and brittle creep in crustal rocks: A review [J]. Journal of Structural Geology, 2013,52: 17―43.
HU Guang-hui , XU Tao , CHEN Chong-feng , YANG Xue-kai
(Centre for Rock Ⅰnstability & Seismicity Research, Northeastern University, Shenyang 110819, China)
Abstract: To analyze the creep process and creep mechanism of rock at a microscopic scale, a numerical model of time-dependent deformation and damage was established by incorporating the parallel-bonded stress corrosion model into the Particle Flow Code 3D (PFC 3D ).The numerical model was validated against laboratory results.The numerical results replicated the typical three creep stages including the primary creep, secondary creep and tertiary creep.Meanwhile, the numerical simulation results indicate that the time-to-failure of creep decreases gradually with the increase of the stress level, while the initial axial strain, initial lateral strain, initial volume strain and creep strain rate gradually increase in this process.The crack propagation pattern in the uniaxial creep test is generally identical to that in the uniaxial compression test and the tension cracks are dominant.The crack growth rate increases with time and sharply increases in the tertiary creep stage.Ⅰn the stress-stepping test, the axial strain, lateral strain, volumetric strain and the final propagation mode of the model are the same as that of the single-stage creep test.Ⅰn addition, the comparisons between two-dimensional and three-dimensional creep simulation show that the three-dimensional modeling is more suitable than the two-dimensional modeling.
Key words: rock creep; PFC 3D ; stress corrosion model; creep strain rate; crack propagation
A MICROSCOPIC STUDY OF CREEP AND FRACTURING OF BRITTLE ROCKS BASED ON DISCRETE ELEMENT METHOD
杨学凯(1993―)，男，河北人，硕士生，主要从事岩石力学研究(E-mail: 1396946033@qq.com).
陈崇枫(1991―)，男，山东人，博士生，主要从事岩石力学研究(E-mail: 871330060@qq.com);
作者简介：
胡光辉(1992―)，男，安徽人，硕士生，主要从事岩石力学研究(E-mail: huguanghuineu@qq.com);
基金项目： 国家重点基础研究发展计划(973计划)项目(2014CB047100)；国家自然科学基金项目(41672301，51474051)；中央高校基本科研业务费项目(N150102002)
通讯作者： 徐 涛(1975―)，男，湖北人，教授，博士，博导，主要从事岩石破裂失稳研究工作(E-mail: xutao@mail.neu.edu.cn).
文章编号： 1000-4750(2018)09-0026-11
收稿日期： 2017-05-15；修改日期：2017-12-08
文献标志码： A
doi: 10.6052/j.issn.1000-4750.2017.05.0356
中图分类号： TD315