岩石锚杆界面滑移破坏的DEM数值模拟

，，,3

(1.武汉大学 土木建筑工程学院，武汉 430072;2.武汉人才市场有限公司，武汉 430000； 3.同济大学 地下建筑与工程系，上海 200092)

1 研究背景

岩体锚固技术是岩土工程的重要分支，在交通工程、矿业工程、隧道工程、水利工程等领域被广泛运用。根据锚固方式，锚杆可分为粘结型锚杆和机械型锚杆[1]。岩石粘结式锚杆主要依靠灌浆材料(砂浆、水泥、树脂等)产生的粘结力使锚杆、锚固体以及岩体结合成为一个整体。在岩体变形过程中，锚杆通过锚固体与岩体产生协调变形，利用锚杆的高强度，从而加强岩体的薄弱环节，提高岩体的整体承载力。尤春安等[2]认为锚固系统失效的主要形式是锚固体与岩体之间的界面产生滑动破坏，因此，开展岩石粘结锚杆的界面力学特性和破坏机制研究，对于正确指导工程实践具有重要意义。

锚杆拉拔试验是揭示锚固体荷载传递及破坏机制的有效手段。国内外许多研究者开展了锚杆拉拔试验研究，取得了较为丰富的成果[3-4]。在荷载传递方面：张永兴等[5]对普通锚杆进行了室内拉拔试验，指出界面剪应力沿锚固段呈逐渐衰减分布； Farmer[6]认为在较小张拉荷载作用下，锚杆轴力和界面剪应力从荷载作用点到锚杆里端以指数形式衰减；杨松林等[7]认为锚杆受到张拉荷载时，传力深度远小于实际工程中锚杆的设置长度。在锚固强度和破坏机制方面：陈妙峰等[8]指出锚杆的抗拔性能很大程度上受粘结材料强度和锚杆几何尺寸的影响；张发明等[9]建议采用综合因素法和岩体BQ法确定界面粘结强度；朱焕春等[10]认为锚固系统的破坏过程首先是锚杆粘结介质的横向破坏，继而是界面的剪切破坏，并且锚固体应力不断向深部传递；Zhang等[11]通过分段分析锚杆上的剪应力分布，得出了锚杆抗拔极限荷载的理论公式。

尽管试验研究取得了不少研究成果，但受到现阶段试验及观测手段的限制，锚固段的应力分布和破坏过程难以准确获得，因此，数值模拟已成为一种重要的研究手段。在岩体锚固的数值模拟方面，过去大多是基于连续介质力学法[12-13]，无法真实模拟锚固系统的变形、破坏全过程。近年来，非连续介质分析方法发展迅速，颗粒离散单元法(DEM)能有效地模拟岩石材料的宏、微观力学响应及微观变形破坏过程[14]。赵同彬等[15]和张思峰[16]采用了PFC2D软件模拟了岩石锚固段受力变形及破坏的特征，得出的结论与试验规律接近，表明采用颗粒离散元研究锚固问题是可行的。但存在的问题主要有：其一，采用PFC自带的接触模型模拟岩石类材料时存在一些缺陷，比如不能同时满足抗拉和抗压强度[14,17]；其二，所建模型较为粗糙，没有精细考虑岩石、锚杆和砂浆相互作用。因此，难以真实揭示岩石锚杆荷载传递的微观机制。

本文采用DEM方法及笔者提出的岩石微观胶结接触模型[18]，着重从微观破坏角度研究了锚固系统界面破坏的机理。首先建立了包含岩石、砂浆和锚杆的DEM分析模型，模拟了岩石锚杆的拉拔试验，研究了荷载-位移关系，探讨了锚固长度对峰值荷载和界面力学特性的影响，然后分析了拉拔破坏过程中轴力和界面剪应力的变化规律，最后从组构的角度研究了锚杆界面破坏的微观机制。

2 DEM计算模型

2.1 模型建立

图1 DEM试样及加载方式Fig.1 Specimen and loading mode in discrete element model (DEM)

参考张永兴等[5]所做的普通锚室内拉拔试验，建立二维DEM模型，见图1。中间部分为锚杆，锚杆左右两边都是水泥砂浆，砂浆两侧为岩石。其中，锚杆直径为4.8 mm，砂浆宽度6.6 mm，岩石试样宽为200 mm，高为250 mm。岩石上下端的竖向位移采用刚性墙约束，左、右侧为自由边界。其中，为了防止试样出现左右滑移，刚性墙的摩擦系数取经验值1，墙的刚度取1×1010N/m，与试样上部紧密贴合，限制竖直方向上的位移。然后在y方向给锚杆端头施加一个向上的速度v来模拟拉拔过程。

表1 岩石无厚度胶结模型中的微观参数Table 1 Microscopic parameters in zero-thickness bond model for rocks

2.2 试样生成

DEM试样采用蒋明镜等[19]提出的分层欠压法制备。锚固系统模型(岩石、砂浆、锚杆)颗粒直径为0.45～2 mm，均匀分布。试样在制备过程中共分8层，所对应的各层孔隙比分别为：ep(1)=0.215，ep(1～2)=0.213，ep(1～3)=0.211，ep(1～4)=0.209，ep(1～5)=0.207，ep(1～6)=0.205，ep(1～7)=0.203，ep(1～8)=0.20。其中，ep(1～2)表示1～2层综合孔隙比，其他类同。然后，保持左、右墙不动，通过控制上、下墙的相向运动，给试样施加一定的竖向固结压力。压力的大小需根据模拟试样的垂直应力大小试算确定，本文选定为529 kPa。待法向应力伺服稳定后，施加胶结，然后继续循环直至稳定，则制样完成。颗粒微观参数的选取详见2.3节。

2.3 微观接触模型及相关参数选取

颗粒接触模型采用蒋明镜等[18,20]提出的岩石类无厚度胶结模型。该模型是总结室内水泥铝棒胶结试验提炼出来的模型，胶结接触物理模型见图2。

图2 微观胶结接触物理模型[18]Fig.2 Physical model of micro-bond contact laws[18]

对比PFC2D自带的Parallel bond模型，该微观接触模型的最大优势是能够同时满足岩石类材料的宏观抗拉强度和抗压强度，因此能更真实地模拟岩石类材料的力学性质，该模型的详细介绍见文献[17,20-21]。该模型由法向接触模型、切向接触模型和转动接触模型3个部分组成，模型参数共11个，详见表1。DEM模型中岩石、水泥砂浆和锚杆均采用上述胶结接触模型，岩石与砂浆以及砂浆与锚杆界面均采用砂浆粘结强度来模拟。各部分的微观参数以如表2所示的宏观力学参数为目标分别进行标定，其中岩石参数取自文献[22]。标定过程为：先模拟各材料的单轴抗压和单轴抗拉试验，通过调整微观参数，使得DEM试样的宏观力学性质与实际的目标相近，最终确定的微观胶结参数如表1所示。

表2 DEM试样的宏观力学参数Table 2 Macroscopic parameters of DEM specimen

注:“/”左侧代表室内试验数据，右侧代表数值模拟结果

2.4 模拟方案

本次研究共模拟4组不同锚固长度的拉拔试验，分别为100，150，200，250 mm。在室内拉拔试验[5]中，普通锚的直径为4.8 mm，锚固长度为250 mm。为了便于与室内拉拔试验结果对比，按锚杆体表面积等效的原则，取二维模型的等效厚度为7.5 mm。将锚固段按40 mm长度均匀分成若干段,分别统计每段的轴力和界面剪应力。轴力统计方法为计算锚杆每段界面颗粒接触法向力之和，界面剪应力根据每段锚杆上下截面轴力之差换算。

3 模拟结果及分析

3.1 锚杆拉拔破坏形态

DEM模拟的锚杆拉拔试验最终破坏形态如图3所示。由图3可以看出，本文模拟的破坏类型是砂浆与岩石的界面破坏，与室内试验[2]破坏类型一致。

图3 DEM试样破坏形态以及胶结破坏Fig.3 Failure form and broken bonds of DEM specimen

3.2 荷载-杆端位移关系

图4为室内试验和数值模拟的荷载-杆端位移关系的对比图。由图4可以看出，拉拔荷载先随着杆端位移的增加而增大，达到峰值之后迅速跌落，最后达到残余强度。模拟荷载-杆端位移关系曲线的变化趋势与室内试验结果基本一致，但是位移和荷载的数值与试验结果还存在一定差异。主要原因是本文建立的是二维离散元模型，而室内锚杆拉拔试验是三维问题，因此两者在定量上存在差异。

图4 荷载-杆端位移关系Fig.4 Relationship between loading and displacement of nodes

3.3 锚固长度的影响

图5给出了拉拔峰值荷载和界面平均粘结强度随锚固长度的变化曲线。界面平均粘结强度定义为

(1)

式中：τy为界面平均粘结强度(MPa)；P为峰值拉拔荷载(N)；b为锚杆等效厚度，取7.5 mm；l为锚固长度(mm)。

图5 峰值荷载和平均粘结强度随锚固长度的变化Fig.5 Variations of peak load and average cohesive strength against anchoring length

由图5可见：①拉拔峰值荷载随着锚固长度的增加而增加，但是并不完全呈线性关系，即随着锚固长度增大，荷载增加幅度逐渐减小；②界面平均粘结强度随着锚固长度的增加而减小，数值在2.7～3.3 MPa范围内，与前人实测的水泥砂浆粘结强度结果较为接近[23]。

3.4 锚杆轴力分布

以锚固长度250 mm为例,不同杆端位移下锚杆轴力分布如图6,其中δp为峰值荷载对应的杆端位移。由图6可以看出：①当杆端位移未超过δp时，锚杆轴力沿长度呈递减分布，且各测点轴力均随杆端位移增加而增大，这与室内试验结果一致[24]。②当杆端位移超过δp后，各测点轴力逐步出现先减后增的现象，与赵同彬等[15]的研究结果一致。例如，当位移达到2δp，靠近杆端(测点1和测点2)的轴力随杆端位移明显减小，而其他测点的轴力继续增大；当杆端位移达到3δp时，所有测点的轴力均出现明显降低；当杆端位移达到4δp时，各测点轴力降幅减缓，且轴力分布趋于均匀，锚杆达到整体破坏状态。

图6 不同杆端位移下的锚杆轴力分布Fig.6 Distribution of axial force of bolt in the presence of varying displacement of nodes

3.5 锚固段剪应力变化

同样以锚固长度250 mm为例, 不同杆端位移下锚固段剪应力的分布曲线如图7所示。

图7 不同杆端位移下锚固段剪应力分布Fig.7 Distribution of shear stress of anchorage segment in the presence of varying displacement of nodes

由图7可知：当杆端位移为0.3δp和0.6δp时，剪应力沿锚固段长度呈递减分布，且随着杆端位移增加，各位置的剪应力也同步增大，与文献[5]中普通锚的分析结果一致；当杆端位移等于δp时，测点1的剪应力出现下降，测点2—测点6的剪应力则保持上升，表明仅测点1达到破坏状态；当杆端位移达到2δp时，测点1、测点2的剪应力均明显下降，除测点4略有增加外，其余测点的剪应力则变化幅度不大，表明测点1、测点2处于破坏状态，测点3—测点6接近临界状态；当杆端位移达到3δp时，所有测点的剪应力均明显下降，表明整个锚固段达到破坏；当杆端位移达到4δp时，各测点的剪应力降幅减缓且数值趋于稳定，表明锚固段处于残余破坏阶段。可见，剪应力的变化规律表明锚杆的拉拔过程是一个渐进破坏的过程。

图8 不同测点剪应力、法向力与剪切位移的关系Fig.8 Variations of shear stress and normal stress of different measuring points against shear displacement

为了进一步揭示锚固段的界面摩擦特性，得到各测点剪应力与剪切位移的关系，如图8(a)所示。由于试样竖向变形受刚性板限制，界面处的剪切位移可取该测点锚杆的位移。由图8(a)可以看出，各测点剪应力先随剪切位移增加而增大，剪应力达到峰值后迅速跌落至某一残余值。由图8(a)还可发现，不同测点的峰值剪应力和对应的峰值剪切位移并不相同。距离加载端越近，峰值剪应力越大，对应的峰值剪切位移越小。其原因是测点靠近端部时，边界刚性墙的约束作用增大了界面法向压力，因此提高了界面的抗剪强度和剪切刚度。各测点对应的法向压力如图8(b)所示，由图8(b)可以看出离加载端越近，所承受的法向压力越大，不同测点的峰值剪应力与法向压力基本呈正相关关系。这表明，即使在同样的剪切位移下，由于受力状态差异，锚固段各段剪应力发挥也不同步，这也可能是造成锚固段渐进破坏的原因之一。

3.6 锚固段界面破坏过程的微观分析

图9为荷载和胶结破坏数目随杆端位移的变化曲线(以锚固长度为250 mm为例)。本文采用的微观胶结模型中，胶结破坏的类型定义为拉伸胶结破坏和剪扭胶结破坏。从图9中可以看出：当杆端位移很小时，荷载与位移基本呈线性关系，胶结破坏数量很少；随着杆端位移增加，荷载与杆端位移逐渐呈非线性关系，胶结破坏数目增加较快，且以拉伸破坏为主；随着杆端位移继续增加，荷载逐渐达到峰值，拉伸胶结破坏数目增长减慢，而剪扭胶结破坏数目增长加快；在荷载逐渐跌落至残余值的过程中，2种胶结破坏数目均趋于稳定，表明不再产生新的胶结破坏，残余拉拔荷载主要由锚杆和砂浆之间的界面机械咬合力提供。

图9 荷载和胶结破坏数目随杆端位移的变化曲线Fig.9 Variations of loading and broken bond number with displacement of bolt

锚固段微裂纹的扩展过程可采用胶结破坏点组构的玫瑰分布图表示，如图10所示。所谓的胶结破坏点组构，是指该破坏点的颗粒接触矢量方向与x轴正方向的夹角(如图1)，可以用来表征微裂纹的分布。由图10可以看出：当杆端位移较小时(0.3δp)，胶结破坏的组构有3个明显的方向，即0°(180°)、30°(210°)和150°(330°)，说明产生的微裂纹分别垂直于这3个方向；随着杆端位移增加(0.6δp)，0°方向胶结破坏的组构所占比例明显增加，说明锚杆受力方向的微裂缝增多；当杆端位移等于δp时，0°方向胶结破坏的组构所占比例继续增加，已成为胶结破坏的主方向；当杆端位移大于δp后，胶结破坏数增幅变缓，胶结破坏点组构的分布形态变化不大；当杆端位移达到4δp时，胶结破坏点组构的分布已基本稳定。可见，在整个拉拔过程中，锚固段微裂纹逐渐以锚杆受力方向的裂缝扩展为主，微观上表现为垂直于锚杆拉拔方向的胶结破坏，宏观上表现为界面的滑移破坏。

4 结 论

本文采用颗粒离散单元法模拟了岩石砂浆锚杆的拉拔试验，分析了荷载-位移关系、轴力分布和界面剪应力分布，并研究了锚固段界面的微观破坏机制。主要研究结论为：

(1)模拟的荷载-位移曲线与普通岩石锚杆试验结果基本一致，拉拔峰值荷载随锚固长度的增加而增大，但增幅逐渐减小，界面平均粘结强度随着锚固长度的增加而减小。

(2)达到峰值荷载之后，轴力和界面剪应力沿杆长方向出现增减不同步的现象，即靠近拉拔端的轴力和剪应力首先下降，其他部位的轴力和剪应力依次出现先增后减的变化，表明锚固段界面产生渐进破坏。

(3)试样微观胶结破坏主要为拉伸破坏，锚固段破坏在宏观上表现为界面滑移，在微观上主要表现为锚杆受力方向上的微裂纹扩展。

参考文献：

[1] 徐祯祥，闫莫明，苏自约.岩土锚固技术与西部开发[M].北京:人民交通出版社, 2002.

[2] 尤春安,战玉宝.预应力锚索锚固段界面滑移的细观力学分析[J].岩石力学与工程学报,2009,28(10):1976-1985.

[3] HOBST L, ZAJIC J. Anchoring in Rock and Soil [M]. New York: Elsevier Scientific Publishing Company, 1983.

[4] ELIGEHAUSEN R, LEHR B,MESZAROS J,etal. Behavior and Design of Anchorage with Bonded Anchors under Tension Load[C]∥Anchoring & Grouting—Proceedings of International Conference on Anchoring & Grouting Towards the New Century, Guangzhou, October 9, 1999. Guangzhou: Zhongshan University Publishing House, 1999:103-115.

[5] 张永兴,饶枭宇,唐树名，等.压花锚锚固性能的试验研究与数值分析[J].岩石力学与工程学报,2008,27(3):607-614.

[6] FARMER I W. Stress Distribution Along a Resin Grouted Rock Anchor[J]. International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts, 1975, 12(11):347-351.

[7] 杨松林,荣 冠,朱焕春.混凝土中锚杆荷载传递机理的理论分析和现场实验[J].岩土力学，2001,22(1):71-74.

[8] 陈妙峰, 唐德高，周早生，等.锚杆锚固机理试验研究[J].建筑技术开发,2003,30(4):21-23.

[9] 张发明,陈祖煜,刘 宁.岩体与锚固体间粘结强度的确定[J]. 岩土力学,2001,22(4):470-473.

[10] 朱焕春,荣 冠,肖 明,等.张拉荷载下全长黏结锚杆工作机理试验研究[J].岩石力学与工程学报,2002,21(3):379-384.

[11] ZHANG B R, BENMOKRANE B, CHENNOUF A. Prediction of Tensile Capacity of Bond Anchorages for FRP Tendons[J]. Journal of Composites for Construction, 2000, 4(2):39-47.

[12] 江文武,徐国元,马长年.FLAC-3D的锚杆拉拔数值模拟试验[J].哈尔滨工业大学学报,2009,41(10):129-133.

[13] 张 凯,杨 庆,蒋景彩,等.全长粘结岩石锚杆拉拔数值模拟[J].大连理工大学学报,2013,53(5):710-714.

[14] POTYONDY D O, CUNDALL P A. A Bonded-particle Model for Rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(8):1329-1364.

[15] 赵同彬,尹延春,谭云亮，等.锚杆界面力学试验及剪应力传递规律细观模拟分析[J].采矿与安全工程学报,2011,28(2): 220-224.

[16] 张思峰.预应力内锚固段作用机理及其耐久性研究[D].上海:同济大学,2007.

[17] 蒋明镜,方 威,司马军.模拟岩石的平行粘结模型微观参数标定[J].山东大学学报(工学版),2015,45(4):50-56.

[18] 蒋明镜,陈 贺,刘 芳.岩石微观胶结模型及离散元数值仿真方法初探[J].岩石力学与工程学报,2013,32(1):15-23.

[19] JIANG M J, KONRAD J M, LEROUEIL S. An Efficient Technique for Generating Homogeneous Specimens for DEM Studies[J]. Computers and Geotechnics, 2003, 30(7): 579-597.

[20] JIANG M J, CHEN H, CROSTA G B. Numerical Modeling of Rock Mechanical Behavior and Fracture Propagation by a New Bond Contact Model[J]. International Journal of Rock Mechanics and Mining Sciences, 2015, 78:175-189.

[21] 陈 贺.岩石宏微观力学特性及高陡研制边坡的离散元数值模拟[D].上海:同济大学,2013.

[22] MARTIN C D. The Strength of Massive Lac du Bonnet Granite Around Underground Openings[D]. Winnipeg, Canada: University of Manitoba, 1993.

[23] 韩 军,陈 强,刘元坤，等.锚杆灌浆体与岩(土)体间的粘结强度[J].岩石力学与工程学报,2005,24(19):3482-3486.

[24] 尤春安,高 明,张利民，等.锚固体应力分布的试验研究[J]. 岩土力学,2004,25(增刊):63-66.