基于TCK模型的非贯通节理岩体动态损伤本构模型*

刘红岩，杨 艳，李俊峰，张力民

(1.中国地质大学(北京)工程技术学院，北京 100083；2.西藏大学工学院，西藏 拉萨 850000；3.黄淮学院建筑工程学院，河南 驻马店 463000；4.北京科技大学土木与环境工程学院，北京 100083；5.河北承德钢铁公司，河北 承德 067000)

基于TCK模型的非贯通节理岩体动态损伤本构模型*

刘红岩1,2，杨 艳3，李俊峰1，张力民4,5

(1.中国地质大学(北京)工程技术学院，北京 100083；2.西藏大学工学院，西藏 拉萨 850000；3.黄淮学院建筑工程学院，河南 驻马店 463000；4.北京科技大学土木与环境工程学院，北京 100083；5.河北承德钢铁公司，河北 承德 067000)

提出在岩体动态损伤本构模型中应同时考虑宏、细观缺陷；基于能量原理和断裂力学理论推导得出了同时考虑节理几何及力学特征的宏观损伤变量(张量)的计算公式；基于综合考虑宏、细观缺陷的复合损伤变量(张量)及完整岩石动态损伤Taylor-Chen-Kuszmaul(TCK)模型，建立了相应的单轴压缩下节理岩体动态损伤本构模型；利用该模型讨论了节理内摩擦角及节理长度对岩体动态力学特性的影响规律。研究表明，试件动态峰值强度随着节理内摩擦角的增大而增大，随着节理长度的增加而减小。

固体力学；动态损伤本构模型；非贯通节理岩体；宏观缺陷；细观缺陷；损伤耦合；应力强度因子

工程岩体都不可避免地含有节理、裂隙等宏观缺陷和微裂隙、微孔洞等细观缺陷，这2类不同尺度缺陷分别以不同的作用机理影响着岩体的动态力学特性。早期人们通常以宏观完整岩石为研究对象，仅考虑微裂纹等细观缺陷的影响，如B.Budiamsky等[1]采用细观力学对岩石中弥散分布的微裂纹群进行了分析，建立了相应的岩石爆破损伤本构模型。D.E.Grady等[2]认为岩石中的原生微裂纹服从双参数Weibull分布，在外载下，其中一些微裂纹被激活并扩展，由此提出了岩石爆破的GK损伤模型。L.M.Taylor[3]引进B.Budianshy等[1]的有效体积模量和泊松比与微裂纹密度的关系表达式和D.E.Grady等[2]给出的碎块尺寸表达式，最终建立了Taylor-Chen-Kuszmaul(TCK)模型。然而由于仅考虑细观缺陷的岩石损伤本构模型中无法考虑节理等宏观缺陷对岩体动态力学特性的影响。而如何描述节理对岩体力学特性的影响则是目前岩体力学研究中的重点和难点问题。T.Kyoya等[4]借助于损伤理论对该问题进行研究，以节理的长度、倾角、间距及组数等几何特征来定义相应的损伤张量。但是该方法种仅考虑了节理几何性质如节理长度及倾角的影响，而未考虑节理力学性质如抗剪强度的影响，也就是说这种定义方法中认为损伤是无法传递应力的，这与拉伸荷载下的岩体力学特性较为接近，而与压缩荷载下的岩体力学特性相差甚远。这主要是由于在压缩荷载下，节理面两侧的岩体将在节理面处发生闭合及摩擦滑动进而沿节理尖端发生扩展，此时节理面将能够传递部分压应力和剪应力，而且传递系数的大小也与节理面的抗剪强度参数(如摩擦角和粘聚力等)密切相关。为此，不少学者又采用不同的方法对仅考虑节理几何特性的损伤模型进行修正，如T.Kawamoto等[5]通过引入节理传压及传剪系数来考虑压缩荷载下节理能够传递部分压应力与剪应力的特点以对上述模型进行修正，但如何准确确定这两个系数又成为新的难题。同时关于宏、细观两类缺陷对岩体动力学特性的影响问题，张力民等[6]、H.Y.Liu等[7]基于损伤理论建立了考虑宏细观缺陷耦合的复合损伤变量(张量)，进而建立了相应的节理岩体动态损伤本构模型，但是其研究对象仅限于贯通节理岩体。

综合目前非贯通节理岩体的动态损伤本构模型研究现状可知，有以下2个问题亟待解决：(1)目前的岩体宏观损伤张量定义没有很好地同时考虑节理几何及力学参数的共同影响，因而需要提出新的节理岩体宏观损伤变量定义方法；(2)目前学术界仅提出了考虑宏细观缺陷耦合的贯通节理岩体动态损伤本构模型，而针对非贯通节理岩体的研究则较为少见。

为此，本文中首先在含非贯通节理岩体在单轴压缩荷载下的扩展机理进行研究的基础上，提出同时考虑节理几何及力学参数的宏观损伤变量计算方法；其次，在张力民等[6]和H.Y.Liu等[7]研究的基础上，提出基于完整岩石TCK模型的非贯通节理岩体单轴压缩动态损伤本构模型；采用该模型对单条非贯通闭合节理的岩体进行模拟，分析其力学特性。

1 考虑宏观缺陷的非贯通节理岩体损伤本构模型

1.1 非贯通节理岩体损伤模型的建立

根据断裂力学，对于平面应变问题(仅存在Ⅰ、Ⅱ型裂纹，不存在Ⅲ型裂纹)，弹性体因节理存在而引起的附加应变能增加U1为[8]：

(1)

式中：E和ν分别为弹性体的弹性模量和泊松比，G为能量释放率，KⅠ、KⅡ分别为节理尖端的Ⅰ、Ⅱ型应力强度因子，A为节理表面积。

对于单个节理，A=Ba(单边节理)或2Ba(中心节理)；对于多个节理，A=NBa(单边节理)或2NBa(中心节理)，其中：N为节理个数，B为节理深度，a为节理半长。

在单轴应力σ下，损伤应变能释放率为：

(2)

式中：D为节理对岩体造成的损伤变量。

令UE表示与应力σ对应的单位体积弹性应变能，在单轴应力状态下该弹性应变能可写为：

(3)

把式(2)代入式(3)可得：

(4)

当岩体内不含节理时，则D=0，此时式(4)可写为：

(5)

因节理存在而引起的单位体积弹性应变能改变量为：

(6)

假设研究对象的体积为V，弹性体因节理存在而引起的弹性应变能改变量为：

(7)

式(7)中的ΔUE和式(1)中的U1都是由节理存在引起的弹性应变能改变量，二者应相等，即:

ΔUE=U1

(8)

由式(1)、(7)～(8)可得：

(9)

下面对节理岩体进行受力分析，求出KⅠ和KⅡ的表达式。

1.2 应力强度因子计算

1.2.1 单条非贯通节理岩体的应力强度因子

图1 翼裂纹扩展模型示意图Fig.1 Sketch of wing crack growth model

如图1所示，节理与水平面的夹角为α，上部岩块在压缩荷载下将会产生沿节理面的滑动趋势。若沿节理面的剪应力超过其摩擦力时，试件将沿节理面发生摩擦滑移，进而引起翼裂纹由节理尖端大约以θ=70.5°的方向扩展[9-10]。此时节理面上的正应力σα和切应力τα分别为：

σα=σcos2α

(10)

(11)

式中：α为节理倾角。

设节理面的摩擦角为φ，则其摩擦因数μ=tanφ。由于节理面上的黏聚力相对于摩擦力来说要小得多，因此可以忽略不计，故由式(10)～(11)可得节理面上的滑移驱动力为：

(12)

在单轴压缩下，参照文献[11]并考虑翼裂纹的扩展方向，节理尖端翼裂纹Ⅰ和Ⅱ型应力强度因子KⅠ和KⅡ可修改为：

(13)

式中：a为节理半长；l为翼裂纹扩展长度；引入l*=0.27a，使l=0时，KI、KⅡ非奇异，θ为节理尖端翼裂纹扩展角。

考虑翼裂纹即将扩展的临界状态，即当l=0时，翼裂纹应力强度因子KⅠ、KⅡ可修正为：

(14)

图2 非贯通裂隙岩体模型Fig.2 Model of intermittently cracked rockmass

翼裂纹长度l=0所对应的状态即为非贯通节理岩体未发生扩展时的初始状态，若求出此时的节理尖端应力强度因子，代入式(9)即可得到岩体由于含有初始非贯通节理时的初始损伤变量，进而通过张量化得到相应的初始损伤张量。可以看出由该方法求出的损伤张量可同时考虑节理的几何及力学特性，由此求出的节理岩体损伤模型也将更符合实际情况。

1.2.2 单排及多排非贯通平行节理的应力强度因子

图2为多排非贯通裂隙岩体模型示意图，其中，b为同排相邻节理中心间距，d为排间距。

对于单排节理的岩体，考虑节理间相互作用，岩体的有效应力强度因子为[12]：

(15)

式中：KΙ0、KⅡ0分别为单个Ⅰ、Ⅱ型节理的应力强度因子；φ为非贯通节理的连通率，φ=2a/b。

若岩体含有多排非贯通节理时，则其有效应力强度因子为[12]：

(16)

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

当岩体含单组单排节理时，把式(10)～(12)、(14)、(15)代入式(9)可得：

(17)

而当岩体含单组多排节理时，把式(10)～(12)、(14)、(16)代入式(9)可得：

(18)

2 考虑宏细观缺陷耦合的非贯通节理岩体单轴压缩动态损伤本构模型

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

由于岩体中的宏、细观两类缺陷均对岩体的动力学特性有着重要影响，因此在采用损伤理论对岩体的动力学特性进行研究时，应建立能够同时考虑上述2类不同尺度缺陷的损伤变量(张量)计算方法。杨更社等[13]、张力民等[6]已对此进行了较为深入的研究，下面就张力民等[6]的研究成果，给出如下考虑宏细观缺陷耦合的损伤变量(张量)计算公式：

(19)

式中：Ω12为上述2种不同尺度缺陷引起的耦合损伤变量；D为微裂纹等细观缺陷引起的损伤变量；Ω为节理等宏观缺陷引起的损伤张量，Ω=DN，其中N是一个二阶对称张量，其计算方法为：

(1)对于含单组非贯通节理的岩体，假定非贯通节理的法向矢量n与x轴夹角为β，则有：

(20)

(2)对于含2组以上非贯通节理的岩体，N的计算方法是：设岩体中有I组非贯通节理，其单位法向矢量分别为n(i)(i=1,2，…，I)，则

(21)

2.2 考虑宏细观缺陷耦合的非贯通节理岩体动态损伤本构模型

由前述可知，L.M.Taylor等[3]已经建立了仅考虑微裂纹等细观缺陷的岩石动态损伤本构模型，即

(22)

那么由损伤理论可知，当岩体内同时含有宏、细观缺陷时，在损伤模型中则体现为损伤变量的变化。因此把式(22)中的损伤变量ω，用式(19)所示的耦合损伤变量Ω12替换即可，而其中的细观损伤变量D则相应地用ω来替换，即可得到非贯通节理岩体的单轴动态压缩损伤本构方程：

(23)

式中：σ、ε和E分别为张量形式的应力、应变和弹性模量。

3 算例分析

3.1 含单条非贯通闭合节理的岩体力学特性

这里借鉴L.M.Taylor等[3]的实验资料，取宏观完整岩石的动态杨氏模量、泊松比、密度及应变率分别为10.8 GPa、0.2、2 270 kg/m3和102s-1时，k和m可分别取为5.12×1022m-3和7。

(24)

模型上部受到垂直向下的动荷载，下部固定。按平面应变问题分析。

由式(23)计算可得岩体单轴压缩动态应力应变曲线如图4所示。可以看出：当有非贯通节理等宏观缺陷存在时，岩体动态应变峰值强度减小。对本算例而言，非贯通节理岩体的动态峰值强度为31.89MPa，约为相应完整岩石的97.5%，这说明非贯通节理的存在削弱了岩体峰值强度，即表明宏、细观2种缺陷对岩体的强度均有影响。

图3 计算模型Fig.3 Calculation model

图4 岩体单轴压缩动态应力应变计算曲线Fig.4 Dynamic stress-strain calculation curve of rock axial compression

3.2 节理内摩擦角不同时的岩体动力学特性

图5 不同节理内摩擦角的试件动态应力应变曲线Fig.5 Dynamic stress-strain curves of the samples with different joint internal friction angles

计算模型仍为图3所示的非贯通节理试件，载荷应变率为100 s-1，节理长度为4 cm，取节理内摩擦角(φ)分别为：0°、15°、30°，其余条件及计算参数同上。图5为不同节理内摩擦角的试件动态应力应变曲线。由图5可知：岩体动态应力应变曲线斜率即岩体动态弹性模量随着节理内摩擦角的增加增加，即节理岩体刚度增加；同时随着节理内摩擦角的增加，试件动态峰值强度也逐渐增加，3种不同工况下的试件动态峰值强度分别为29.85、31.89和33.98 MPa。这是因为随着节理内摩擦角的增加，节理抗剪强度提高，相应地由节理对试件造成的宏观损伤降低，进而导致试件的复合损伤程度降低，强度相应提高。

3.3 节理长度不同时的岩体动力学特性

图6 不同节理长度的试件动态应力应变曲线Fig.6 Dynamic stress-strain curves of the samples with different joint length

计算模型仍为图4所示倾角为45°的非贯通节理试件，载荷应变率为100 s-1，节理内摩擦角为15°，取节理长度分别为2、4和6 cm，其余条件及计算参数同上。图6给出了不同节理长度的时间动态应力应变虚线。由图6可知：岩体动态应力应变曲线的斜率(即岩体动态弹性模量)随着节理长度的增加而减小，即节理岩体刚度降低；同时试件动态峰值强度也随着节理长度增加的而逐渐减小，这是因为随着节理长度增加，由节理对试件造成的宏观损伤增加，进而导致试件的复合损伤程度提高，动态峰值强度降低。当节理长度由2 cm分别增加到4和6 cm时，试件动态峰值强度则由33.79 MPa分别下降到31.89和29.14 MPa，下降幅度分别为5.62%和13.76%，即随着节理长度增加，试件动态峰值强度的下降幅度比较明显。

4 结 语

(1)基于能量原理和断裂力学理论推导出了能够同时考虑节理几何及力学参数的宏观损伤变量(张量)计算公式，它很好地克服了目前大多数岩体损伤变(张)量无法考虑节理力学参数的不足。

(2)基于完整岩石动态损伤TCK模型及宏细观缺陷耦合的观点，建立了同时考虑宏细观缺陷的非贯通节理岩体动态单轴压缩损伤本构模型。

(3)根据模型计算，结果表明非贯通节理的存在将导致岩体强度降低、刚度弱化。随着节理内摩擦角的增加，试件动态峰值强度及弹性模量均随之增加；而随着节理长度增加，则反之。

[1] Budiansky B, O’Connell R J. Elastic moduli of a cracked solid[J]. International Journal of Solids Structures, 1976,12(2):81-97.

[2] 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, 1980,17(3):174-157.

[3] Taylor L M, Chen E-P, Kuszmaul J S. Microcrack induced damage accumulation in brittle rock under dynamic loading[J]. Computer Method in Applied Mechanics and Engineering,1986,55(3):301-320.

[4] Kyoya T, Ichikawa Y, Kawamoto T. Damage mechanics theory for discontinuous rock mass[C]∥Proceedings of the 5th International Conference on Numerical Methods in Geomechanics. Nagoya, Japan, 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 of Numerical Analysis Method in Geomechanics, 1988,12(1):1-30.

[6] 张力民,吕淑然,刘红岩.综合考虑宏细观缺陷的岩体动态损伤本构模型[J].爆炸与冲击,2015,35(3):428-436. Zhang Limin, Lü Shuran, Liu Hongyan. A dynamic damage constitutive model of rock mass by comprehensively considering macroscopic and mesoscopic flaws[J]. Explosion and Shock Waves, 2015,35(3):428-436.

[7] Liu H Y, Lv S R, Zhang L M, et al. A dynamic damage constitutive model for a rock mass with persistent joints[J]. International Journal of Rock Mechanics & Mining Sciences, 2015,75(4):132-139.

[8] 李世愚,和泰名,尹祥础.岩石断裂力学导论[M].合肥：中国科学技术大学出版社,2010:89-98.

[9] 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.

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

[11] Lee S, Ravichandran G. Crack initiation in brittle solids under multiaxial compression[J]. Engineering Fracture Mechanics, 2003,70(13):1645-1658.

[12] 李建林,哈秋瓴.节理岩体拉剪断裂与强度研究[J].岩石力学与工程学报,1998,17(3)：259-266. Li Jianlin, Ha Qiuling. A study of tensile-shear crack and strength related to jointed rock mass[J]. Chinese Journal of Rock Mechanics and Engineering, 1998,17(3):259-266.

[13] 杨更社,谢定义.岩体宏观细观损伤的耦合计算分析[C]∥第六次全国岩石力学与工程学术大会论文集.湖北武汉，2000：327-329.

(责任编辑 王小飞)

Dynamic damage constitutive model for rock mass with non-persistent joints based on the TCK model

Liu Hongyan1,2, Yang Yan3, Li Junfeng1, Zhang Limin4,5

(1.CollegeofEngineering&Technology,ChinaUniversityofGeosciences(Beijing),Beijing100083,China;2.SchoolofEngineering,TibetUniversity,Lhasa850000,Xizang,China;3.ArchitectureEngineeringCollege,HuanghuaiUniversity,Zhumadian463000,Henan,China; 4.CivilandEnvironmentalEngineeringSchool,UniversityofScienceandTechnologyBeijing,Beijing100083,China;5.HebeiChengdeIronandSteelCorporation,Chengde067002,Hebei,China)

This paper proposes that both macroscopic and mesoscopic flaws should be considered in the dynamic damage constitutive model for the jointed rock mass. Firstly, the calculation formula of the macroscopic damage variable (tensor) of the jointed rock mass is deduced based on the energy principle and the fracture mechanics theory so that the geometrical and mechanical parameters can be considered at the same time. Secondly, the compound damage variable (tensor) including both macroscopic and mesoscopic flaws and the Taylor-Chen-Kuszmaul model (TCK) for the intact rock are adopted to establish the corresponding dynamic damage constitutive model for the jointed rock mass under uniaxial compression. Finally, the effect law of the joint internal friction angle and the joint length on rock mass's dynamic mechanical property is examined with this model. The results show that the dynamic climax strength of the samples increases and decreases respectively with the increase of the joint internal friction and the joint length.

solid mechanics; dynamic damage constitutive model; rock mass with non-persistent joints; macroscopic flaw; mesoscopic flaw; damage coupling; stress intensity factor

10.11883/1001-1455(2016)03-0319-07

2014-09-22；

2015-02-22

国家自然科学基金项目(41002113，41162009)； 中央高校基本科研业务费专项基金项目(2652014019，2652015263)

刘红岩(1975- )，男，博士，教授，lhyan1204@126.com。

O341国标学科代码：13015

A