单轴压缩下断续节理岩体动态损伤本构模型*

刘红岩,李俊峰,裴小龙

(1.中国地质大学(北京)工程技术学院,北京 100083;2.武警黄金部队第四支队,辽宁 辽阳 111000)

目前岩体工程如隧道、边坡及地下厂房等在爆炸、地震等动载下的力学响应及破坏机理已引起了极大的关注。而岩体的强度、破坏模式及变形等力学特性主要受控于节理等天然缺陷[1-3]。由于节理发育规模一般远小于断层，延伸长度大多为数米，但是其数目众多，因而无法一一考虑其对岩体力学性质的弱化效应。为此，20世纪80年代Kyoya等[4]最早提出了岩体损伤理论的概念，视节理为岩体的初始损伤，从而建立了节理岩体损伤模型。目前岩体损伤张量定义方法主要有2种。一种是Kawamoto等[5]和Swoboda等[6]提出的几何损伤张量定义方法，他们采用节理密度、表面积、法向矢量等参数来定义损伤张量。但是该方法中仅考虑了节理的长度和方位等几何参数对岩体力学特性的影响，而未考虑节理摩擦角及黏聚力等强度参数的影响，认为应力无法在节理位置传递，因而该方法更适用于岩体拉伸破坏的模拟，而无法用于岩体压缩力学特性的分析。这是因为在受压条件下，节理面将产生压缩闭合或剪切滑移，从而传递部分压、剪应力[7-8]。针对这一问题，Li等[9]和Liu等[10]基于损伤及断裂力学中的能量原理，提出了一个新的能够同时考虑节理的几何参数(如节理长度、倾角、组数等)及强度参数(如节理摩擦角)的岩体损伤张量计算方法。相对于前述的几何损伤理论而言，第2种方法减小了由于节理传压及传剪系数的人为选取所带来的误差，因而具有较强的合理性[10-11]。

然而，尽管相对于几何损伤理论而言，第2种方法已经有了较大的发展，但是仍不够完善。这主要是由于目前学术界为了更准确地刻画节理的物理力学性质，提出了描述节理物理力学性质的3类参数，即几何参数、强度参数和变形参数(如法向和切向刚度)。而Prudencio等[12]通过对节理岩体的压缩实验，发现岩体的力学特性同样受节理变形参数的影响。而节理变形参数在目前的节理岩体损伤张量计算中均未得到很好地考虑，这显然会导致一定的计算误差。因此，本文中拟基于Li等[9]和Liu等[10]的研究方法，在现有的、能够考虑节理几何及强度的损伤张量计算方法的基础上，考虑节理变形参数的影响，最终在经典岩石细观动态损伤本构模型—Taylor-Chen-Kuszmaul (TCK)模型的基础上，提出能够考虑宏观缺陷(如节理)与细观缺陷(如微裂纹)耦合的断续节理岩体单轴压缩动态损伤本构模型，并采用相关算例研究节理几何参数(如节理长度)、强度参数(如节理摩擦角)和变形参数(如节理法向及切向刚度)对岩体单轴压缩动态应力应变曲线的影响。

1 岩石动态损伤TCK模型

1986年，Taylor等[13]提出了著名的岩石动态损伤本构模型，即TCK模型，其相应的单轴动态压缩条件下轴向应力σ与应变ε的关系[14]为：

(1)

2 断续节理岩体损伤本构模型

2.1 单条断续节理引起的岩体损伤张量

Li等[9]推导了当岩体含有一条如图1所示的倾角为α的断续节理时，平面应力条件下，由其引起的加载方向上的岩体损伤变量：

(2)

式中：KⅠ、KⅡ分别为节理尖端的Ⅰ、Ⅱ型应力强度因子；V为试件体积；A为节理表面积；单个节理时，A=Ba(单边节理)或2Ba(中心节理)；多个节理时，A=NBa(单边节理)或2NBa(中心节理);N为节理个数，B为节理深度，a为节理半长。

2.2 单条断续节理的应力强度因子

节理在压缩荷载下往往会在已有节理尖端产生翼裂纹并近似沿外力方向扩展[15-17]，如图1所示。

对于完整岩石，在单轴压缩下，倾角为α的斜面上的正应力和切应力分别为：

(3)

(4)

而对于节理岩体，即当在倾角为α的斜面上存在一条长度为2a的断续节理时，该面上的正应力σα和切应力τα将明显受到节理力学性质的影响，可表示为：

σα=(1-Cn)σcos2α

(5)

τα=(1-Ct)σsinαcosα

(6)

由于节理面黏聚力对试件强度影响较小，因此可忽略不计，同时若设节理摩擦角为φ，节理面摩擦系数为tanφ，那么在单轴压缩下，由式(5)、(6)可得节理处的下滑力为：

(7)

根据Lee等[19]的研究，则可得：

(8)

式中：引入l*=0.27a，以保证l=0时，KI和KII不存在奇异；其余参数[15-17]如图1所示。

对于节理未扩展的初始状态，即翼裂纹长度l=0时，则有：

(9)

把式(9)代入式(2)即可得到由初始节理所引起的岩体损伤变量，进而通过张量化处理即可得到相应的损伤张量。可知由该方法求出的损伤张量同时考虑了节理的几何参数(如节理长度、倾角等)、强度参数(如节理摩擦角)和变形参数(如节理法向及切向刚度)，因此将更加符合实际情况。

2.3 单排及多排断续平行节理的应力强度因子

当岩体含有如图2所示的单排及多排断续平行节理时(其中，b、d分别为相邻两节理中心距及排间距)，节理尖端应力强度因子的计算方法已在相关文献中有较多阐述[20-21]，这里不再重述。

2.4 含单组节理的岩体损伤变量

当岩体含单组节理时，把式(5)～(6)、(9)及相应的节理尖端应力强度因子计算公式代入式(2)，进而通过求积分可得相应的损伤变量。然而由于被积函数比较复杂，式(2)没有相应的解析解，因此这里采用Matlab软件求其近似解。

2.5 损伤变量的张量化

由2.4节求得的损伤变量仅是加载方向上的损伤值，为反映节理对岩体损伤的奇异性，必须对其进行张量化处理。这里采用Chen等[22]提出的方法，即引入损伤张量：

(10)

式中：D0即为由式(2)求得的损伤值。

3 考虑宏细观缺陷耦合的断续节理岩体单轴压缩动态损伤本构模型

3.1 宏细观缺陷耦合的损伤变量计算方法

这里引用Liu等[10]提出的宏细观缺陷耦合损伤张量计算公式：

(11)

3.2 考虑宏细观缺陷耦合的断续节理岩体动态损伤本构模型

基于连续介质力学损伤理论，把式(1)所示的完整岩石损伤变量ω(t)，用式(11)所示的宏细观缺陷耦合损伤变量Ω12替换，即可得到相应的单轴压缩断续节理岩体动态损伤本构方程：

(12)

式中：E为完整岩石的弹性矩阵，其余参数同前。

4 算例分析

4.1 含单组断续节理的岩体动态力学特性

根据上述参数，计算结果如图4所示。可以看出：(1)采用本文中提出的模型所计算出的断续节理岩体动态应力应变曲线也包含了应力线性上升段、极值点及应力下降段等，很好地反映了岩体变形全过程；(2)断续节理岩体的峰值强度及弹性模量分别为25.65 MPa和8.32 GPa，而相应完整岩石的动态峰值强度及弹性模量则分别为32.7MPa和10.8GPa，前者分别约为后者的78.44%和77.04%，均低于完整岩石的峰值强度，即断续节理会显著削弱岩体的动态强度及刚度；(3)另外，采用同样的算例，刘红岩等[11]计算出的相应岩体峰值强度及弹性模量分别为10.28 MPa和3.25 GPa，分别约为相应完整岩石的动态峰值强度及弹性模的31.44%和30.09%，如图4所示，即远小于本文计算结果。刘红岩等[11]在计算由断续节理对岩体造成的损伤时考虑了节理倾角的影响，认为式(12)中的E应为断续节理岩体的弹性矩阵，进而视岩体为岩块与节理组成的复合体系来计算E，在此过程中，同时又引入了节理倾角对弹性模量的影响，因此可以认为刘红岩等[11]考虑了2次节理倾角对断续岩体动态应力应变曲线的影响，即重复计算了一次节理倾角对断续岩体动态应力应变曲线的影响，因此导致其计算出的岩体动态应力偏小，因而计算结果不甚合理。而相比而言，在本文所提出的模型中，对节理的3类参数，即几何参数、强度参数及变形参数等均只考虑了一次，不存在重复计算的问题。

4.2 节理摩擦角对岩体动力学特性的影响

计算模型如图3所示，设节理摩擦角分别为0°、15°和30°，其余计算参数同4.1节，计算结果如图5所示。由图5可知：(1)由岩体动态应力应变曲线斜率所反映的动态弹性模量是随着节理摩擦角的增大而增大的，即3种不同节理摩擦角条件下的岩体动态弹性模量分别为7.24、8.32和9.61 GPa；(2)当节理摩擦角逐渐增大时，相应试件的动态峰值强度也由22.34 MPa分别增大到25.65和29.42 MPa。这是因为试件宏观损伤程度随节理摩擦角的增大而降低，进而导致试件强度提高，即节理摩擦角对试件的动态力学特性有较大影响。

4.3 节理法向及切向刚度对岩体动力学特性的影响

计算模型仍为图3所示，首先取节理法向刚度Kn分别为2、20和200 GPa/cm，其余计算参数同4.1节；其次取节理切向刚度Ks分别为0.8、8和80 GPa/cm，其余计算参数仍同4.1节，计算结果如图6所示。可以看出：(1)当节理法向刚度由2 GPa/cm依次增大到20和200 GPa/cm时，即节理法向刚度依次增大到原来的10倍及100倍时，试件的峰值强度则由25.65 MPa依次降低到23.35和22.34 MPa，这说明随着节理法向刚度的增大，试件动态峰值强度是降低的。这是因为在同样的垂直应力下，随着节理法向刚度的增大，节理面上的法向应力会相应减小，进而节理面上抵抗节理发生滑移的摩擦力就会减小，最终导致节理更易发生滑移破坏，进而峰值强度降低。(2)当节理切向刚度由0.8 GPa/cm依次增大到8和80 GPa/cm时，即节理切向刚度依次增大到原来的10倍及100倍时，试件的峰值强度则由25.65 MPa依次增大到31.82和32.32 MPa，这说明节理切向刚度对试件的动态峰值强度有较大影响，且随着节理切向刚度逐渐增大，节理切向刚度对岩体动态峰值强度的影响则逐渐趋于定值。(3)如前所述，由图6(a)、(b)分别看出，随着节理法向及切向刚度的增大，试件的动态峰值强度分别减小及增大。然而Prudencio等[12]研究发现随着节理法向及切向刚度的增大，试件动态峰值强度则是增大的。分析认为对图6(a)、(b)两种工况，节理的法向及切向刚度的选取是任意的。而节理法向及切向刚度则首先是由Goodman等[23-24]提出用以描述法(切)向应力与法(切)向位移之间的关系。Bandis等[25]发现节理的切向刚度远小于其法向刚度，并研究了Kn/Ks的值随法向应力的变化关系。这说明二者应该是按比例变化的。然而图6(a)、(b)中所给出的结果却并没有考虑这一因素。为此，下面就研究二者按比例变化时岩体的动态应力应变曲线。取节理法向及切向刚度的初始值分别为2和0.8 GPa/cm，然后分别增大到10倍及100倍，结果如图6(c)所示。可以看出，岩体的动态峰值强度则由26.65 MPa分别增大到30.36和32.64 MPa。这说明随着节理法向及切向刚度的增大，岩体动态峰值强度是增大的，这与Prudencio等[12]的研究成果则是一致的。(4)从应力应变曲线的斜率还可以看出，节理法向及切向刚度对试件的动态弹性模量也均有一定的影响。总之，节理法向及切向刚度对试件动态峰值强度及弹性模量均有一定影响。

4.4 节理长度对岩体动力学特性的影响

计算模型如图3所示，取节理长度2a分别为10、20和30 mm，其余计算参数同4.1节,计算结果如图7所示。由图7可知,当2a由10 mm依次增大到20和30 mm时，试件峰值强度则由30.57 MPa依次降低到25.65和20.63 MPa，即分别减小为完整岩石峰值强度的93.49%、78.44%和63.09%，即随节理长度的增大，试件动态峰值强度明显下降。这是由于节理长度的增大引起岩体损伤变量的增大及节理传压和传剪系数的变化，进而导致强度降低。同时岩体动态弹性模量也随节理长度的增大而减小。

4.5 载荷应变率对岩体动力学特性的影响

计算模型如图3所示，取载荷应变率分别为100、200和300 s-1，其余计算参数同4.1节，计算结果如图8所示。由图8可知，当载荷应变率由100 s-1依次升高到200和300 s-1时，试件峰值强度则由25.65 MPa依次增大到31.65和35.71 MPa，即试件具有明显的应变率效应。关于应变率对岩体动态力学特性的影响机理已有较多阐述[11-26]。

5 结 论

(1)针对目前节理岩体宏观损伤张量计算中仅考虑了节理几何及强度参数的不足，根据前人的研究成果，提出了能够同时考虑节理几何参数(如节理长度及倾角等)、强度参数(如节理摩擦角)及变形参数(如节理法向及切向刚度)的宏观损伤张量计算方法，这将使得宏观岩体损伤张量的计算更符合实际。

(2)基于完整岩石动态TCK模型及宏细观损伤耦合的观点，建立了考虑宏细观缺陷耦合的非贯通节理岩体动态单轴压缩损伤本构模型。该模型中同时考虑了节理几何、强度及变形参数对非贯通节理岩体动态力学特性的影响。

(3)算例分析表明，节理的存在降低了岩体的动态弹性模量，削弱了岩体的峰值强度。同时利用该模型讨论了节理摩擦角、节理法向及切向刚度、节理长度及应变率等对岩体动态力学特性的影响，计算结果与已有研究成果吻合较好，说明该模型是合理性的。

[1] AKIN M. Slope stability problems and back analysis in heavily jointed rock mass: A case study from Manisa, Turkey[J]. Rock Mechanics and Rock Engineering, 2013,46(2):359-371.

[2] 张力民,吕淑然,刘红岩.综合考虑宏细观缺陷的岩体动态损伤本构模型[J].爆炸与冲击,2015,35(3):428-436.

ZHANG Limin, LÜ Shuran, LIU Hongyan. A dynamic damage constitutive model for rock mass by comprehensively considering macroscopic and mesoscopic flaws[J]. Explosion and Shock Waves, 2015,35(3):428-436.

[3] 刘红岩,杨艳,李俊峰,等.基于TCK模型的断续节理岩体动态损伤本构模型[J].爆炸与冲击,2016,36(3):319-325.

LIU Hongyan, YANG Yan, LI Junfeng, et al. Dynamic damage constitutive model for rock mass with non-persistent joints based on the TCK model[J]. Explosion and Shock Waves, 2016,36(3):319-325.

[4] KYOYA T, ICHIKAWA Y, KAWAMOTO T. A damage mechanics theory for discontinuous rock mass[C]∥Proceedings of the 5th International Conference on Numerical Methods in Geomechanics. Nagoya, 1985:469-480.

[5] KAWAMOTO T, ICHIKAWA Y, KYOYA T. Deformation and fracturing behavior of discontinuous rock mass and damage mechanics theory[J]. International Journal for Numerical Analysis Method in Geomechanics, 1988,12(1):1-30.

[6] SWOBODA G, SHEN X P, ROSAS L. Damage model for jointed rock mass and its application to tunneling[J]. Computers and Geotechnics, 1998,22(3/4):183-203.

[7] YUAN X P, LIU H Y, WANG Z Q. An interacting joint-mechanics based model for elastoplastic damage model of rock-like materials under compression[J]. International Journal of Rock Mechanics & Mining Sciences, 2013,58(9):92-102.

[8] SWOBODA G, YANG Q. An energy-based damage model of geomaterials:Ⅰ: Formulation and numerical results[J]. International Journal of Solids and Structures, 1999,36(9):1719-1734.

[9] LI N, CHEN W, ZHANG P, et al. The mechanical properties and a fatigue-damage model for jointed rock mass subjected to dynamic cyclical loading[J]. International Journal of Rock Mechanics & Mining Sciences, 2001,38(7):1071-1079.

[10] LIU Hongyan, ZHANG Limin. A damage constitutive model for rock mass with non-persistently closed joints under uniaxial compression[J]. Arabian Journal for Science and Engineering, 2015,40(1):3107-3117.

[11] 刘红岩,王新生,张力民,等.断续节理岩体单轴压缩动态损伤本构模型[J].岩土工程学报,2016,38(3):426-436.

LIU Hongyan, WANG Xinsheng, ZHANG Limin, et al. A dynamic damage constitutive model for rock mass with non-persistent joints under uniaxial compression[J]. Chinese Journal of Geotechnical Engineering, 2016,38(3):426-436.

[12] PRUDENCIO M, JAN M V S. Strength and failure modes of rock mass models with non-persistent joints[J]. International Journal of Rock mechanics & Mining Sciences, 2007,46(6):890-902.

[13] TAYLOR L M, CHEN E P, KUSZMAUL J S. Microjoint induced damage accumulation in brittle rock under dynamic loading[J]. Computer Method in Applied Mechanics & Engineering, 1986,55:301-320.

[14] GRADY D E, KIPP M E. Continuum modeling of explosive fracture in oil shale[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1987,17(3):147-157.

[15] HUANG C, SUBHASH G, VITTON S J. A dynamic damage growth model for uniaxial compressive response of rock aggregates[J]. Mechanics of Materials, 2002,34(5):267-277.

[16] HUANG C, SUBHASH G. Influence of lateral confinement on dynamic damage evolution during uniaxial compressive response of brittle solids[J]. Journal of the Mechanics and Physics of Solids, 2003,51(6):1089-1105.

[17] PALIWAL B, RAMESH K T. An interacting micro-joint damage model for failure of brittle materials under compression[J]. Journal of the Mechanics and Physics of Solids, 2008,56(3):896-923.

[18] LIU Taoying, CAO Ping, LIN Hang. Damage and fracture evolution of hydraulic fracturing in compression-shear rock cracks[J]. Theoretical and Applied Fracture Mechanics, 2014,74:55-63.

[19] LEE S, RAVICHANDRAN G. Joint initiation in brittle solids under multiaxial compression[J]. Engineering Fracture Mechanics, 2003,70(13):1645-1658.

[20] 李建林,哈秋瓴.节理岩体拉剪断裂与强度研究[J].岩石力学与工程学报,1998,17(3):259-266.

LI Jianlin, HA Qiuling. A study of tensile-shear joint and strength related to jointed rock mass[J]. Chinese Journal of Rock Mechanics and Engineering, 1998,17(3):259-266.

[21] 范景伟,何江达.含定向闭合断续节理岩体的强度特性[J].岩石力学与工程学报,1992,11(2):190-199.

FAN Jingwei, HE Jiangda. The strength behavior of rockmasses containing oriented and closed intermittent joints[J]. Chinese Journal of Rock Mechanics and Engineering, 1992,11(2):190-199.

[22] CHEN W, LA BORDERIE C, MAUREL O, et al. Simulation of damage-permeability coupling for mortar under dynamic loads[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2014,38(5):457-474.

[23] GOODMAN R E. The mechanical properties of joints[C]∥Proceeding of the 3rd Congress ISRM. Denver, 1974,I(A):127-140.

[24] GOODMAN R E, TAYLOR R L, BREKKE T. A model for the mechanics of jointed rock[J]. Journal of Soil Mechanics and Foundations Division, 1968,94:637-659.

[25] BANDIS S C, LUMSDEN A C, BARTON N R. Fundamentals of rock joint deformation[J]. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 1983,20:249-268.

[26] KUMAR A. The effect of stress rate and temperature on the strength of basalt and granite[J]. Geophysics, 1968,33(3):501-510.