堰塞體形成全過程的連續(xù)離散耦合數(shù)值模擬

王 葉，周 偉，馬 剛，陳 遠(yuǎn)，常曉林

(1. 武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 430072；2. 武漢大學(xué)水工巖石力學(xué)教育部重點實驗室，武漢 430072)

自然界中堰塞體的形成因素很多，根據(jù)統(tǒng)計，有90%以上是地震誘發(fā)滑坡形成的[1]，滑坡體堆積河谷形成堰塞體，一旦失穩(wěn)將引發(fā)下游洪水災(zāi)害，因此深入研究地震作用下的堰塞體形成機制，對于評估堰塞體的穩(wěn)定性、堰塞湖潛在次生災(zāi)害可能性及其損失具有重要意義。

目前對于堰塞體形成機制的研究方法主要有有限元法(FEM)[2,3]、離散元法(DEM)[4-7]、非連續(xù)變形分析方法(DDA)[8,9]，連續(xù)-離散耦合分析方法(FDEM)[11-17]等。徐文杰等[2]在對肖家橋滑坡區(qū)詳細(xì)的地質(zhì)調(diào)查基礎(chǔ)上，結(jié)合動力有限元分析進行，再現(xiàn)滑坡三維空間失穩(wěn)過程；劉傳正等[3]通過有限差分方法研究了紅石巖山體在地震的動力響應(yīng)特征和崩塌的形成機理。在強震作用下，巖土體發(fā)生局部大變形甚至形成強間斷面(巖石拉斷、沿結(jié)構(gòu)面滑移)，此時基于連續(xù)介質(zhì)力學(xué)的傳統(tǒng)有限元方法和有限差分法受小變形假設(shè)的限制而很難適用。曹琰波等[4]采用離散單元法，對唐家山滑坡由變形累積到破壞滑動的全過程進行了模擬，研究了地震作用下順層巖質(zhì)滑坡全過程；Wu等[5]采用離散元方法分析一個地震誘發(fā)的大規(guī)?；?，研究滑坡體滑動過程和滑坡后的形態(tài)等；焦玉勇等[6]采用離散單元法模擬了邊坡失穩(wěn)過程；鄔愛清等[8]采用非連續(xù)變形分析方法，以唐家山滑坡形成的堰塞壩形態(tài)和位置作為目標(biāo)函數(shù)，對唐家山滑坡過程進行復(fù)演，探討地震荷載下高速滑坡的形成機制；Zhang等[9]提出附加接觸模型的改進的DDA方法，研究有黏聚力材料的邊坡模型，評估邊坡節(jié)理面處的抗剪強度；Huang等[10,11]采用考慮節(jié)理動摩擦退化的改善非連續(xù)變形分析方法(IDDA)研究地震誘發(fā)的滑坡體的速度、位移等動力學(xué)行為。

DEM和DDA能模擬邊坡失穩(wěn)及滑動全過程，但是它們對裂紋、應(yīng)力、應(yīng)變的表征不夠直觀，因此有學(xué)者開始采用連續(xù)離散耦合分析方法來研究邊坡失穩(wěn)及滑動全過程。Barla等[12]采用FDEM對位于意大利Alpetto礦山的2處邊坡進行數(shù)值模擬，與早前的研究進行比較，驗證了FDEM模擬邊坡失穩(wěn)的可行性；Antolini等[13]采用FDEM方法對Torgiovannetto di Assisi滑坡體進行數(shù)值模擬，研究滑坡體的啟動機制和演化過程，為滑坡體的風(fēng)險評估提供參考。Piovano等[14]研究了意大利Aosta山谷的Beauregard滑坡，表明運用FDEM方法能夠預(yù)測深層滑坡的典型不穩(wěn)定面和破壞機理；常曉林等[15]采用FDEM模擬了在暴雨等情況下邊坡失穩(wěn)滑移的過程，較好地模擬了抗剪斷參數(shù)降低時邊坡從小變形-大變形-局部滑動-整體滑移的全過程失效；Zhou等[16]基于FDEM方法研究了節(jié)理傾角對邊坡破壞模式、運動過程、破壞后的堆積形態(tài)等的影響。O. K. Mahabadi等[17]應(yīng)用基于FDEM的Y-Geo代碼較好地模擬了山崖侵蝕的過程。但是FDEM進行滑坡全過程模擬大多限于二維，對三維滑坡體研究甚少。

在連續(xù)-離散耦合分析方法中引入界面單元，采用內(nèi)聚力模型能有效地捕捉加載過程中裂紋的擴展路徑，模擬巖體漸進破壞過程，充分發(fā)揮連續(xù)介質(zhì)力學(xué)方法和非連續(xù)介質(zhì)力學(xué)方法各自的優(yōu)勢[15-18]，實現(xiàn)巖體從連續(xù)狀態(tài)到非連續(xù)狀態(tài)的轉(zhuǎn)化過程。本文采用FDEM，結(jié)合紅石巖堰塞體，對三維滑坡由變形破壞到整體滑動的全過程進行模擬，以研究地震作用下巖質(zhì)滑坡的失穩(wěn)破壞過程，直觀地捕捉滑坡體從開裂破壞到最終的堆積狀態(tài)、滑坡體失穩(wěn)過程中的速度分布情況等。

1 基于內(nèi)聚力模型的巖石開裂擴展模擬

在巖體連續(xù)-離散耦合分析方法中[12-18]，巖石被離散為由實體單元和無厚度界面單元組成的系統(tǒng)，每個實體單元均為單獨的離散塊體。界面單元失效之前，離散塊體之間通過界面單元“黏結(jié)”成一個整體，以模擬巖石材料的連續(xù)變形。實體單元只發(fā)生彈性變形，損傷和斷裂僅發(fā)生在界面單元上。界面單元的破壞準(zhǔn)則為帶拉斷的Mohr-Coulomb準(zhǔn)則，界面單元的應(yīng)力狀態(tài)滿足破壞準(zhǔn)則后，采用基于斷裂能的線性損傷演化模型，損傷因子達(dá)到1.0后完全失效。界面單元失效后從模型中刪除，兩側(cè)實體單元發(fā)生接觸關(guān)系，采用罰函數(shù)法進行接觸分析。

采用內(nèi)聚力模型模擬巖石材料的界面開裂行為，通過界面單元的失效模擬裂紋的萌生、擴展、交匯和貫通。在張拉應(yīng)力狀態(tài)下，簡化的內(nèi)聚力區(qū)應(yīng)力分布見圖1。真實裂縫的尖端應(yīng)力為零，內(nèi)聚力區(qū)尖端應(yīng)力為材料的抗拉強度，從真實裂縫尖端到內(nèi)聚力區(qū)尖端應(yīng)力逐漸提高。二維界面單元具有2個積分點，均位于單元的中平面上。采用共節(jié)點的方式實現(xiàn)界面單元與相鄰實體單元之間力、位移的傳遞。圖2給出局部模型的有限元網(wǎng)格拓?fù)浣Y(jié)構(gòu)圖，裂縫可從布置界面單位的位置擴展。

圖1 準(zhǔn)脆性材料的內(nèi)聚力區(qū)應(yīng)力分布Fig.1 Illustration of cohesive zone for quasi brittle material

圖2 網(wǎng)格拓?fù)鋱D(實體單元面積收縮便于顯示界面單元，實際界面單元厚度為0)Fig.2 Mesh topology (the elastic elements are shrunk for illustration purpose, the actually element thickness is zero)

在界面單元出現(xiàn)損傷之前，假定其本構(gòu)關(guān)系是線彈性：

(1)

式中：tn、ts分別為界面單元的法向應(yīng)力、切向應(yīng)力；Kij為剛度矩陣，本文計算時未考慮應(yīng)力、應(yīng)變的法向分量與切向分量之間的耦合關(guān)系，因此Kns=0 ；εn、εs分別為界面單元的法向應(yīng)變、切向應(yīng)變。

法向應(yīng)變、切向應(yīng)變定義為：

(2)

式中：δn、δt分別為界面單元沿著法向和切向產(chǎn)生的位移；T0為界面單元的本構(gòu)厚度。

本文采用無厚度的界面單元模擬材料開裂，倘若取T0為幾何厚度，將會導(dǎo)致應(yīng)變奇異。為了消除應(yīng)變奇異，在進行界面單元的應(yīng)變和應(yīng)力計算時采用的是單元的本構(gòu)厚度。通常定義T0=1，則界面單元的應(yīng)變與相應(yīng)方向上的位移大小相等。界面單元的線彈性本構(gòu)方程可寫成：

(3)

式中：kn、ks分別為界面單元的法向剛度、切向剛度。

為了模擬裂紋的萌生、擴展過程，需要定義合理的裂縫起裂準(zhǔn)則，考慮巖石等準(zhǔn)脆性材料的破壞往往是由拉應(yīng)力和剪應(yīng)力共同作用導(dǎo)致的，本文采用二次應(yīng)力準(zhǔn)則判斷界面單元是否開裂，采用帶拉斷的Mohr-Coulomb準(zhǔn)則計算巖石材料的抗剪強度：

(4)

式中：c為巖石材料的黏聚力；φ為巖石材料的內(nèi)摩擦角；t0n為抗拉強度；t0s為抗剪強度。

二次應(yīng)力準(zhǔn)則：當(dāng)各個方向的應(yīng)力與其相應(yīng)的臨界應(yīng)力的比值的平方和達(dá)到1時，界面單元即開始出現(xiàn)損傷：

(5)

在損傷演化階段，界面單元的本構(gòu)方程為：

(6)

ts=(1-D)ksδs

式中：D為無量綱的損傷因子，當(dāng)D=0時，表明界面單元未出現(xiàn)損傷；當(dāng)D=1時，表明界面單元完全失效，失去承載能力。

損傷因子的表達(dá)式如下：

(7)

式中：t0eff為達(dá)到破壞準(zhǔn)則時的等效應(yīng)力；δ0m為達(dá)到破壞準(zhǔn)則時的等效位移；δmaxm為加載歷史中的最大等效位移。

在大多數(shù)情況下，界面單元處于在混合加載狀態(tài)下，同時發(fā)生法向和切向的變形，因此需要定義界面單元的等效應(yīng)力teff和等效位移δm：

(8)

通過定義界面單元在損傷過程中耗散的能量，即可控制界面單元的損傷演化過程。在本文中，耗散能量也就是通常所說的斷裂能，它等于應(yīng)力-分離曲線下面所包含的面積大小。本文模擬時采用基于線性軟化的Benzeggagh-Kenane準(zhǔn)則，界面單元的應(yīng)力-分離曲線見圖3。

圖3 界面單元的應(yīng)力-分離曲線(Ⅰ、Ⅱ)Fig.3 Constitutive relations of the cohesive elements

Benzeggagh-Kenane準(zhǔn)則的表達(dá)式如下：

GT=Gs+GⅠ

(9)

式中：Gshear為界面單元在剪切荷載作用下的斷裂能；GT為界面單元在混合荷載作用下的斷裂能；η為材料常數(shù)，通常由彎曲試驗測得，本文取2。

2 紅石巖滑坡體FDEM三維模擬

2.1 工程概況

紅石巖堰塞壩形成于2014年發(fā)生的云南魯?shù)椤?·03”地震?！?·03”地震發(fā)生后，左岸滑坡堆積物表層松動并向河床滑動，右岸山體產(chǎn)生大規(guī)?；?，在河床形成堰塞湖。根據(jù)現(xiàn)場調(diào)查，堰塞體頂部呈馬鞍形，頂部左岸高，右岸低，右岸邊緣為滑坡巖石堆積體。堰塞體組成松散，最大厚度約103 m。頂部橫河向最低高程點1 222 m，堰塞體左岸最高點為1 270 m，上游迎水面平均坡比約1∶2.5，下游面平均坡比約1∶5.5。順河向底寬約910 m，沿1 222 m高程壩軸線長度約307 m。右岸滑坡后，地形發(fā)生了較大變化，上部滑床后緣為陡崖，高度估計約150～200 m，中部形成一個橫河向近水平且傾向下游的斜面地形，在此斜面以下為陡崖。左岸坡雖有表面松散并向下滾落，但體量不大，總的地形未有大的改變。堰塞體物質(zhì)組成較為復(fù)雜，主要為右岸邊坡崩滑后形成的崩塌堆積物。

2.2 FDEM模擬所需參數(shù)率定

采用連續(xù)-離散耦合分析方法進行數(shù)值模擬時，需要對所用參數(shù)進行率定。采用二維平面應(yīng)變模型進行雙軸數(shù)值壓縮試驗來率定模擬參數(shù)，通過不同圍壓下的偏應(yīng)力應(yīng)變曲線得到應(yīng)力莫爾圓從而計算出巖土體的宏觀抗剪強度參數(shù)。用于參數(shù)率定的數(shù)值試樣尺寸為80 mm×160 mm(見圖4)。剛性加載板與試樣之間定義接觸關(guān)系，采用位移控制方式加載施加于上、下剛性加載板，圍壓施加在左右加載板上。剛性加載板與數(shù)值試樣之間的摩擦系數(shù)取為0.1。

圖4 雙軸壓縮試驗試樣尺寸及加載示意圖Fig.4 Size of numerical experiments of biaxial compression and illustration of loading

紅石巖滑坡體各土層的物理力學(xué)參數(shù)見表1。在FDEM數(shù)值模擬中，界面單元的參數(shù)包括剛度、抗拉強度、黏聚力c、內(nèi)摩擦角φ、斷裂能。通過反復(fù)試算，使得數(shù)值試驗得到的力學(xué)參數(shù)與表1所示的室內(nèi)試驗值接近。表2對比了數(shù)值試驗得到的各土層宏觀黏聚力、內(nèi)摩擦角與室內(nèi)試驗各土層抗剪強度，相對誤差都在10%以內(nèi)。最終確定的FDEM模擬參數(shù)見表3。

表1 紅石巖巖體室內(nèi)試驗物理力學(xué)參數(shù)Tab.1 The physical and mechanical parameters of laboratory tests for Hongshiyan rocks

表2 紅石巖巖體力學(xué)參數(shù)與數(shù)值試樣計算力學(xué)參數(shù)結(jié)果對比Tab.2 The comparison of mechanical parameters ofrocks for Hongshiyan rocks between laboratorytests and numerical experiments

表3 紅石巖巖體數(shù)值模擬界面單元細(xì)觀參數(shù)Tab.3 The mesoscopic parameters of cohesive interface elements for numerical simulation of Hongshiyan rocks

2.3 計算模型

根據(jù)震前邊坡原地形圖進行滑坡前的還原。由于地震影響范圍較廣，為減小邊界反射效應(yīng)，將模型邊界范圍取得足夠大。三維計算模型中，x軸方向為橫河向，由左岸指向右岸為正；y軸方向為順河向，向下游為正；z軸方向為豎直向，向上為正。橫河向長1 427.1 m，縱河向1 620.5 m，豎直向970 m。

地震荷載作用下，只考慮邊坡底部的彈性作用，消除邊坡底部對地震的放大作用，采用有剛度無質(zhì)量方案進行分析，故底部約束為固定約束，側(cè)向邊界為法向約束。在盡量考慮實際地形地貌和地質(zhì)結(jié)構(gòu)的前提下，同時考慮計算效率，三維整體模型共剖分為341 651 個單元和376 002 個節(jié)點，其中實體單元為215 271個，界面單元為126 380 個。本文重點關(guān)注滑坡形成堰塞體的整個過程，所以只在滑坡體部分插入界面單元。整個FDEM模擬分2步進行，首先進行地應(yīng)力平衡，得到初始應(yīng)力場，再施加地震加速度荷載，進行地震分析。

在計算模型底部輸入3個方向的地震加速度，模擬滑坡體在地震作用下的啟動、滑動到最終堆積直至形成堰塞體的全過程。加速度時程曲線見圖5(a)、5(b)、5(c)，地震作用過程歷時20 s(10～30 s)，總的作用時間為100 s。

圖5 地震加速度動力時程曲線Fig.5 Time-history curves of earthquake acceleration

3 堰塞體形成全過程分析

圖6為堰塞體最大剖面對比圖(對比范圍相同)，實際堰塞體最大堆積厚度約為103 m，模擬的對應(yīng)縱剖面最大堆積厚度為101.97 m，由于模型的簡化和模擬范圍限制，2者的堆積形態(tài)稍有不同，但大致相似。

圖6 堰塞體最大縱剖面對比Fig.6 Comparison for longitudinal section of landside dam

定義堰塞體堆積最大寬度為Lmax，堆積最大高度為Hmax，寬高比Lmax/Hmax，用這3個參數(shù)表征堆積形態(tài)。表4統(tǒng)計了4個典型剖面堆積形態(tài)實測值和FDEM計算值及2者的相對誤差(百分制)。由于模型簡化和網(wǎng)格劃分關(guān)系，表4中個別計算值與實測值有較大差別，故其相對誤差較大。圖7對比了4個典型橫剖面實測與FDEM模擬的堆積形態(tài)。從表4和圖7可以看出模擬得到的橫剖面堆積形態(tài)與實測堆積形態(tài)基本一致，進一步驗證了本文FDEM模擬邊坡失穩(wěn)、滑動形成堰塞體的合理性。

表4 典型橫剖面堆積形態(tài)實測值與計算值對比Tab.4 The comparison of accumulation state of typical sections between measured value and calculated value

滑坡體不同時刻的滑動狀態(tài)見圖8。選取4個時刻滑坡典型剖面g-g′的形態(tài)[見圖9(a)～圖9(d)] 研究地震誘發(fā)滑坡失穩(wěn)過程運動特征。從圖8和圖9可以看出，在地震初期由于受地震作用，表層巖體首先出現(xiàn)裂縫然后貫通形成整體滑裂面；隨著地震作用，滑坡體內(nèi)部損傷不斷累積，t=20 s左右時，整個滑坡體與下覆基巖出現(xiàn)明顯的運動分離[見圖9(a)]；然后在地震和重力共同作用下滑坡體高速下滑，且在下滑過程中滑坡體不斷破碎、解體，受對岸山體的阻擋，與之發(fā)生高速撞擊，進一步破碎、解體、堆積、堵塞河道，在60 s后滑坡體逐漸穩(wěn)定，逐漸堆積密實，100 s時基本靜止，形成堰塞體。

圖7 典型橫剖面實測與計算堆積形態(tài)對比Fig.7 The comparison of accumulation state of typical sections between actual measurement and calculation

圖8 滑坡體不同時刻滑動形態(tài)Fig.8 Sliding patterns of landslide at different time

圖9 滑坡體典型剖面g-g′不同時刻滑坡形態(tài)Fig.9 Sliding patterns for a typical section named “g-g′” at different time

圖10為滑坡體典型剖面g-g′不同時刻速度云圖。t=20 s左右滑坡體速度開始增大，滑塊底部首先達(dá)到8 m/s左右，在地震和重力作用下，滑坡體的速度持續(xù)增大，t=30 s左右高速下滑，達(dá)到20 m/s左右，速度的分布受邊坡坡度的影響，由于上部坡度較陡而中部坡度緩，下部與中部坡度相近，但下部是臨空面，而中部有下部巖體阻擋，故巖體失穩(wěn)時滑動速度中部較小，上部和下部速度較大；t=40 s滑坡體已經(jīng)解體，滑坡體在下滑過程中積聚很大的動能，故此時分布有最大滑坡速度，撞擊對岸巖體后速度迅速改變方向，不斷堆積河谷，直至逐漸穩(wěn)定，60 s左右速度逐漸減小為零。

圖10 滑坡體典型剖面g-g′不同時刻速度云圖Fig.10 Cloud pictures of velocity of a typical section named “g-g′” at different time

圖11顯示了系統(tǒng)的能量演化曲線(E為能量值，Emax為相應(yīng)能量的最大值)，包括系統(tǒng)總動能歷時曲線、系統(tǒng)摩擦耗散能歷時曲線和破碎耗散能歷時曲線。在歷時20 s左右，系統(tǒng)的總動能和摩擦耗散能出現(xiàn)突變，邊坡開始失穩(wěn)，滑坡體高速下滑，在40 s左右由于撞擊對岸巖體，動能開始迅速減小。在t=60 s，各能量基本保持不變，滑坡堆積體達(dá)到穩(wěn)定狀態(tài)。從摩擦耗散能和破碎耗散能歷時曲線可知，隨著滑坡體的失穩(wěn)和滑動，失效界面單元不斷增加，導(dǎo)致滑坡體由連續(xù)體向非連續(xù)體轉(zhuǎn)換，并最終形成松散的堆積體。

圖11 系統(tǒng)能量演化曲線Fig.11 Evolution curves of system energy

4 結(jié) 論

(1)將FDEM應(yīng)用于紅石巖的三維邊坡失穩(wěn)模擬，重現(xiàn)了堰塞體形成的全過程。在地震作用10 s左右，系統(tǒng)動能出現(xiàn)突變，邊坡開始失穩(wěn)，地震總的作用時間20 s后，滑坡體高速下滑，進而受到對岸撞擊，不斷解體并堆積在河谷。通過FDEM分析能夠直觀地捕捉滑坡體各個時刻的滑動狀態(tài)、速度分布和滑坡體的最終堆積形態(tài)。

(2)將最大縱剖面和4個典型橫剖面的數(shù)值模擬的堆積形態(tài)與現(xiàn)場實際堆積形態(tài)進行對比，2者吻合程度較高。對比了表征堆積形態(tài)的參數(shù)，基本一致，進一步驗證了FDEM模擬邊坡失穩(wěn)和滑動堆積形成堰塞體全過程的合理性。

(3)FDEM三維模擬結(jié)果表明，紅石巖堰塞體形成機制為地震誘發(fā)巖土體內(nèi)部損傷和變形累積、結(jié)構(gòu)面貫通、滑坡體高速下滑并伴隨破壞、解體、撞擊對面山體，堆積河道、不斷密實，形成堰塞體。

□

[1] Costa J E, Schuster R L. The formation and failure of natural dams [J]. Geological Society of America Bulletin, 1988,100(7):1 054-1 068.

[2] 徐文杰，陳祖煜，何秉順，等. 肖家橋滑坡堵江機制及災(zāi)害鏈效應(yīng)研究[J]. 巖石力學(xué)與工程學(xué)報，2010，29 (5)：933-942.

[3] 劉傳正，葛永剛，江興元，等. 魯?shù)榈卣鸺t石巖崩塌觸發(fā)機理分析[J]. 防災(zāi)減災(zāi)工程學(xué)報，2016，36(4)：601-607.

[4] 曹琰波，戴福初，許 沖，等. 唐家山滑坡變形運動機制的離散元模擬[J]. 巖石力學(xué)與工程學(xué)報，2011，30(增1)：2 879-2 887.

[5] Wu Jianhong, Lin Jeenshang, Chen Chaoshi. Dynamic discrete analysis of an earthquake-induced large-scale landslide[J]. International Journal of Rock Mechanics & Mining Sciences, 2009,(46):397-407.

[6] 焦玉勇，葛修潤，劉泉聲,等. 三維離散單元法及其在滑坡分析中的應(yīng)用[J]. 巖土工程學(xué)報，2000，20(1)：101-104.

[7] Zhou Wei, Lai Zhiqiang, Ma Gang. Effect of base roughness on size segregation in dry granular flows[J].Granular Matter, 2016,(18):83.

[8] 鄔愛清，林紹忠，馬貴生，等. 唐家山堰塞壩形成機制DDA模擬研究[J]. 人民長江，2008，39(22)：91-95.

[9] 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):533-545.

[10] Huang Da, Song Yixiang, Cen Duofeng, et al. Numerical modeling of earthquake-induced landslide using an improved discontinuous deformation analysis considering dynamic friction degradation of joints[J]. Rock Mech Rock Eng, 2016,49:4 767-4 786.

[11] Barrett P J. The shape of rock particles, a critical review [J]. Sedimentology, 1980,27(3):291-303.

[12] Marco Barla, Giovanna Piovano, Giovanni Grasselli. Rock slide simulation with the combined finite-discrete element method [J]. Geomechanics, 2012,12(6):711-721.

[13] Francesco Antolini, Marco Barla, Andrea Giorgetti, et al. Combined finite-discrete numerical modeling of runout of the torgio-vannetto di assisi rockslide in central italy [J]. International Journal of Geo-mechanics, 2016,16(6):1-16.

[14] Giovanna Piovano, Francesco Antolini, Marco Barla, et al. Continuum-discontinuum modelling of failure and evolution mechanisms of deep seated land-slides[C]∥6th Int. Conf. on Discrete Element Method, Colorado School of Mines, Golden, CO, USA, 2013:295-300.

[15] 常曉林，胡 超，馬 剛，等. 模擬巖體失效全過程的連續(xù)-非連續(xù)變形體離散元方法及應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報，2011，30(10)：2 004-2 011.

[16] Zhou Wei, Yuan Wei, Ma Gang, et al. Combined finite-discrete element method modeling of rockslides [J]. Engineering Computations, 2016,33(5):1 530-1 559.

[17] O K Mahabadi, A Lisjak, A Munjiza, et al. Y-Geo: new combined finite-discrete element numerical code for geome-chanical applications[J]. Geomechanics, 2012,12(6):676-688.

[18] 馬 剛，周創(chuàng)兵，常曉林，等. 巖石破壞全過程的連續(xù)-離散耦合分析方法[J]. 巖石力學(xué)與工程學(xué)報，2011，30 (12)：2 444-2 454.