DDA方法中節(jié)理參數(shù)定義方式的改進

王金梅, 張迎賓, 趙興權(quán),2, 余鵬程, 王 潘

(1. 西南交通大學(xué)土木工程學(xué)院, 四川 成都 610031; 2. 山東建筑大學(xué)土木工程學(xué)院, 山東 濟南 250101; 3. 四川省地質(zhì)工程勘察院, 四川 成都 610072)

DDA (discontinuous deformation analysis)方法[1]是以研究非連續(xù)塊體系統(tǒng)不連續(xù)位移和變形為目的的一種數(shù)值方法,能夠很好地模擬塊體間的滑動、張開和閉合,自提出以來,已廣泛地應(yīng)用于滑坡[2-6]、落石[7-9]、隧洞坍塌[10-11]、開挖爆破[12-13]和傾倒破壞[14]等許多工程領(lǐng)域.

許多學(xué)者對DDA進行了更加深入的研究.作為一種非連續(xù)介質(zhì)力學(xué)數(shù)值分析方法,DDA方法的最大特點是可以反映巖石塊體之間接觸面的滑移、分離與傾覆旋轉(zhuǎn)等大變形大位移問題,因而,塊體間的接觸相關(guān)問題成為DDA方法的核心問題之一.針對DDA方法在接觸上存在的一些缺陷,許多學(xué)者做了進一步的修正和改進.焦玉勇等[15]在DDA方法中加入塊體自適應(yīng)剖分算法,并引入虛實節(jié)理轉(zhuǎn)換來模擬裂紋的擴展、貫通到碎裂、坍塌的全過程.夏才初等[16]利用表征虛節(jié)理向?qū)嵐?jié)理轉(zhuǎn)化程度的虛節(jié)理連通率對虛節(jié)理力學(xué)性質(zhì)參數(shù)的弱化規(guī)律進行研究,基于Jennings強度準則和最大抗拉強度準則導(dǎo)出了節(jié)理力學(xué)性質(zhì)參數(shù)弱化函數(shù),以模擬斷續(xù)節(jié)理的開裂.初始的DDA程序由于其內(nèi)在的接觸判斷、接觸處理方式導(dǎo)致材料粘聚力不能準確地被考慮,使DDA的穩(wěn)定性分析功能被大打折扣,對于此問題,Wang等[17]采用接觸界面相對位移量作為界面粘聚力的失效判定,提出了與接觸界面相對位移相關(guān)的剪切強度準則;Zhang等[18]針對邊邊接觸提出了一種用邊邊接觸上的兩個角邊接觸共同判定整個邊邊接觸上的接觸狀態(tài)的處理方式,以提高剪切作用下考慮粘聚力時節(jié)理上的計算精度;馬永政等[19]在DDA中引入剛塑性模型,補償接觸點在滑面上的黏聚力作用,即令粘聚力在整個滑移過程中均保持不變,使得用DDA方法分析邊坡的安全系數(shù)明顯提高.Jiao等[20]在接觸力的法向方向引入分段力-位移關(guān)系本構(gòu),切向采用摩爾庫倫準則,以模擬巖石開裂的整個過程.為了防止應(yīng)力波在模型邊界上的反射,Jiao[21]、張秀麗[22]等在DDA程序中加入了對黏性邊界的模擬,研究了節(jié)理面對應(yīng)力波傳播的影響.Bao等[23]針對DDA中用最短路徑法搜索角角接觸時接觸邊不唯一的缺陷,提出了兩種改進的最短路徑搜索方法:通過施加一個臨時的角角接觸彈簧去決定接觸塊體的運動方向;當接觸角發(fā)生嵌入時,用接觸角在上一時步的軌跡來決定進入邊.Fan等[24]對角角接觸的預(yù)處理提出了兩個臨時彈簧法,對角角接觸向角邊接觸的轉(zhuǎn)換提出了方向角法.

目前,學(xué)者們對DDA方法在其接觸上的改進工作主要集中在接觸的搜索、虛實節(jié)理的轉(zhuǎn)換以及節(jié)理強度的處理等方面,對DDA方法及其程序計算過程中其節(jié)理上的參數(shù)賦值卻少有提及.在原DDA方法中,節(jié)理參數(shù)的定義是根據(jù)分析初始時刻塊體系統(tǒng)中塊體的分布和塊體間的節(jié)理性質(zhì)將其定義到塊體邊上的,計算時按接觸算法選擇某一接觸邊上的參數(shù)進行計算.這種節(jié)理強度由接觸對中單邊決定的節(jié)理參數(shù)定義,對于含有多種節(jié)理的塊體系統(tǒng),在計算過程中任意塊體間發(fā)生接觸形成節(jié)理時的節(jié)理參數(shù)選取是不適用的,會導(dǎo)致賦有不同節(jié)理參數(shù)的兩條塊體邊發(fā)生接觸時節(jié)理參數(shù)選取的非一致性.對此,本文提出一種由節(jié)理兩側(cè)接觸邊共同決定的節(jié)理參數(shù)定義方式以期對其進行改進.

1 DDA方法概述

1.1 DDA基本原理

DDA將研究對象看成不同塊體組成的塊體系統(tǒng),塊體之間的相互作用通過接觸彈簧來實現(xiàn),荷載按增量方法分時步施加到塊體,每一時步塊體的變形和位移滿足小變形和小位移假定,塊體的大變形和大位移是由眾多時步的小變形和小位移疊加的結(jié)果.設(shè)每一塊體通體有常應(yīng)力和常應(yīng)變,每一時步內(nèi),塊體i內(nèi)部任意一點(x,y)的總位移可用6個位移不變量表示,即

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

(1)

式中：(u0,v0)是塊體內(nèi)特殊點(x0,y0)的剛體位移;r0是塊體繞轉(zhuǎn)動中心(x0,y0)的轉(zhuǎn)動角;(εx,εy,γxy)是該塊體的正應(yīng)變和剪應(yīng)變.

塊體內(nèi)任意點(x,y)的位移(u，v)等于平動、轉(zhuǎn)動和變形共同作用的總和,可表示為

(2)

式中:Ti為塊體i的位移轉(zhuǎn)換矩陣，

各個塊體是連接的,并靠塊體間的接觸和對單個塊體的位移約束形成一個塊體系統(tǒng).通過勢能最小原理求得整個塊體系統(tǒng)的基本方程為

(3)

式中:

k為系統(tǒng)塊體總數(shù),因每個塊體有6個自由度(u0,v0,r0,εx,εy,γxy),式(3)給出的系數(shù)矩陣中每個元素Kij(i,j=1,2,…,k)是一個6×6的子矩陣,Di、Fi是6×1子矩陣；

子矩陣Kii與塊體材料特性和幾何尺寸有關(guān)；

Kij(i≠j)由塊體i和塊體j間的接觸條件決定；

Fi是在塊體i上分配給6個變形變量的荷載.

1.2 DDA接觸理論

塊體系統(tǒng)中,任一塊體可能會與其他塊體發(fā)生接觸.DDA方法中,塊體之間存在3種接觸方式:角對角接觸、角對邊接觸、邊對邊接觸.運算中將一對邊邊接觸處理成兩對角邊接觸,相當于塊體之間僅有角對角、角對邊兩種接觸,如圖1所示，圖中，節(jié)點編號即為角的編號，同時兩個角的編號代表一條邊，如，12代表角1、角2及邊12.若采用Mohr-Coulomb強度準則,則涉及3個節(jié)理參數(shù):摩擦角φ、粘聚力c、抗拉強度σt,這些參數(shù)可根據(jù)相應(yīng)的材料試驗得到.

對于角角接觸(點對點的接觸)、單純的角邊接觸(點線接觸),計算中只需用到φ這一個參數(shù).對于由邊邊接觸轉(zhuǎn)化而來的角邊接觸(線接觸),根據(jù)接觸的性質(zhì),計算中除涉及到φ外,還可能用到c、σt參數(shù).由于計算中將邊邊接觸處理成兩個角邊接觸,邊邊接觸上由c提供的抗剪能力和σt提供的抗拉能力則根據(jù)接觸長度均分到兩個等價的角邊接觸上.塊體的運動過程中,不允許塊體之間受拉和嵌入.若塊體在某個接觸上發(fā)生嵌入,則在該接觸的法向和切向施加兩個給定變形狀態(tài)的剛硬彈簧,使得迭代計算后塊體系統(tǒng)在該部位不再發(fā)生嵌入.DDA方法中每一時步塊體的接觸狀態(tài)包括:閉合、滑動、張開3種狀態(tài),塊體接觸狀態(tài)由式(4)決定.

(4)

式中:

Kg、Ks分別為法向和切向彈簧剛度;

dg、ds分別為法向和切向嵌入距離;

l為1/2接觸長度.

若接觸狀態(tài)為閉合,則在相應(yīng)接觸部位法向和切向施加兩個彈簧;若接觸狀態(tài)為滑動,則在接觸部位施加法向彈簧,切向加摩擦力;若接觸狀態(tài)為張開,則不施加彈簧,直至系統(tǒng)中所有的接觸部位都滿足不嵌入和無張拉準則,則本時步計算完畢,進入下一時步.

(a) 角-角接觸(b) 角-邊接觸(c) 邊-邊接觸圖1 DDA塊體接觸類型Fig.1 Contact types of blocks in DDA

2 原DDA中接觸參數(shù)定義存在的問題

巖體中由各種成因形成的有規(guī)律的破裂面稱為節(jié)理,天然巖體是由節(jié)理和巖塊組成.從巖體力學(xué)屬性來看,可認為完整巖體屬連續(xù)介質(zhì)力學(xué)范疇,節(jié)理巖體因受節(jié)理切割的影響,認為其部分屬非連續(xù)介質(zhì)力學(xué)范疇.從巖體的力學(xué)強度來看,巖石的強度與組成巖體的巖塊和節(jié)理力學(xué)性質(zhì)有很大不同,通常,節(jié)理強度低于巖石強度,節(jié)理巖體的強度處在完整巖石強度和節(jié)理強度之間.因此,節(jié)理的存在造成了巖體介質(zhì)的不連續(xù)性,并使巖體與完整巖石的力學(xué)特性之間有很大差異,因而節(jié)理面的性質(zhì)將大大影響整個巖體的穩(wěn)定性.

節(jié)理面亦稱結(jié)構(gòu)面,是由兩個表面組成.節(jié)理強度主要受到節(jié)理面的粗糙度和起伏度、面壁強度、開度與填充物等因素影響,與構(gòu)成節(jié)理的兩個表面息息相關(guān).

用DDA方法進行數(shù)值模擬時,需要根據(jù)節(jié)理分布對系統(tǒng)進行塊體切割,對完整巖石切割后的塊體,塊體之間以虛節(jié)理表示,計算時用完整巖石的強度進行處理,虛節(jié)理破壞后即成為真實節(jié)理,節(jié)理強度隨之降低.原DDA方法程序中在處理節(jié)理強度參數(shù)時默認將節(jié)理參數(shù)賦值到節(jié)理兩側(cè)塊體的塊體邊上,沿塊體外側(cè)將節(jié)理參數(shù)編號保存到塊體邊左側(cè)的節(jié)點中,如圖2所示.塊體i中每條邊上的節(jié)理參數(shù)保存到了與邊左側(cè)顏色一致的節(jié)點信息中(包括塊體的自由邊或臨空邊,程序中也需要對其賦值節(jié)理參數(shù)).計算時,選擇接觸邊(一組邊邊接觸視為兩個角邊接觸)所賦的節(jié)理參數(shù)進行接觸計算,如圖2中的塊體j的角7與塊體i的邊56形成的角邊接觸VC7-56(VC表示角邊接觸,其下標分表為接觸角和接觸邊),則以塊體i的邊56上所賦的節(jié)理參數(shù)作為計算標準.

圖2 節(jié)理參數(shù)對應(yīng)賦值點示意Fig.2 The corresponding vertex of the joint parameters diagram

計算初始時,根據(jù)已知的節(jié)理參數(shù)和巖石參數(shù)對DDA塊體系統(tǒng)中節(jié)理接觸兩側(cè)的塊體邊賦上相同的節(jié)理參數(shù).在塊體系統(tǒng)的運動過程中,隨著初始接觸的破壞、新接觸的生成,當任意兩個塊體由于接觸而形成新的接觸時,若新接觸的兩條接觸邊上所賦的節(jié)理參數(shù)不一致時,其接觸節(jié)理面上的計算參數(shù)的選取則存在多種情況,并給計算帶來極大的誤差.例如:n1、n2為二維塊體系統(tǒng)中的兩個任意塊體，如圖3所示.塊體n1的邊14和塊體n2的邊58 定義的節(jié)理參數(shù)編號分別為J1、J2,當前時步下(第s1步),兩塊體尚未發(fā)生接觸.計算中某時刻(第s2步,s2>s1)假設(shè)塊體n1的邊14與塊體n2的邊58發(fā)生邊邊接觸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接觸上的參數(shù)取J2節(jié)理參數(shù)編號下的參數(shù)值,VC5-14接觸上的參數(shù)取J1節(jié)理參數(shù)編號下的參數(shù)值;圖4(b)中接觸均取J2節(jié)理參數(shù)編號下的參數(shù)值;圖4(c)中接觸均取J1節(jié)理參數(shù)編號下的參數(shù).圖4(a)～(c)中,14、58兩種結(jié)構(gòu)面形成的同一性質(zhì)的節(jié)理應(yīng)該有相同的節(jié)理參數(shù).在原DDA方法模擬中,由于節(jié)理兩側(cè)接觸邊參數(shù)的不同,其在不同情況下接觸計算參數(shù)可能不一致,從而使得計算結(jié)果與實際情況不相符.

(a) (b) (c) 圖3 第s1步塊體接觸狀態(tài)Fig.3 Contact state of blocks in the s1th step

(a) (b) (c) 圖4 第s2步塊體狀態(tài)(邊邊接觸)Fig.4 Contact state of blocks in the s2 th step(edge to edge contact)

3 接觸參數(shù)定義的改進

原DDA方法中將節(jié)理強度參數(shù)賦值到塊體的每一條邊上,計算時只選取接觸的某一邊上的強度參數(shù)進行計算,即節(jié)理強度由單邊決定.當某一對邊邊接觸的兩個接觸邊上的強度參數(shù)不一致時,同一性質(zhì)的節(jié)理在不同情況下接觸計算參數(shù)就可能不一致.造成這種錯誤的主要原因是原DDA方法中節(jié)理強度由單邊參數(shù)決定的節(jié)理參數(shù)賦值方式與實際物理情況是不相符的,實際物理情況下,節(jié)理強度是由組成節(jié)理的兩個結(jié)構(gòu)面共同決定的.對此,本文提出一種改進方法,即在一個分析初始時具有m種實節(jié)理的塊體系統(tǒng)中,本文額外定義(m2-m)/2種節(jié)理參數(shù),來處理計算過程中當某一對新形成的邊邊接觸的兩個接觸邊上的強度參數(shù)不一致時節(jié)理強度參數(shù)取值的問題.在石根華老師編寫的原DDA程序[1]基礎(chǔ)上改進后的程序計算流程圖如圖5所示.

圖5 改進后的DDA計算流程Fig.5 Flow chart of the improved DDA

圖5中：Jvi為虛節(jié)理的節(jié)理參數(shù)編號；φvi、cvi、σvti為虛節(jié)理相應(yīng)的節(jié)理力學(xué)參數(shù)；cn為當前計算時步內(nèi)搜索到的接觸數(shù)；ng為當前計算時步.

設(shè)節(jié)理參數(shù)編號為Ja、Jb(a,b=1,2,..,m,其相應(yīng)的節(jié)理參數(shù)分別為(φa、ca、σta)、(φb、cb、σtb))的兩條接觸邊形成的節(jié)理，其節(jié)理參數(shù)編號定義為Jab或Jba(相應(yīng)的節(jié)理參數(shù)為(φab、cab、σtab),Jab=Jba),節(jié)理參數(shù)編號均為Ja的兩條接觸邊形成的節(jié)理，其節(jié)理參數(shù)編號仍為Ja,如表1所示.(m2-m)/2種節(jié)理參數(shù)可通過試驗獲得.通過節(jié)理參數(shù)編號即可取得該節(jié)理的參數(shù),若采用摩爾庫倫強度準則和最大抗拉強度準則,則節(jié)理強度參數(shù)包括:摩擦角φab、粘聚力cab、抗拉強度σtab.

表1 節(jié)理接觸參數(shù)定義Tab.1 The definition of joint parameters

4 算例驗證

以塊體的滑動為例,如圖6所示.塊體n2、n3固定,初始時,塊體n1在塊體n3上,并以一定速度v從塊體n3上滑動到塊體n2上.圖中,1、2、5、6、9、10為塊體的部分節(jié)點編號.在原DDA方法中,由于初始時塊體n1只與塊體n3接觸,可預(yù)先給定塊體n1、n3間的節(jié)理強度參數(shù),并賦值到塊體n1、n3接觸邊相應(yīng)的節(jié)點上;而塊體n1、n2尚未發(fā)生接觸.計算過程中,若其發(fā)生接觸,則其接觸上的節(jié)理強度參數(shù)只能根據(jù)其接觸面上所賦的參數(shù)按接觸算法進行給定.設(shè)塊體n1剛好完全滑到塊體n2上時的速度為v0(v0=5 m/s),如圖7所示.

圖7中分別給出了計算的兩種工況:工況1中,塊體n1、n2的材料編號分別為R1、R2,邊12、56上的節(jié)理參數(shù)編號分別為J1、J2;工況2中,塊體n1、n2的材料編號分別為R2、R1,邊12、56上的節(jié)理參數(shù)編號分別為J2、J1;其力學(xué)參數(shù)和界面力學(xué)參數(shù)如表2～3所示.兩種工況下,當塊體n1剛好完全滑倒塊體n2上后,塊體n1的速度v、位移隨時間的變化曲線如圖8所示.

兩種工況下,形成塊體n1、n2間的接觸節(jié)理的兩個結(jié)構(gòu)面性質(zhì)相同,應(yīng)該有相同的節(jié)理參數(shù),塊體n1應(yīng)該有相同的速度和位移.然而,從圖8可見,兩種情況下,塊體n1以相同的速度從塊體n2上開始滑動,經(jīng)過同一性質(zhì)的節(jié)理,其速度、位移并不相同.這是因為在原DDA方法中,當邊12與邊56發(fā)生邊邊接觸時,程序自動將其處理為兩個角邊接觸VC1-65、VC2-65,計算中選擇角邊接觸中接觸邊上的參數(shù)進行計算(這里兩個角邊接觸的接觸邊均為56),因而當兩個塊體交換材料時,原DDA方法均選擇接觸邊56上的參數(shù)進行計算,這就導(dǎo)致同一性質(zhì)的節(jié)理間的接觸節(jié)理參數(shù)不唯一,從而引起計算結(jié)果的不一致.

圖6 塊體滑動算例Fig.6 Example for block sliding

(a) 工況1

(b) 工況2圖7 兩種塊體滑動Fig.7 Two working condition of the block sliding

表2 塊體力學(xué)參數(shù)Tab.2 The mechanical parameters of blocks

表3 界面力學(xué)參數(shù)Tab.3 The mechanical parameters of joints

實際上節(jié)理性質(zhì)由其兩個表面共同決定,即使塊體1、2交換接觸邊12、56上的節(jié)理參數(shù)和接觸EC12-56仍為同一性質(zhì)的節(jié)理,其節(jié)理參數(shù)應(yīng)該不變.根據(jù)本文提出的節(jié)理參數(shù)定義方法,設(shè)原DDA方法中由節(jié)理參數(shù)編號為J1、J2的兩條接觸邊構(gòu)成的接觸其節(jié)理參數(shù)編號為J12,由本文改進的DDA程序算得兩種工況下塊體n1的速度、位移變化如圖9所示.

由圖9可見,塊體n1的速度、位移變化一致,說明了本文的節(jié)理參數(shù)定義方式的有效性,該定義方式能在計算中反映更加貼近實際的節(jié)理參數(shù).

(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 結(jié)束語

塊體間的接觸處理相關(guān)問題是DDA方法計算中的核心問題,在原DDA方法中,程序在處理節(jié)理強度參數(shù)時默認將節(jié)理參數(shù)賦值到節(jié)理兩側(cè)塊體的塊體邊上,計算時選取某一側(cè)接觸邊上的節(jié)理參數(shù)進行計算.這種處理方式,對于塊體系統(tǒng)運動過程中不同參數(shù)性質(zhì)的塊體間接觸時接觸節(jié)理參數(shù)的選取是不適宜的,將造成塊體間同一節(jié)理上的節(jié)理參數(shù)取值的不一致.

針對上述問題,本文從節(jié)理的實際物理意義出發(fā),提出了一種由節(jié)理兩側(cè)接觸邊共同決定的節(jié)理參數(shù)定義方式,編寫了相應(yīng)的程序整合到原DDA程序中,并給出了具體算例.本文的定義方式更符合物理實際的節(jié)理參數(shù)定義,可以在計算中反映更加貼近實際的節(jié)理參數(shù);克服了原DDA方法中節(jié)理參數(shù)選取時出現(xiàn)的非一致情況;程序簡潔、思路清晰,可以很好的整合到原程序中去,能夠方便地實現(xiàn)了分析初期和分析過程中塊體間節(jié)理參數(shù)的賦值.

[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數(shù)值方法在巖質(zhì)邊坡傾倒破壞分析中的應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報, 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] 焦玉勇,張秀麗,劉泉聲,等. 用非連續(xù)變形分析方法模擬巖石裂紋擴展[J]. 巖石力學(xué)與工程學(xué)報,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] 夏才初,許崇幫. 非連續(xù)變形分析(DDA)中斷續(xù)節(jié)理擴展的模擬方法研究和試驗驗證[J]. 巖石力學(xué)與工程學(xué)報,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法計算邊坡安全系數(shù)的黏聚力影響分析[J]. 巖土工程學(xué)報,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] 張秀麗,焦玉勇,劉泉聲,等. 節(jié)理對爆炸波傳播影響的數(shù)值研究[J]. 巖土力學(xué),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.