基于复合损伤的节理岩体动态本构模型研究

邓正定，王 桢，刘红岩

（1.中国铁道科学研究院，北京 100081；2.中铁西北科学研究院有限公司，甘肃 兰州 730000；3.中国地质大学（北京）工程技术学院，北京 100083）

1 引言

自然界的岩体中存在大量的节理、层理等不连续面，这些含有不连续结构面的岩体称为节理岩体。节理岩体的定性不仅与完整岩石的强度有关，而且还受岩体内发育的节理面所控制。宏观上节理岩体由于节理面的方位、倾角不同而呈现出各向异性的特点，同时也是非均匀分布的。另外，节理岩体内部还存在许多微小的裂纹分布，这是节理岩体细观上的缺陷，由于岩体存在这些宏观和微观的缺陷，如何建立有效的节理岩体动态本构模型显得较为困难。目前，对完整岩石的动态本构模型已经有了广泛的研究，对含微观节理的动态本构模型也有一定研究，但对于存在各向异性的宏观节理的岩体在高应变率下的动态本构模型研究比较少，如Seamen等[1－2]提出NAG-FRAG模型，该模型利用显微镜观测方法，基于统计学的原理知识，建立了岩体的动态本构模型。Burlion等[3]提出一个损伤与塑性耦合的本构模型来描述混凝土、岩石等脆性材料的动态破坏，提出压缩的弱化效应由塑性应变控制，拉伸损伤由弹性应变所控制，将压缩损伤和塑性应变耦合起来。刘文韬等[4]基于连续介质力学的观点提出了能够描述岩石塑性变形的连续损伤本构模型。单仁亮等[5]对大理岩和花岗岩的本构特性进行深入的研究，根据其应力-应变曲线的特点，考虑损伤的影响将统计损伤模型和黏弹性体结合，建立冲击时效损伤模型，利用模型对大理岩和花岗岩进行曲线拟合。但这些研究都没有充分考虑节理岩体中起控制性作用的宏观节理及各向异性的特点，因存在一定的偏差，需要对模型进行改进，而实际工程中岩体常处于地震、爆破及开挖扰动等动力荷载作用之下，岩体的破坏失稳与节理岩体的动态力学特性密切相关[6]。因此，对高应变率下节理岩体的本构关系研究具有巨大的实际工程应用价值。

在前人研究基础上，运用分离式霍普金森压杆装置（SHPB）对几种典型几何特征的节理岩体在不同高应变率下进行了大量试验，综合考虑节理几何特征、应变率效应对试验结果进行合理分析的基础上，基于复合损伤的分析方法，构造一个能够描述节理岩体动态响应的复合损伤型本构模型，模型中各个参数都有明确的物理意义，以便于工程实际中应用。

2 节理岩体动态复合损伤分析

2.1 损伤张量

岩体中存在对岩体力学特性起控制性作用的节理面，研究的主要难点是节理岩体的节理面是各向异性的、不均匀的。由于节理面的存在，使得节理岩体的力学特征与完整岩石存在非常大的差异。因此，建立能够有效描述节理岩体动态力学特性的本构关系模型，最主要的是考虑节理各向异性的特点。

节理等宏观缺陷的存在将导致岩体力学性质的弱化及各向异性，为了反映这种影响，目前常采用损伤力学的方法进行描述，假定损伤后的岩体本构关系仍服从虎克定律[7]，则节理对岩体的损伤则体现为弹性常数的弱化，即节理岩体的弹性常数与损伤张量之间的关系，可以表示为

式中：E0、E 分别为完整岩块和节理岩体的弹性张量；I为单位张量；Ω为节理岩体损伤张量。

下面以平面问题为例，讨论节理岩体损伤张量的计算。目前，常用二阶张量描述节理岩体的各向异性损伤，假如节理岩体中存在n 组节理，则节理岩体的初始损伤张量为

式中：V为岩体体积；l为节理平均间距(节理厚度），n为试件中节理组数；mi为样本中第i 条节理表面上的单位法向矢量；ai为样本中第i 条节理表面积，⊗指前、后两个矢量的张量积。mi⊗mi为二阶张量。

在自然界的节理岩体中，节理间并不是完全没有物质的，而是往往夹杂着一定的充填物，通过试验可知，不同的充填物对节理岩体的损伤有一定的影响，为此这里引入充填物影响系数c，式（2）可以改为

式中：c为充填物系数。

由此可知，初始损伤张量可以很好地反映出节理几何特征对节理岩体损伤的影响。

现假设有一条节理与y 方向呈θ 角，见图1，则节理的外法线张量为

式中：θ为节理面与y 方向夹角。

图1 典型节理面Fig.1 A typical joint surface

则初始损伤张量为

当有n 条平行节理时，式（5）可推广为

2.2 损伤变量

前面所述损伤张量也称初始损伤张量，从式中可以看出，该损伤张量可以很好地表示节理岩体各向异性损伤的特点，但岩体在加载的过程中，裂纹会继续扩展，而该损伤可认为是各向同性的。针对该损伤引入损伤因子D，显然，该损伤因子D 是与应变和应变率有关的变量。利用这一方法可以得到材料的损伤型本构关系，但最重要的还在于找到损伤变量与本构状态参量之间的函数关系，即损伤演化规律。有的损伤模型将裂纹成核密度作为损伤演化的参数，但裂纹密度难以定量描述；有的损伤模型将某一特定应变作为损伤演化的阀值，但它没有考虑应变率对损伤演化的贡献。其实，不管宏观损伤还是微观损伤，应变可对损伤有强化作用，应变率对损伤有弱化作用，即共同影响损伤的发展，因此，将损伤演化方程的一般函数形式表示为

对式（7）作无量纲化处理，则有

由于在动态荷载加载之前，因为有节理的存在，节理岩体已经存在损伤。为了处理方便，将其定量用0D 表示，则损伤演化方程表示为

式中：D0为初始损伤值；εth为损伤阀值应变，通常取最大应变的0.2～0.4倍；a 和b 均为待定系数，用于描述试验所得损伤演化曲线的形状；KD为损伤变量系数。需要指出的是，这里的损伤变量是指岩体初始微细观裂纹以及在加载过程中的裂纹扩展等造成的损伤，并不是前面所述的宏观节理而造成的初始损伤，而初始微细观裂纹相对宏观节理可忽略不计，因此，可取D0=0 。

2.3 一维动荷载下的有效应力

由前述可知，节理岩体同时存在宏微观两种损伤，二者都起到了弱化岩体刚度和强度的作用，因此，在节理岩体力学分析中应综合考虑两种不同损伤的共同作用[7]。由于岩石中的微损伤是各向同性，而宏观节理、裂隙造成的损伤却是明显各向异性的，由式（1）可知，E0为完整无损岩石的弹性模量，结合式（9）可以得到综合考虑宏微观两种损伤的节理岩体本构模型为

式中：σ为应力张量；ε为应变张量。

用二阶损伤张量描述岩体的节理裂隙时，针对岩体受压应力情况，考虑到节理裂隙能够传递部分压应力和剪应力的特点，引入下列有效应力张量：

式中：λ、η 均为与各参数相关的量。

实际上，在一维动荷载下，因为不考虑围压的影响，故只需考虑法线方向与x 轴方向相同的应力即可，故只需取下式：

3 节理岩体动态本构模型

3.1 本构关系一般描述

一般来说，任何材料本构模型可以描述为应力、应变、应变率和温度等物理、力学参量之间的函数关系，其形式可表示为

式中：σ为表观应力张量，表示含损伤的材料应力。

大量研究结果表明，黏性效应和损伤滞后是产生岩石材料非线性和率敏感性的两个主要原因。如果不考虑温度效应，岩体类材料本构方程一般可表示为

当D=0时，表示材料无损伤；当D=1时，表示材料完全丧失承载能力。

由于节理岩体是有初始损伤的，并是各向异性的，并且初始损伤是以张量的形式表示的，因此，上式都用张量表示：

式中：σi为材料无损伤时的应力张量；I为单位张量；损伤变量D是应变ε、应变率的函数。

在一维冲击应力下的节理岩体本构关系模型，即不考虑围压的作用，所以式（16）的应力-应变都只考虑加载方向的变化情况，而损伤张量也可以只考虑加载方向，只需方向对应即可，这样可以减少很多不必要的运算。

3.2 损伤型节理岩体动态本构模型

Kinoshita等[8]在进行试验研究的基础上认为，岩石的动力学本构特性可以用Bingham模型来描述，见图2，该模型通常又称为过应力模型。其本构方程为

图2 Bingham模型Fig.2 Bingham model

原始的过应力模型本构方程认为，E0、S0、S分别为静态压缩强度后的弹性模量、静态压缩破坏强度和静态压缩强度后的应力，经过修正的过应力模型认为，E0对应于动载应力-应变曲线线弹性部分的斜率，S0和S一样，都为动载应力-应变曲线线弹性部分的极限[9]，α、β 均为不同岩体材料的固有常数。

于亚伦等[10]曾用修正的过应力模型对几种磁铁矿和花岗岩进行曲线拟合，得到了相应的过应力模型本构方程，使得修正后的过应力模型得到很好的应用。但修正后的过应力模型存在明显的缺点，即不能够完全反映弹性模量随应变率变化的特性。同时，对式（17）中的σ ≥ S0的情况，进行积分转换会发现，修正的过应力模型使用范围有局限性，其积分转换得到：

进行移项简化后得到

节理岩体在加载过程中的损伤，主要表现在弹性模量的变化，而加载过程中动态弹性模量是因节理岩体动态损伤引起，即随损伤变量D而改变，其关系如下：

联立式（9）、（20）、（21）可得

式中：σth、εth分别指应力阀值和应变阀值，通常取峰值的0.2～0.4倍。

式（21）中考虑了加载过程中损伤的影响，同时考虑边界条件ε=εth时，S=σth即可。

式（21）仅仅考虑了加载过程中的动态损伤，并没有考虑节理岩体的初始损伤，故用节理岩体初始损伤的有效应力 σ*来代替σ，式（21）可变为

由式（22）可以看出，除去试验可得参数之外，还有α、β、KD、a、b、c这6个参数是要通过数据拟合得到的，通过观察发现，一些参数可以合并成一个参数，具体如下：

式中：K为损伤系数；A为应变率损伤指数；B为应变损伤指数；c为充填物系数，可得最终简化的节理岩体动态损伤本构关系模型为

4 试验验证及分析

4.1 试验方案

运用SHPB装置，通过改变节理岩体试件的节理面倾角、节理贯通程度、节理组数、节理厚度、节理间填充物、长径比及加载率对试件进行冲击荷载试验，综合对比分析各因素的影响规律。试验方案按不同因素对试样动荷载力学特性影响可归纳为表1。

表1 SHPB试验方案设计Table 1 Scheme of SHPB tests

4.2 试验与理论比较分析

通过SHPB试验数据来确定节理岩体动态损伤本构模型待定参数。冲击荷载下岩体材料的应变率范围在(101～105) s-1之间，取其加载段应变率的平均值为试验响应应变率，这样试验近似可以看作是恒应变率的，式（26）中包含了13个材料参数，每个参数都有明确的物理意义。其中σ、E0、l、V、ε、θ、ai、可以通过试验数据或通过节理岩体的几何形状取值确定，损伤阀值εth取为最大应变的0.3倍，K、A、B、c 通过数值拟合来确定。表2为通过数值拟合得到最优模型拟合参数。

表2 最优模型拟合参数Table 2 Fitting parameters for the optimal model

利用表2给出的模型参数计算得到了模型预测曲线，并与SHPB试验所测试得到的曲线进行比较。图3～5给出了部分试验测试曲线和模型预测曲线的比较，其中实线为理论预测曲线，虚线为试验曲线。

从图3～5中不难发现，该模型可较好地描述节理岩体动荷载下初始弹性变形阶段、稳态塑性变形阶段和加速变形破坏阶段。随着冲击速度的提高，节理岩体材料内部产生了较大的塑性变形，模型曲线与实测结果相比略有偏低，这表明损伤型黏弹性本构模型在描述强冲击条件下节理岩体的动态力学特性方面略有不足。

图3 不同倾角节理岩体试验与理论预测曲线比较Fig.3 Comparison between the theoretical and experimental curves of rock masses with joints of different inclination angles

图4 不同贯通度节理岩体试验与理论预测曲线比较Fig.4 Comparison oftestand theoretical curves at different transfixion degrees of joint rock mass

图5 不同应变率节理岩体试验和理论预测曲线比较Fig.5 Comparison of theoretical and experimental curves of joint rock masses at different strain rates

在SHPB试验中，由于受到试验技术条件和数据选取主观性的限制，动态测试曲线与本构模型曲线在塑形变形阶段存在一定的偏差，这也造成了模型预测和试验结果的不符。通常认为，测试曲线中的稳定发展阶段为材料的实体弹性模量，然而相近应变率下，实体弹性模量的离散性仍然较大，这也造成了模型曲线与试验曲线在稳定发展阶段的偏差。不过，由于该模型参数较少，且相对容易确定，因而在工程实践中仍有较大的应用价值。

5 结论

（1）基于复合损伤力学的观点，引入节理岩体各向异性损伤张量，加载时的损伤变量，综合考虑不同应变率、不同几何特征等多种因素研究了节理岩体的动态力学特性，构造了一个能够描述节理岩体材料在冲击荷载下动态力学响应的损伤型黏弹性本构模型。

（2）由于该本构模型所取得参数是采用SHPB试验装置对节理岩体动态加载试验所得，所以建立的节理岩体动态本构关系模型具有较高的工程实用价值。

（3）本文提出的本构模型中需确定的材料参数较少，且不需要很复杂的试验即可确定。

（4）本文试验中考虑了长径比的效应，发现长径比对强度特性有一定影响影响，建立的节理岩体动态本构模型中并没有考虑长径比的因素，建议后续加强这方面的研究。

（5）在建立的节理岩体动态本构模型中，没有考虑围压的影响，而自然界的岩体大多处于有围压状态，建议在后续的节理岩体本构关系研究中，考虑围压的因素，建立更符合工程实际的节理岩体动态本构模型。

[1]SEAMEN L，CURRAN D R，SHOCKEY D A.Computational models for ductile and brittle fracture[J].Journal of Applied Physics，1976，47(11):4814－4826.

[2]SEAMEN L，CURRAN D R，MURRI W J.A continuum model for dynamic tensile micro-fracture and fragmentation[J].Journal of Applied Mechanics ——Transactions of the ASME，1985，52(3):593－600.

[3]BURLION N，GATUINGT F，GILLES P C，et al.Compaction and tensile damage in concrete:Constitutive modeling and application to dynamics[J].Computers Methods in Applied Mechanics and Engineering，2000，183(3):291－308.

[4]刘文韬，王肖钧，周钟.岩石的连续损伤本构模型及在地下爆炸波数值计算中的应用[J].岩石力学与工程学报，2004，23(13):2149－2156.LIU Wen-tao，WANG Xiao-jun，ZHOU Zhong.Continuously damaged constitutive model of rock and its application to numerical simulation for underground strong explosion[J].Chinese Journal of Rock Mechanics and Engineering，2004，23(13):2149－2156.

[5]单仁亮，薛友松，张倩.岩石动态破坏的时效损伤本构模型[J].岩石力学与工程学报，2003，22(11):1771－1776.SHAN Ren-liang，XUE You-song，ZHANG Qian.Time dependent damage model of rock under dynamic loading[J].Chinese Journal of Rock Mechanics and Engineering，2003，22(11):1771－1776.

[6]张平，李宁，贺若兰.不同应变速率下非贯通裂隙介质的力学特性研究[J].岩土工程学报，2006，28(6):750－755.ZHANG Ping，LI Ning，HE Ruo-lan.Mechanical properties of fractured media containing intermittent fractures at different strain rates[J].Chinese Journal of Geotechnical Engineering，2006，28(6):750－755.

[7]张吉宏，刘红岩.综合考虑宏微观复合损伤的节理岩体本构模型[J].煤田地质与勘探，2013，41(6):49－52.ZHANG Ji-hong，LIU Hong-yan.Constitutive model of jointed rock mass by combining macroscopic and microcopic composite damage[J].Coal Geology and Exploration，2013，41(6):49－52.

[8]KINOSHITA S，SATO K，KAWAKITA M.On the mechanical behavior of rocks under impulsive loading[J].Bulletin of the Faculty of Engineering，Hokkaido University，1977，(83):51－62.

[9]谢理想，赵光明，孟祥瑞.岩石在冲击载荷下的过应力本构模型研究[J].岩石力学与工程学报，2013，32(增刊1):2772－2781.XIE Li-xiang，ZHAO Guang-ming，MENG Xiang-rui.Research on excess stress constitutive model of rock under impact load[J].Chinese Journal of Rock Mechanics and Engineering，2013，32(Supp.1):2772－2781.

[10]于亚伦，金科学.高应变率下的矿岩特性研究[J].爆炸与冲击，1990，10(3):266－271.YU Ya-lun，JIN Ke-xue.Studies on the dynamic behavior of magnetites and granite under high strain rates[J].Explosion and Shock Waves，1990，10(3):266－271.

[11]良昌玉，李晓，李守定，等.岩石静态和准动态加载应变率的界限值研究[J].岩石力学与工程学报，2012，31(6):1156－1160.LIANG Chang-yu，LI Xiao，LI Shou-ding，et al.Study of strain rates threshold value between static loading and quasi-dynamic loading of rock[J].Chinese Journal of Rock Mechanics and Engineering，2012，31(6):1156－1160.

[12]WONG L N Y，EINSTEIN H H.Crack coalescence in molded gypsum and Carrara marble:Part 1——Macroscopic observations and interpretation[J].Rock Mechanics and Rock Engineering，2009，42(3):475－511.

[13]LEE H，JEON S.A test and numerical study of fracture coalescence in pre-cracked specimens under uniaxial compression[J].International Journal of Solids and Structures，2011，48:979－999.

[14]MOUTHEREAU F，FILLON C，MA K F.Distribution of strain rates in the Taiwan orogenic wedge[J].Earth and Planetary Science Letters，2009，284(3－4):361－385.

[15]王文，朱维申，马海萍，等.不同倾角节理组和锚固效应对岩体特性的影响[J].岩土力学，2013，34(3):887－893.WANG Wen，ZHU Wei-shen，MA Hai-ping.Influence of joints set with different dip angles and anchorage effect on rock mass behaviour[J].Rock and Soil Mechanics，2013，34(3):887－893.

[16]孙旭曙，李建林，王乐华.单一预制节理试件各向异性力学特性试验研究[J].岩土力学，2014，35(增刊1):29－34.SUN Xu-shu，LI Jian-lin，WANG Le-hua.Experimental research on anisotropic mechanical characteristic of samples with single prefabricated joint[J].Rock and Soil Mechanics，2014，35(Supp.1):29－34.

[17]翟越.岩石类材料的动态性能研究[博士学位论文D].西安:长安大学，2008.

[18]崔少东.岩石力学参数的时效性及非定常流变本构模型研究[博士学位论文D].北京:北京交通大学，2010.