DDA方法中节理参数定义方式的改进

王金梅, 张迎宾, 赵兴权,2, 余鹏程, 王 潘

(1. 西南交通大学土木工程学院, 四川 成都 610031; 2. 山东建筑大学土木工程学院, 山东 济南 250101; 3. 四川省地质工程勘察院, 四川 成都 610072)

DDA (discontinuous deformation analysis)方法[1]是以研究非连续块体系统不连续位移和变形为目的的一种数值方法,能够很好地模拟块体间的滑动、张开和闭合,自提出以来,已广泛地应用于滑坡[2-6]、落石[7-9]、隧洞坍塌[10-11]、开挖爆破[12-13]和倾倒破坏[14]等许多工程领域.

许多学者对DDA进行了更加深入的研究.作为一种非连续介质力学数值分析方法,DDA方法的最大特点是可以反映岩石块体之间接触面的滑移、分离与倾覆旋转等大变形大位移问题,因而,块体间的接触相关问题成为DDA方法的核心问题之一.针对DDA方法在接触上存在的一些缺陷,许多学者做了进一步的修正和改进.焦玉勇等[15]在DDA方法中加入块体自适应剖分算法,并引入虚实节理转换来模拟裂纹的扩展、贯通到碎裂、坍塌的全过程.夏才初等[16]利用表征虚节理向实节理转化程度的虚节理连通率对虚节理力学性质参数的弱化规律进行研究,基于Jennings强度准则和最大抗拉强度准则导出了节理力学性质参数弱化函数,以模拟断续节理的开裂.初始的DDA程序由于其内在的接触判断、接触处理方式导致材料粘聚力不能准确地被考虑,使DDA的稳定性分析功能被大打折扣,对于此问题,Wang等[17]采用接触界面相对位移量作为界面粘聚力的失效判定,提出了与接触界面相对位移相关的剪切强度准则;Zhang等[18]针对边边接触提出了一种用边边接触上的两个角边接触共同判定整个边边接触上的接触状态的处理方式,以提高剪切作用下考虑粘聚力时节理上的计算精度;马永政等[19]在DDA中引入刚塑性模型,补偿接触点在滑面上的黏聚力作用,即令粘聚力在整个滑移过程中均保持不变,使得用DDA方法分析边坡的安全系数明显提高.Jiao等[20]在接触力的法向方向引入分段力-位移关系本构,切向采用摩尔库伦准则,以模拟岩石开裂的整个过程.为了防止应力波在模型边界上的反射,Jiao[21]、张秀丽[22]等在DDA程序中加入了对黏性边界的模拟,研究了节理面对应力波传播的影响.Bao等[23]针对DDA中用最短路径法搜索角角接触时接触边不唯一的缺陷,提出了两种改进的最短路径搜索方法:通过施加一个临时的角角接触弹簧去决定接触块体的运动方向;当接触角发生嵌入时,用接触角在上一时步的轨迹来决定进入边.Fan等[24]对角角接触的预处理提出了两个临时弹簧法,对角角接触向角边接触的转换提出了方向角法.

目前,学者们对DDA方法在其接触上的改进工作主要集中在接触的搜索、虚实节理的转换以及节理强度的处理等方面,对DDA方法及其程序计算过程中其节理上的参数赋值却少有提及.在原DDA方法中,节理参数的定义是根据分析初始时刻块体系统中块体的分布和块体间的节理性质将其定义到块体边上的,计算时按接触算法选择某一接触边上的参数进行计算.这种节理强度由接触对中单边决定的节理参数定义,对于含有多种节理的块体系统,在计算过程中任意块体间发生接触形成节理时的节理参数选取是不适用的,会导致赋有不同节理参数的两条块体边发生接触时节理参数选取的非一致性.对此,本文提出一种由节理两侧接触边共同决定的节理参数定义方式以期对其进行改进.

1 DDA方法概述

1.1 DDA基本原理

DDA将研究对象看成不同块体组成的块体系统,块体之间的相互作用通过接触弹簧来实现,荷载按增量方法分时步施加到块体,每一时步块体的变形和位移满足小变形和小位移假定,块体的大变形和大位移是由众多时步的小变形和小位移叠加的结果.设每一块体通体有常应力和常应变,每一时步内,块体i内部任意一点(x,y)的总位移可用6个位移不变量表示,即

Di=(u0,v0,r0,εx,εy,γxy)T,

(1)

式中：(u0,v0)是块体内特殊点(x0,y0)的刚体位移;r0是块体绕转动中心(x0,y0)的转动角;(εx,εy,γxy)是该块体的正应变和剪应变.

块体内任意点(x,y)的位移(u，v)等于平动、转动和变形共同作用的总和,可表示为

(2)

式中:Ti为块体i的位移转换矩阵，

各个块体是连接的,并靠块体间的接触和对单个块体的位移约束形成一个块体系统.通过势能最小原理求得整个块体系统的基本方程为

(3)

式中:

k为系统块体总数,因每个块体有6个自由度(u0,v0,r0,εx,εy,γxy),式(3)给出的系数矩阵中每个元素Kij(i,j=1,2,…,k)是一个6×6的子矩阵,Di、Fi是6×1子矩阵；

子矩阵Kii与块体材料特性和几何尺寸有关；

Kij(i≠j)由块体i和块体j间的接触条件决定；

Fi是在块体i上分配给6个变形变量的荷载.

1.2 DDA接触理论

块体系统中,任一块体可能会与其他块体发生接触.DDA方法中,块体之间存在3种接触方式:角对角接触、角对边接触、边对边接触.运算中将一对边边接触处理成两对角边接触,相当于块体之间仅有角对角、角对边两种接触,如图1所示，图中，节点编号即为角的编号，同时两个角的编号代表一条边，如，12代表角1、角2及边12.若采用Mohr-Coulomb强度准则,则涉及3个节理参数:摩擦角φ、粘聚力c、抗拉强度σt,这些参数可根据相应的材料试验得到.

对于角角接触(点对点的接触)、单纯的角边接触(点线接触),计算中只需用到φ这一个参数.对于由边边接触转化而来的角边接触(线接触),根据接触的性质,计算中除涉及到φ外,还可能用到c、σt参数.由于计算中将边边接触处理成两个角边接触,边边接触上由c提供的抗剪能力和σt提供的抗拉能力则根据接触长度均分到两个等价的角边接触上.块体的运动过程中,不允许块体之间受拉和嵌入.若块体在某个接触上发生嵌入,则在该接触的法向和切向施加两个给定变形状态的刚硬弹簧,使得迭代计算后块体系统在该部位不再发生嵌入.DDA方法中每一时步块体的接触状态包括:闭合、滑动、张开3种状态,块体接触状态由式(4)决定.

(4)

式中:

Kg、Ks分别为法向和切向弹簧刚度;

dg、ds分别为法向和切向嵌入距离;

l为1/2接触长度.

若接触状态为闭合,则在相应接触部位法向和切向施加两个弹簧;若接触状态为滑动,则在接触部位施加法向弹簧,切向加摩擦力;若接触状态为张开,则不施加弹簧,直至系统中所有的接触部位都满足不嵌入和无张拉准则,则本时步计算完毕,进入下一时步.

(a) 角-角接触(b) 角-边接触(c) 边-边接触图1 DDA块体接触类型Fig.1 Contact types of blocks in DDA

2 原DDA中接触参数定义存在的问题

岩体中由各种成因形成的有规律的破裂面称为节理,天然岩体是由节理和岩块组成.从岩体力学属性来看,可认为完整岩体属连续介质力学范畴,节理岩体因受节理切割的影响,认为其部分属非连续介质力学范畴.从岩体的力学强度来看,岩石的强度与组成岩体的岩块和节理力学性质有很大不同,通常,节理强度低于岩石强度,节理岩体的强度处在完整岩石强度和节理强度之间.因此,节理的存在造成了岩体介质的不连续性,并使岩体与完整岩石的力学特性之间有很大差异,因而节理面的性质将大大影响整个岩体的稳定性.

节理面亦称结构面,是由两个表面组成.节理强度主要受到节理面的粗糙度和起伏度、面壁强度、开度与填充物等因素影响,与构成节理的两个表面息息相关.

用DDA方法进行数值模拟时,需要根据节理分布对系统进行块体切割,对完整岩石切割后的块体,块体之间以虚节理表示,计算时用完整岩石的强度进行处理,虚节理破坏后即成为真实节理,节理强度随之降低.原DDA方法程序中在处理节理强度参数时默认将节理参数赋值到节理两侧块体的块体边上,沿块体外侧将节理参数编号保存到块体边左侧的节点中,如图2所示.块体i中每条边上的节理参数保存到了与边左侧颜色一致的节点信息中(包括块体的自由边或临空边,程序中也需要对其赋值节理参数).计算时,选择接触边(一组边边接触视为两个角边接触)所赋的节理参数进行接触计算,如图2中的块体j的角7与块体i的边56形成的角边接触VC7-56(VC表示角边接触,其下标分表为接触角和接触边),则以块体i的边56上所赋的节理参数作为计算标准.

图2 节理参数对应赋值点示意Fig.2 The corresponding vertex of the joint parameters diagram

计算初始时,根据已知的节理参数和岩石参数对DDA块体系统中节理接触两侧的块体边赋上相同的节理参数.在块体系统的运动过程中,随着初始接触的破坏、新接触的生成,当任意两个块体由于接触而形成新的接触时,若新接触的两条接触边上所赋的节理参数不一致时,其接触节理面上的计算参数的选取则存在多种情况,并给计算带来极大的误差.例如:n1、n2为二维块体系统中的两个任意块体，如图3所示.块体n1的边14和块体n2的边58 定义的节理参数编号分别为J1、J2,当前时步下(第s1步),两块体尚未发生接触.计算中某时刻(第s2步,s2>s1)假设块体n1的边14与块体n2的边58发生边边接触EC14-58,(EC表示边边接触,其下标分别为接触的两条边),如图4所示.由于DDA方法中一个边边接触等价于两个角边接触,接触EC14-58将被分别处理成图4(a)VC1-58和VC5-14、图4(b)VC8-14和VC5-14、图4(c)VC1-58和VC4-58两个角边接触.计算时,图4(a)中VC1-58接触上的参数取J2节理参数编号下的参数值,VC5-14接触上的参数取J1节理参数编号下的参数值;图4(b)中接触均取J2节理参数编号下的参数值;图4(c)中接触均取J1节理参数编号下的参数.图4(a)～(c)中,14、58两种结构面形成的同一性质的节理应该有相同的节理参数.在原DDA方法模拟中,由于节理两侧接触边参数的不同,其在不同情况下接触计算参数可能不一致,从而使得计算结果与实际情况不相符.

(a) (b) (c) 图3 第s1步块体接触状态Fig.3 Contact state of blocks in the s1th step

(a) (b) (c) 图4 第s2步块体状态(边边接触)Fig.4 Contact state of blocks in the s2 th step(edge to edge contact)

3 接触参数定义的改进

原DDA方法中将节理强度参数赋值到块体的每一条边上,计算时只选取接触的某一边上的强度参数进行计算,即节理强度由单边决定.当某一对边边接触的两个接触边上的强度参数不一致时,同一性质的节理在不同情况下接触计算参数就可能不一致.造成这种错误的主要原因是原DDA方法中节理强度由单边参数决定的节理参数赋值方式与实际物理情况是不相符的,实际物理情况下,节理强度是由组成节理的两个结构面共同决定的.对此,本文提出一种改进方法,即在一个分析初始时具有m种实节理的块体系统中,本文额外定义(m2-m)/2种节理参数,来处理计算过程中当某一对新形成的边边接触的两个接触边上的强度参数不一致时节理强度参数取值的问题.在石根华老师编写的原DDA程序[1]基础上改进后的程序计算流程图如图5所示.

图5 改进后的DDA计算流程Fig.5 Flow chart of the improved DDA

图5中：Jvi为虚节理的节理参数编号；φvi、cvi、σvti为虚节理相应的节理力学参数；cn为当前计算时步内搜索到的接触数；ng为当前计算时步.

设节理参数编号为Ja、Jb(a,b=1,2,..,m,其相应的节理参数分别为(φa、ca、σta)、(φb、cb、σtb))的两条接触边形成的节理，其节理参数编号定义为Jab或Jba(相应的节理参数为(φab、cab、σtab),Jab=Jba),节理参数编号均为Ja的两条接触边形成的节理，其节理参数编号仍为Ja,如表1所示.(m2-m)/2种节理参数可通过试验获得.通过节理参数编号即可取得该节理的参数,若采用摩尔库伦强度准则和最大抗拉强度准则,则节理强度参数包括:摩擦角φab、粘聚力cab、抗拉强度σtab.

表1 节理接触参数定义Tab.1 The definition of joint parameters

4 算例验证

以块体的滑动为例,如图6所示.块体n2、n3固定,初始时,块体n1在块体n3上,并以一定速度v从块体n3上滑动到块体n2上.图中,1、2、5、6、9、10为块体的部分节点编号.在原DDA方法中,由于初始时块体n1只与块体n3接触,可预先给定块体n1、n3间的节理强度参数,并赋值到块体n1、n3接触边相应的节点上;而块体n1、n2尚未发生接触.计算过程中,若其发生接触,则其接触上的节理强度参数只能根据其接触面上所赋的参数按接触算法进行给定.设块体n1刚好完全滑到块体n2上时的速度为v0(v0=5 m/s),如图7所示.

图7中分别给出了计算的两种工况:工况1中,块体n1、n2的材料编号分别为R1、R2,边12、56上的节理参数编号分别为J1、J2;工况2中,块体n1、n2的材料编号分别为R2、R1,边12、56上的节理参数编号分别为J2、J1;其力学参数和界面力学参数如表2～3所示.两种工况下,当块体n1刚好完全滑倒块体n2上后,块体n1的速度v、位移随时间的变化曲线如图8所示.

两种工况下,形成块体n1、n2间的接触节理的两个结构面性质相同,应该有相同的节理参数,块体n1应该有相同的速度和位移.然而,从图8可见,两种情况下,块体n1以相同的速度从块体n2上开始滑动,经过同一性质的节理,其速度、位移并不相同.这是因为在原DDA方法中,当边12与边56发生边边接触时,程序自动将其处理为两个角边接触VC1-65、VC2-65,计算中选择角边接触中接触边上的参数进行计算(这里两个角边接触的接触边均为56),因而当两个块体交换材料时,原DDA方法均选择接触边56上的参数进行计算,这就导致同一性质的节理间的接触节理参数不唯一,从而引起计算结果的不一致.

图6 块体滑动算例Fig.6 Example for block sliding

(a) 工况1

(b) 工况2图7 两种块体滑动Fig.7 Two working condition of the block sliding

表2 块体力学参数Tab.2 The mechanical parameters of blocks

表3 界面力学参数Tab.3 The mechanical parameters of joints

实际上节理性质由其两个表面共同决定,即使块体1、2交换接触边12、56上的节理参数和接触EC12-56仍为同一性质的节理,其节理参数应该不变.根据本文提出的节理参数定义方法,设原DDA方法中由节理参数编号为J1、J2的两条接触边构成的接触其节理参数编号为J12,由本文改进的DDA程序算得两种工况下块体n1的速度、位移变化如图9所示.

由图9可见,块体n1的速度、位移变化一致,说明了本文的节理参数定义方式的有效性,该定义方式能在计算中反映更加贴近实际的节理参数.

(a) 位移(b) 速度图8 原DDA计算的块体n1速度、位移曲线Fig.8 The velocity and displacement time history curves of block n1 by original DDA

(a) 位移(b) 速度图9 改进DDA计算的块体n1速度、位移曲线Fig.9 The velocity and displacement time history curves of block n1 by improved DDA

5 结束语

块体间的接触处理相关问题是DDA方法计算中的核心问题,在原DDA方法中,程序在处理节理强度参数时默认将节理参数赋值到节理两侧块体的块体边上,计算时选取某一侧接触边上的节理参数进行计算.这种处理方式,对于块体系统运动过程中不同参数性质的块体间接触时接触节理参数的选取是不适宜的,将造成块体间同一节理上的节理参数取值的不一致.

针对上述问题,本文从节理的实际物理意义出发,提出了一种由节理两侧接触边共同决定的节理参数定义方式,编写了相应的程序整合到原DDA程序中,并给出了具体算例.本文的定义方式更符合物理实际的节理参数定义,可以在计算中反映更加贴近实际的节理参数;克服了原DDA方法中节理参数选取时出现的非一致情况;程序简洁、思路清晰,可以很好的整合到原程序中去,能够方便地实现了分析初期和分析过程中块体间节理参数的赋值.

[1] SHI Genhua. Discontinuous deformation analysis: a new numerical model for the statics and dynamics of block system[D]. Berkeley: University of California, 1988.

[2] SITAR N, MACLAUGHLIN M M, DOOLIN D M. Influence of kinematics on landslide mobility and failure mode[J]. Journal of Geotechnical & Geoenvir onmental Engineering, 2005, 131(6): 716-728.

[3] WU Jianhong. Applying discontinuous deformation analysis to assess the constrained area of the unstable Chiu-fen-erh-shan landslide slope[J]. International Journal for Numerical & Analytic Methods in Geomechanics, 2010, 31(5): 649-666.

[4] KVELDSVIK V, EINSTEIN H H, NILSEN B, et al. Numerical analysis of the 650 000 m2Åknes rock slope based on measured displacements and geotechnical data[J]. Rock Mechanics and Rock Engineering, 2009, 42(5): 689-728.

[5] WU Jianhong, CHEN Chunhwa. Application of DDA to simulate characteristics of the Tsaoling landslide[J]. Computers & Geotechnics, 2011, 38(5): 741-750.

[6] ZHANG Yingbin, CHEN Guangqi, ZHENG Lu, et al. Effects of near-fault seismic loadings on run-out of large-scale landslide: a case study[J]. Engineering Geology, 2013, 166(8): 216-236.

[7] CHEN Guangqi. Numerical modelling of rockfall using extended DDA[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(6): 926-931.

[8] MATSUYAMA H, NISHIYAMA S, OHNISHI Y. Practical studies on rockfall simulation by DDA[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2011, 3(1): 57-63.

[9] CHEN Guangqi, ZHENG Lu, ZHANG Yingbin, et al. Numerical simulation in rockfall analysis: a close comparison of 2-D and 3-D DDA[J]. Rock Mechanics and Rock Engineering, 2013, 46(3): 527-41.

[10] WU Jianhong, OHNISHI Y, NISHIYAMA S. Simulation of the mechanical behavior of inclined jointed rock masses during tunnel construction using discontinuous deformation analysis (DDA)[J]. International Journal of Rock Mechanics and Mining Sciences, 2004, 41(5): 731-743.

[11] TSESARSKY M, HATZOR Y. Tunnel roof deflection in blocky rock masses as a function of joint spacing and friction-A parametric study using discontinuous deformation analysis (DDA)[J]. Tunnelling & Underground Space Technology, 2006, 21(1): 29-45.

[12] NING Youjun, YANG Jun, AN Xinmei, et al. Simulation of blast induced crater in jointed rock mass by discontinuous deformation analysis method[J]. Frontiers of Structural and Civil Engineering, 2010, 4(2): 223-232.

[13] NING Youjun, YANG Jun, AN Xinmei, et al. Modelling rock fracturing and blast-induced rock mass failure via advanced discretisation within the discontinuous deformation analysis framework[J]. Computers & Geotechnics, 2011, 38(1): 40-49.

[14] 孙东亚,彭一江,王兴珍. DDA数值方法在岩质边坡倾倒破坏分析中的应用[J]. 岩石力学与工程学报, 2002,121(1): 39-42.

SUN Dongya, PENG Yijiang, WANG Xingzhen. Application of DDA method in stability analysis of topple rock slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 121(1): 3942.

[15] 焦玉勇,张秀丽,刘泉声,等. 用非连续变形分析方法模拟岩石裂纹扩展[J]. 岩石力学与工程学报,2007,26(4): 682-691.

JIAO Yuyong, ZHANG Xiuli, LIU Quansheng, et al. Simulation of rock crack propagation using discontinuous deformation analysis method[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(4): 682-691.

[16] 夏才初,许崇帮. 非连续变形分析(DDA)中断续节理扩展的模拟方法研究和试验验证[J]. 岩石力学与工程学报,2010,29(10): 2027-2033.

XIA Caichu, XU Chongbang. Atudy of fracturing algorithm of intermitteng joint by DDA and experimental validation[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(10): 2027-2033.

[17] WANG Lizhong, JIANG Hongyi, YANG Zhongxuan, et al. Development of discontinuous deformation analysis with displacement-dependent interface shear strength[J]. Computers and Geotechnics. 2013, 47(1): 91-101.

[18] ZHANG Yingbin, XU Qiang, CHEN Guangqi, et al. Extension of discontinuous deformation analysis and application in cohesive-frictional slope analysis[J]. International Journal of Rock Mechanics & Mining Sciences, 2014, 70(9): 533-545.

[19] 马永政,郑宏,朱合华,等. DDA法计算边坡安全系数的黏聚力影响分析[J]. 岩土工程学报,2009,31(7): 1088-1093.

MA Yongzheng, ZHEN Hong, ZHU Hehua, et al. Effect of cohesion on evaluating slope stability factor of safety by DDA method[J]. Chinese Journal of Geotechnical Engineering, 2009, 31(7): 1088-1093.

[20] JIAO Yuyong, ZHANG Xiuli, ZHAO Jian. Two-dimensional DDA contact constitutive model for simulating rock fragmentation[J]. Journal of Engineering Mechanics, 2011, 138(2): 199-209.

[21] JIAO Yuyong, ZHANG Xiuli, ZHAO Jian. Viscous boundary of DDA for modeling stress wave propagation in jointed rock[J]. International Journal of Rock Mechanics and Mining Sciences, 2007, 44(7): 1070-1076.

[22] 张秀丽,焦玉勇,刘泉声,等. 节理对爆炸波传播影响的数值研究[J]. 岩土力学,2008,29(3): 717-721.

ZHANG Xiuli, JIAO Yuyong, LIU Quansheng, et al. Numerical study on effect of joints on blasting wave progation in rock mass[J]. Rock and Soil Mechanics, 2008, 29(3): 717-721.

[23] BAO Huirong, ZHAO Zhiye. An alternative scheme for the corner-corner contact in the two-dimensional Discontinuous Deformation Analysis[J]. Advances in Engineering Software, 2010, 41(2): 206-212.

[24] FAN Huo, HE Siming. An angle-based method dealing with vertex-vertex contact in the two-dimensional discontinuous deformation analysis[J]. Rock Mechanics and Rock Engineering, 2015, 48(5): 2031-2043.