不同采空區(qū)開挖滑坡災(zāi)害與結(jié)構(gòu)體動力響應(yīng)研究

畢鈺璋 孫新坡 何思明 單 雨 王安輝

（1.東南大學(xué)巖土工程研究所，江蘇南京210096；2.四川理工學(xué)院土木工程學(xué)院，四川自貢643000；3.中國科學(xué)院山地災(zāi)害與地表過程重點實驗室，四川成都610041；4.中國科學(xué)院水利部山地災(zāi)害與環(huán)境研究所，四川成都610041；5.中國科學(xué)院青藏高原地球科學(xué)卓越創(chuàng)新中心，北京100101；6.西南交通大學(xué)土木工程學(xué)院，四川成都610041）

滑坡災(zāi)害的發(fā)生通常伴隨著大量勢能轉(zhuǎn)化為動 能的過程，并且以其強沖擊力對周圍的結(jié)構(gòu)體造成極大的破壞［1-4］?；聻?zāi)害的發(fā)生在采礦工程中更為常見，人為的擾動會打破邊坡的力學(xué)平衡，并且造成局部的應(yīng)力集中，進而造成滑坡災(zāi)害的產(chǎn)生［5-8］。尤其是露天開采和地下開采的轉(zhuǎn)換過程中，采空區(qū)會觸發(fā)滑坡或者造成古滑坡的復(fù)活，對邊坡的穩(wěn)定性造成不可估量的影響，并對周圍的交通路線、村莊、各種民用公用結(jié)構(gòu)造成難以估量的威脅［9-10］。

關(guān)于采空區(qū)對邊坡穩(wěn)定性的影響，國內(nèi)外學(xué)者都作了大量的工作。Fujisawa等人［11］建立了相關(guān)的數(shù)學(xué)模型并模擬了二維條件下的采空區(qū)對邊坡穩(wěn)定性的影響；Shao Yong等［6］基于實際工況，采用三維有限元的方法分析了某礦采空區(qū)對周圍山體穩(wěn)定性的影響，并給出了不同開挖條件下的位移云圖；柴保紅等［12］基于強度折減法采用FLAC3D對不同開采條件的邊坡進行了詳細的研究，分析了開挖地點在邊坡下部、邊坡中部、邊坡上部等工況下對邊坡穩(wěn)定性的影響；宋偉東等［13］基于有限差分方法和極限平衡法研究了大冶鐵礦東露天邊坡的工況，分析了露采轉(zhuǎn)地采應(yīng)該滿足的安全條件；丁桂伶［14］基于現(xiàn)場實地的調(diào)查并采用有限元方法對大村澗村采空區(qū)上方的邊坡穩(wěn)定性進行了分析，不僅得到了天然狀態(tài)和飽和狀態(tài)的安全系數(shù)，而且分析了房屋結(jié)構(gòu)的沉降情況。雖然大量的研究都是關(guān)于采空區(qū)和滑坡的，然而前人的研究主要側(cè)重于采空區(qū)對邊坡穩(wěn)定性的影響，對災(zāi)害產(chǎn)生后對周邊結(jié)構(gòu)體的潛在影響并沒有分析，而這個恰恰是防災(zāi)減災(zāi)工程中重要的一個指標(biāo)。

關(guān)于災(zāi)害體和結(jié)構(gòu)體之間的動力機理研究屬于災(zāi)害動力學(xué)的范疇，相關(guān)的研究工作主要集中在研究巖崩、泥石流等災(zāi)害對周邊結(jié)構(gòu)體的影響。Teufelsbauer等［15］基于二維離散元方法研究了顆粒流和結(jié)構(gòu)體之間的相互作用，并得出了邊坡形狀等因素對于顆粒流沖擊力的影響；Zanuttigh等［16］基于淺水方程研究了泥石流和災(zāi)害體的相互作用，揭示了流深等因素對于促進泥石流沖擊力的關(guān)系；畢鈺璋等研究了災(zāi)害和新型結(jié)構(gòu)體的相互的作用［2］，以及約束條件下粗細混合顆粒和結(jié)構(gòu)體之間的動力響應(yīng)［17］，揭示了顆粒分選對于災(zāi)害破壞力的影響?？梢娊沂緸?zāi)害體和結(jié)構(gòu)體之間的動力響應(yīng)規(guī)律對于指導(dǎo)工程應(yīng)用有著十分重要的意義。

本研究采用二維離散元方法，通過和前人工況算例對比來確定相應(yīng)的數(shù)值模擬參數(shù)；進而建立物理模型對不同開挖工況進行模擬，分析其崩塌規(guī)律；最后分析不同采空區(qū)條件下滑坡災(zāi)害與結(jié)構(gòu)體之間的動力響應(yīng)規(guī)律。旨在得到相關(guān)的規(guī)律性結(jié)論，并對實際的采礦工程有一定指導(dǎo)意義。

1 數(shù)值模型及參數(shù)選取

離散單元法（discrete element method）由Cundall和Strack在1979年提出［18-19］。離散單元法的基本組成元素包括二維圓盤單元、三維圓球單元、和墻單元。該方法通過圓盤或者圓球的粘結(jié)進而研究固體材料的力學(xué)特性。本次研究采用PFC2D（二維顆粒流方法）來搭建數(shù)值模型，進而進行一系列的數(shù)值分析。

1.1 顆粒接觸本構(gòu)模型

本次研究中對顆粒間的接觸模型采用平行粘結(jié)模型。平行粘結(jié)模型是顆粒流中常用的一種模型，它是顆粒間相互粘結(jié)的模型，可以通過對顆粒體的膠結(jié)來模擬巖土體［19］。平行粘結(jié)模型的粘結(jié)只能發(fā)生在顆粒接觸點很小的范圍內(nèi)，通常用于模擬粘性摩擦質(zhì)地的材料。從物理的角度來看，通常粘結(jié)巖石和膠結(jié)土的顆粒間的距離是一定的，膠結(jié)體相當(dāng)于“顆粒鍵”的作用，不僅可以傳遞力還可以傳遞彎矩［20］。在本次研究中，DEM主要作為數(shù)值工具來模擬實際工況并對宏觀的滑坡災(zāi)害進行分析，而微觀顆粒間的力學(xué)性能并不是本文分析的重點?；诖?，本研究中“顆粒鍵”模型的選取盡量簡化，從而減少微觀參數(shù)的數(shù)量。

一般而言，平行粘結(jié)模型主要由顆粒密度、顆粒形狀、顆粒尺寸、顆粒分布等微觀特性共同決定。目前關(guān)于顆粒微觀參數(shù)和實際材料的宏觀參數(shù)之間沒有確切的數(shù)學(xué)模型來描述。因此關(guān)于顆粒的微觀參數(shù)的確定一般采用反演的方法，即用顆粒構(gòu)成的宏觀材料和實際材料相對比，使二者宏觀力學(xué)特性相吻合，進而確定用于離散元數(shù)值模擬的微觀力學(xué)參數(shù)。數(shù)值雙軸試驗和剪切試驗通常被用于評價材料的微觀參數(shù)［21-22］。在PFC2D中使用以下關(guān)系計算顆粒剛度以及接觸鍵參數(shù)：

式中，kn是顆粒法向剛度；t是顆粒沿著平面的厚度；Ec是顆粒接觸的楊氏模量；φn和φs分別表示接觸鍵的法向強度和切向強度；σn和σs分別表示材料的法向強度和切向強度；R為2個顆粒的平均半徑。

方程（1）展示了在二維情況下，單個接觸鍵的強度與顆粒半徑成比例，并且當(dāng)其承載的任一應(yīng)力等于或超過材料強度時，接觸鍵將斷裂。

1.2 數(shù)值實驗的參數(shù)確定

離散元模擬中的對象的宏觀力學(xué)行為取決于微觀顆粒的力學(xué)參數(shù)的設(shè)定。目前關(guān)于離散元模擬的參數(shù)選取主要采用參數(shù)反演的方法。有些學(xué)者通過物理試驗和數(shù)值模擬的堆積形態(tài)的擬合來確定最終的參數(shù)［23］；有些學(xué)者通過實際的工程數(shù)據(jù)和數(shù)值模擬的結(jié)果來反演參數(shù)［1-2］。本次研究的目的是采空區(qū)對邊坡動力響應(yīng)的影響，因此如果選取堆積形態(tài)或者力學(xué)參數(shù)的話就不具備很強的說服力。基于此，通過對實際的隧道開挖工程進行模擬，進而研究開挖后產(chǎn)生的邊坡擾動，對比前人的計算結(jié)果來反演本次研究所需的參數(shù)。

Koizumi等［24］在研究隧道開挖對邊坡失穩(wěn)的過程中，采用了不同的計算模型來研究邊坡的最大水平位移和“隧道——軟弱巖層面”距離的關(guān)系。如圖1（a）所示是Koizumi的二維有限元模型，整個邊坡分成3個部分組成：沙土沉積層（1）、軟弱巖層（2）、以及基巖層（3）。其中，隧道距離軟弱巖層面的距離分別表示為：0.5D、1.0D、1.5D、2.0D、3.0D，其中D表示的是隧道半徑。圖1（b）為相同尺寸的離散元數(shù)值模擬的二維模型。

圖2顯示的是邊坡表面最大水平位移和“隧道——軟弱巖層面”距離的關(guān)系曲線圖。Koizumi分別用雙線性模型和應(yīng)變軟化模型分別模擬了不同“隧道—軟弱巖層面”距離值對邊坡表面最大水平位移的影響。中間曲線代表本次離散元數(shù)值模擬的結(jié)果。從圖中可以很好地反映出本次數(shù)值模擬的結(jié)果和Koizumi對工況模擬的結(jié)果不僅在規(guī)律上、數(shù)值上都很接近，進而可以確定本次研究所用的數(shù)值模擬參數(shù)見表1。

?

1.3 本次研究的物理模型

本次研究采用二維圓盤粘結(jié)并模擬滑坡體?；麦w模型的構(gòu)建采用“雨滴法”建模，如圖3（a）所示。相對于滑坡體400 m的長度和75 m的高度，顆粒的半徑設(shè)置成0.6 m和0.8 m 2種，這樣使得顆粒尺寸相對于滑坡體而言很小，有利于更好地去研究邊坡的失穩(wěn)和破壞。而2種以上不同尺寸的顆?？梢杂行У胤乐诡w粒有序性，而顆粒有序性會限制邊坡破壞面的產(chǎn)生。如圖3（b）所示，一共生成的初始顆粒數(shù)量是50 000個。邊坡模型分條添加一種顏色，方便數(shù)值模擬直觀地得到沉降、崩塌的結(jié)果對比。當(dāng)“雨滴法”完成圖3（b）的模型后，在對其一側(cè)進行開挖，刪除多余的顆粒，進而得到如圖3（c）所示的最終的邊坡模型。顆粒間被設(shè)置了初始的粘結(jié)力，進而保證邊坡在沒有外力條件下是穩(wěn)定的。然后通過對模型中所有顆粒鍵施加較低的拉伸強度來降低粘結(jié)強度，進而在重力影響下顆粒鍵斷裂，顆粒因此破碎并造成滑動，形成崩塌、滑坡。在本次研究中還設(shè)定了采空區(qū)的工況（圖3（d）），在邊坡的坡腳處設(shè)置了采空區(qū)，研究不同位置采空區(qū)造成滑坡的最大沖擊力的結(jié)果?！安煽諈^(qū)—自由面”距離為圖中所示的S。相關(guān)的物理模型幾何參數(shù)如表2所示。

2 “災(zāi)害—結(jié)構(gòu)體”動力響應(yīng)結(jié)果分析

?

研究“災(zāi)害—結(jié)構(gòu)體”之間的相互作用對于實際的防災(zāi)減災(zāi)有著十分重要的意義。而災(zāi)害和結(jié)構(gòu)體之間的動力響應(yīng)又以力學(xué)的性能反應(yīng)出來。因此本次研究著重研究不同采空區(qū)開挖條件下滑坡災(zāi)害對結(jié)構(gòu)體的沖擊力，對采空區(qū)高度、采空區(qū)—自由面之間的距離、結(jié)構(gòu)體與坡腳的距離等一系列參數(shù)進行詳細的分析。

圖4顯示的是邊坡滑移面的生成過程。從圖4可以看出，邊坡的破壞從坡腳a點處開始，進而擴散到b點和c點，最終貫通整個坡面形成一個完整的邊坡滑移面。由于滑坡內(nèi)部的破壞機理十分復(fù)雜，本次研究不考慮開挖對邊坡內(nèi)部所造成的破壞，而是直接考慮滑移面產(chǎn)生后，滑坡災(zāi)害體運動過程中對結(jié)構(gòu)體造成的破壞。評判指標(biāo)主要為災(zāi)害體與結(jié)構(gòu)體之間的沖擊力。

2.1 “災(zāi)害體—結(jié)構(gòu)體”堆積結(jié)果影響分析

圖5顯示了不同坡腳條件下的滑坡初始情況和崩塌結(jié)果，從圖5可以看出圖5（a2）（坡角90°）情況下的崩塌量要明顯大于圖5（b2）（坡角70°）情況下的崩塌量。圖5（a2）條件下的平均下沉量也大于圖5（b2）條件下的平均下沉量。不僅如此，圖5（a2）情況下的滑動距離L同樣也大于圖5（b2）的情況。這是因為坡角的增大會使得潛在的滑坡量會隨之增加，當(dāng)顆粒鍵力逐漸減小時，邊坡在重力的作用下會隨之產(chǎn)生崩塌，而邊坡面越陡，潛在滑動塊的勢能就越大，因此其滑動的距離L也大于坡面不陡的情況。

圖5的（c1）顯示的是坡角70°，擋墻距離坡腳5 m時的災(zāi)害堆積結(jié)果，α為最終堆積角度，Hd為災(zāi)害體在結(jié)構(gòu)體接觸面的堆積高度。通過對比圖5（c1）和圖5（c2）以及圖5（d1）和圖5（d2）可以得出，隨著擋墻距離的增大，堆積角度α隨之減小，災(zāi)害體在結(jié)構(gòu)體接觸面上的堆積高度Hd也隨之減小。通過對比圖5（c1）和圖5（d1）以及圖5（c2）和圖5（d2）可以得出，隨著邊坡角度的增大，堆積角度α隨之增大，災(zāi)害體在結(jié)構(gòu)體接觸面上的堆積高度Hd也隨之增大。

2.2 采空區(qū)到邊坡自由面距離對災(zāi)害沖擊力的影響

如圖6（a）所示，顯示了在不同采空區(qū)離邊坡自由面距離的情況下，災(zāi)害體對結(jié)構(gòu)體的沖擊力隨時間步長的演化規(guī)律。從圖中同樣可以看出主滑體的沖擊力要遠遠大于崩塌和滾石所造成的沖擊力，此外S=20 m時的沖擊力的值幾乎在各個時段（除了1.7×106～2.1×106這段時步）都要大于其他各個工況的沖擊力值。

圖6（b）中給出了沖擊力隨不同采空區(qū)距邊坡自由面距離的演化規(guī)律，分別給出了災(zāi)害體與結(jié)構(gòu)體之間的最大沖擊力以及平均沖擊力的演化規(guī)律。從圖中可以很明顯地看出，最大沖擊力的值先是隨著S值的增加而增加，當(dāng)S值增加到某個值之后，繼續(xù)增加S的值會導(dǎo)致最大沖擊力的略微下降。同樣的規(guī)律可以反映到其平均沖擊力上?？梢钥闯鲈O(shè)置不同的采空區(qū)距邊坡自由面的距離會對災(zāi)害的動力性能造成不同的結(jié)果。而且，S值的變化會影響災(zāi)害的崩塌、滾石的沖擊力的演化規(guī)律，使其呈現(xiàn)主滑體的沖擊力的演化規(guī)律（最大沖擊力）。

2.3 滑坡和結(jié)構(gòu)體相互作用動力分析

圖7顯示的是70°坡角條件下的滑坡的最終堆積結(jié)果，其中采空區(qū)距離坡面的距離為10 m，采空區(qū)的采高為5 m。相對于沒有采空區(qū)的情況（圖5），圖7條件下的災(zāi)害體更加碎屑化，破壞區(qū)域更大，堆積范圍更廣。從圖中可以看出，堆積區(qū)B的面積基本和滑源區(qū)A的面積相吻合。除此之外，還存在孤立的堆積區(qū)C，與堆積區(qū)B之間存在一定的距離。

圖8顯示的是災(zāi)害體對結(jié)構(gòu)體的沖擊力（包括最大沖擊力和平均沖擊力）隨采空區(qū)高度增加而變化的曲線圖。從圖8（a）可以看出，最大沖擊力隨著采空區(qū)高度的增加而增長，并在15 m時達到最大值，之后基本穩(wěn)定在最大值附近；平均沖擊力的變化規(guī)律類似于最大沖擊力的規(guī)律，其峰值同樣出現(xiàn)在采高15 m處，之后盡管增加采高，其平均沖擊力值也是維持在該值附近。圖8（b）的最大沖擊力和平均沖擊力的演化規(guī)律和圖8（a）中的相似，同樣是隨著采高的增加而增加，并且維持在一個峰值附近不再變化。不同的是圖8（b）中的最大沖擊力和平均沖擊力的峰值都要小于圖8（a）中的相應(yīng)的峰值，這是因為隨著結(jié)構(gòu)體離坡腳的距離變遠，災(zāi)害體在滑動過程中會存在能量的損耗。

圖8（c）中的變化規(guī)律和圖8（a）和圖8（b）中的略有不同，其最大沖擊力維持在一個較大的值附近波動，與此同時，其平均沖擊力的規(guī)律也是如此。這是因為結(jié)構(gòu)體所處的距離決定的。如圖7所示，當(dāng)結(jié)構(gòu)體位于堆積區(qū)A以內(nèi)，可以攔截到所有的災(zāi)害體；而當(dāng)結(jié)構(gòu)體位于堆積區(qū)A之外，只有部分的崩塌體與結(jié)構(gòu)體相撞擊。該部分的崩塌體只占滑坡體災(zāi)害的極少比例，因而對結(jié)構(gòu)體的沖擊力也小得多。

3 結(jié)語

本次研究使用了二維離散元方法來對不同采高條件下的邊坡失穩(wěn)、滑動后的堆積、與結(jié)構(gòu)體撞擊等過程進行了分析。著重研究了不同采高條件下，滑坡災(zāi)害與結(jié)構(gòu)體之間的相互作用。研究過程分為3個部分進行，首先，通過對前人的工作的總結(jié)、對比來確定本次研究將要使用的微觀力學(xué)參數(shù)；其次，通過“雨滴法”在PFC2D中構(gòu)建二維離散元模型；最后，基于所選的參數(shù)和搭建的物理模型對所提出的問題進行詳細的研究。

通過二維離散元數(shù)值模擬可以得出，在無結(jié)構(gòu)體攔截的條件下，災(zāi)害體的堆積結(jié)果和邊坡的角度有關(guān)，角度越大，其滑程越遠；當(dāng)有結(jié)構(gòu)體攔截時，災(zāi)害體的堆積結(jié)果不僅和邊坡的角度有關(guān)，而且和其結(jié)構(gòu)體距離坡腳的距離有關(guān)。開挖條件下的滑坡體的堆積區(qū)主要分成2個部分（圖6），即主要的堆積區(qū)B和次要的堆積區(qū)C。當(dāng)結(jié)構(gòu)體位于主要的堆積區(qū)B內(nèi)時，災(zāi)害體對結(jié)構(gòu)體的最大沖擊力、平均沖擊力都會隨著采空區(qū)高度的增加而增加，并且最終穩(wěn)定在一個固定值不變；而當(dāng)結(jié)構(gòu)體位于主要的堆積區(qū)C內(nèi)時，災(zāi)害體對結(jié)構(gòu)體的沖擊力則幾乎維持一個范圍內(nèi)波動，不會出現(xiàn)之前的變化規(guī)律。

［1］ 畢鈺璋，付躍升，何思明，等.牛眠溝地震滑坡碎屑化全過程離散元模擬［J］.中國地質(zhì)災(zāi)害與防治學(xué)報，2015（3）:17-25.Bi Yuzhang，F(xiàn)u Yuesheng，He Siming，et al.Simulation of the whole process of Niumiangou creek rock avalanche triggered by the earthquake using a distinct element method［J］.The Chinese Journal of Geological Hazard and Control，2015（3）:17-25.

［2］ 畢鈺璋，何思明，付躍升，等.基于離散元方法的高速遠程滑坡碎屑流新型防護結(jié)構(gòu)［J］.山地學(xué)報，2015（5）:560-570.Bi Yuzhang，He Siming，F(xiàn)u Yuesheng，et al.Simulation of the dynamic response of new type rock avalanche impact defense structure and the mechanism of energy dissipation base on DEM［J］.Journal of Mountain Science，2015（5）:560-570.

［3］Mcsaveney M J.Recent rockfalls and rock avalanches in Mount Cook national park，New Zealand［J］.Reviews in Engineering Geology，2002，15:35-70.

［4］Valentino R，Barla G，Montrasio L.Experimental analysis and micromechanical modelling of dry granular flow and impacts in laboratory flume tests［J］.Rock Mechanics and Rock Engineering，2008，41（1）:153-177.

［5］ 崔 杰，王蘭生，王 衛(wèi)，等.采空區(qū)邊坡變形破裂演化機制研究［J］.采礦與安全工程學(xué)報，2008（4）:409-414.Cui Jie，Wang Lansheng，Wang Wei，et al.Deformation and fracturing mechanism of goaf slope［J］.Journal of Mining&Safety Engineering，2008（4）:409-414.

［6］Shao Yong，Yan Changhong，Xu Baotian.The stability analysis of mountain under intensive goaf［J］.Electronic Journal of Geotechnical Engineering，2016，20（4）:1181-1196.

［7］Shikai F A N.A discussion on the slope stability on the goaf［J］.Resources Environment and Engineering，2006，20:617-627.

［8］Dong L，Li X.An application of grey-general regression neural network for predicting landslide deformation of Dahu Mine in China［J］.Advanced Science Letters，2012，6（1）:577-581.

［9］Yang X J，Hou D G，Hao Z L，et al.Fuzzy comprehensive evaluation of landslide caused by underground mining subsidence and its monitoring［J］.International Journal of Environment and Pollution，2016，59（2-4）:284-302.

［10］Chan N W.Responding to landslide hazards in rapidly developing Malaysia:a case of economics versus environmental protection［J］.Disaster Prevention and Management:An International Journal，1998，7（1）:14-27.

［11］Fujiaswa K，Kasai M，Okuda S，et al.Construction of a numerical analysis model to evaluate interaction between a landslide and a tunnel［C］//The Japan Landslide Society 48th Annual Congress.Niigata，Japan［.s.n.］2009:71-72.

［12］柴紅保，曹 平，柴國武，等.采空區(qū)對邊坡穩(wěn)定性的影響［J］.中南大學(xué)學(xué)報:自然科學(xué)版，2010，41（4）:1528-1534.Chai Hongbao，Cao Ping，Chai Guowu，et al.Influence of goaf on slope stability［J］.Journal of Central South University:Science and Technology ，2010，41（4）:1528-1534.

［13］宋偉東，杜建華，謝正平.大冶鐵礦深凹露天開采最終邊坡穩(wěn)定性分析［J］.北京科技大學(xué)學(xué)報，2005，27（4）:385-389.Song Weidong，Du Jianhua，Xie Zhengpinig.Stability of the finished slope of Daye iron open-pit-mine［J］.Journal of University and Technology Beijing，2005，27（4）:385-389.

［14］丁桂伶.采空區(qū)上方邊坡穩(wěn)定性及滑坡危險性分析［J］.中國地質(zhì)災(zāi)害與防治學(xué)報，2012，23（2）:38-43.Ding Guiling.Stability and landslide risk analysis of slope above goaf［J］.The Chinese Journal of Geological Hazard and Control，2012，23（2）:38-43.

［15］Teufelsbauer H，Wang Y，Pudasaini S P，et al.DEM simulation of impact force exerted by granular flow on rigid structures［J］.Acta Geotechnica，2011，6（3）:119.

［16］Zanuttigh B，Lamberti A.Experimental analysis of the impact of dry avalanches on structures and implication for debris flows［J］.Journal of Hydraulic research，2006，44（4）:522-534.

［17］畢鈺璋，何思明，李新坡，等.約束條件下粗細混合顆粒動力機理分析［J］.巖土工程學(xué)報，2016，38（3）:529-536.Bi Yuzhang，He Siming，Li Xinpo，et al.Kinetic mechanism of mixed particles under constraint conditions［J］.Chinese Journal of Geotechnical Engineering，2016，38（3）:529-536.

［18］Cundall P A，Strack O D.A discrete numerical model for granular assemblies［J］.Geotechnique，1979，29（1）:47-65.

［19］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.

［20］Utili S，Nova R.DEM analysis of bonded granular geomaterials［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2008，32（17）:1997-2031.

［21］Villard P，Chevalier B，Le Hello B，et al.Coupling between finite and discrete element methods for the modelling of earth structures reinforced by geosynthetic［J］.Computers and Geotechnics，2009，36（5）:709-717.

［22］Zhang L，Thornton C.A numerical examination of the direct shear tes［tJ］.Geotechnique，2007，57（4）:343-354.

［23］Bi Y，He S，Li X，et al.Geo-engineered buffer capacity of two-layered absorbing system under the impact of rock avalanches based on Discrete Element Method［J］.Journal of Mountain Science，2016，13（5）:917-929.

［24］Koizumi Y，Lee J，Yokota Y，et al.Numerical analysis of landslide behavior induced by tunnel excavation［C］//ISRM International Symposium-EUROCK 2010.Lausanne，Switzerland:International Society for Rock Mechanics，2010.