一種有限元轉(zhuǎn)化為顆粒離散元的方法及其應(yīng)用研究

馮 春，李世海，劉曉宇

（中國(guó)科學(xué)院力學(xué)研究所 流固耦合系統(tǒng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100190）

1 引 言

地質(zhì)災(zāi)害的形成過程是地質(zhì)體從連續(xù)介質(zhì)到非連續(xù)介質(zhì)的漸進(jìn)破壞過程。地質(zhì)體內(nèi)部裂紋的萌生、擴(kuò)展、交匯、貫通，表征了地質(zhì)災(zāi)害發(fā)展演化的不同階段[1]。地質(zhì)體中多裂紋擴(kuò)展過程的數(shù)值模擬是當(dāng)前國(guó)內(nèi)外學(xué)術(shù)界的研究熱點(diǎn)及難點(diǎn)。其中，有限元法及離散元法是兩類主要的數(shù)值模擬方法。

有限元法以連續(xù)介質(zhì)力學(xué)為基礎(chǔ)，能夠很好地描述地質(zhì)體的彈塑性變形，適用于判斷地質(zhì)體從穩(wěn)定到不穩(wěn)定的臨界點(diǎn)。為了分析地質(zhì)體中裂縫的擴(kuò)展演化過程，需在有限元方法中引入新的計(jì)算思路，如自適應(yīng)網(wǎng)格及擴(kuò)展有限元。自適應(yīng)網(wǎng)格根據(jù)計(jì)算獲得的裂紋起裂方向?qū)Τ跏季W(wǎng)格進(jìn)行局部細(xì)化及調(diào)整，使細(xì)化后的網(wǎng)格邊界沿著裂紋擴(kuò)展的方向，從而實(shí)現(xiàn)裂紋擴(kuò)展過程的模擬[2－3]。此類方法在計(jì)算復(fù)雜工程問題時(shí)，涉及大量網(wǎng)格的重新劃分及調(diào)整，算法穩(wěn)定性差，計(jì)算效率低下，且新網(wǎng)格下各類場(chǎng)信息不能完全有效繼承，往往會(huì)出現(xiàn)數(shù)值振蕩現(xiàn)象。擴(kuò)展有限元法通過引入跳躍函數(shù)實(shí)現(xiàn)有限元單元內(nèi)部的非連續(xù)變形及隱式裂紋的擴(kuò)展[4－6]，該方法具有嚴(yán)格的數(shù)學(xué)推導(dǎo)，但只能描述有限條裂紋的擴(kuò)展過程，且在模擬裂紋交匯、分叉及三維裂縫擴(kuò)展方面不盡如人意。

離散元法以非連續(xù)介質(zhì)力學(xué)為基礎(chǔ)，較為適合模擬地質(zhì)體破壞后的運(yùn)動(dòng)解體過程，塊體離散元及顆粒離散元是兩類主要的離散元方法。塊體離散元法用塊體表征連續(xù)介質(zhì)單元（剛體單元或可變形單元），用塊體間的交界面表征非連續(xù)單元，通過塊體邊界的斷裂實(shí)現(xiàn)裂紋擴(kuò)展過程的模擬[7－10]。此類方法的裂縫僅能沿著塊體邊界擴(kuò)展，因此，網(wǎng)格依賴性嚴(yán)重；且塊體界面剛度的物理意義并不明確。當(dāng)然，通過單元切割法可在一定程度上解決網(wǎng)格依賴性問題[11－12]，利用結(jié)構(gòu)層也可部分解決界面剛度的選取問題[1]。顆粒離散元法通過細(xì)觀顆粒簇的集合表述連續(xù)介質(zhì)的宏觀特性，通過細(xì)觀顆粒的斷裂、滑移反映宏觀裂縫的萌生、擴(kuò)展過程[13－15]。由于顆粒排布的隨機(jī)性，顆粒離散元中的裂縫可向任意方向擴(kuò)展。然而，此類方法無法準(zhǔn)確描述宏觀介質(zhì)的力學(xué)特性，細(xì)觀與宏觀參數(shù)的對(duì)應(yīng)關(guān)系往往需要通過大量的數(shù)值試驗(yàn)進(jìn)行標(biāo)定[16]。

為了充分發(fā)揮有限元與顆粒離散元各自的優(yōu)勢(shì)，將有限元與顆粒離散元進(jìn)行耦合計(jì)算是一條有效途徑。O?ate等[17]對(duì)有限元與顆粒離散元耦合的計(jì)算方法進(jìn)行了詳細(xì)論述，并將該方法成功應(yīng)用于單軸壓縮、巴西劈裂、樁體刺入等典型巖土工程問題中；Rojek等[18]基于文獻(xiàn)[17]提出了一種有限元與顆粒離散元局部區(qū)域重疊耦合的方法，通過漸變函數(shù)實(shí)現(xiàn)有限元與顆粒離散元的自然過渡，并將該方法成功應(yīng)用于隧道巖體的切割計(jì)算中。然而，上述方法中有限元與顆粒離散元的耦合屬于分區(qū)耦合，變形破壞區(qū)用顆粒離散元進(jìn)行描述，弱變形區(qū)用有限元進(jìn)行描述。由于在強(qiáng)變形區(qū)采用了顆粒離散元，該區(qū)域初始狀態(tài)下的宏觀應(yīng)力-應(yīng)變依然無法準(zhǔn)確描述。

因此，本文提出了一種有限元轉(zhuǎn)化為顆粒離散元的方法，數(shù)值模型首先用較粗的有限元網(wǎng)格進(jìn)行離散，當(dāng)達(dá)到破壞條件，刪除有限元單元，引入顆粒離散元，并利用顆粒離散元實(shí)現(xiàn)地質(zhì)體漸進(jìn)破壞過程的模擬。

2 有限元轉(zhuǎn)化為顆粒離散元的方法

2.1 主要思路

首先將研究區(qū)域用較粗的有限元單元進(jìn)行離散，每個(gè)單元按照虎克定律進(jìn)行應(yīng)力-應(yīng)變的計(jì)算。當(dāng)某個(gè)單元的應(yīng)力狀態(tài)滿足Mohr-Coulomb準(zhǔn)則或最大拉應(yīng)力準(zhǔn)則時(shí)，刪除該單元，并在該單元所在位置創(chuàng)建具有一定數(shù)目、隨機(jī)分布且微嵌入的顆粒系統(tǒng)。所創(chuàng)建顆粒的材料性質(zhì)、速度、位移、接觸力等信息根據(jù)插值函數(shù)從單元中繼承。而后利用顆粒間接觸力的演化自然實(shí)現(xiàn)宏觀裂紋的萌生、擴(kuò)展及貫通過程。具體的轉(zhuǎn)化方法如圖1所示。

圖1 有限元轉(zhuǎn)化為顆粒離散元的基本思路Fig.1 Basic idea for transiting FEM into particle DEM

2.2 有限元的求解

本文在有限元計(jì)算時(shí)采用顯式求解策略，主要包含節(jié)點(diǎn)合力計(jì)算及節(jié)點(diǎn)運(yùn)動(dòng)計(jì)算兩個(gè)部分。節(jié)點(diǎn)合力計(jì)算為

式中：F為節(jié)點(diǎn)合力；Fe為節(jié)點(diǎn)外力；Fd為節(jié)點(diǎn)變形力（由單元應(yīng)力貢獻(xiàn)）；Fc為阻尼力。

節(jié)點(diǎn)運(yùn)動(dòng)計(jì)算式為

式中：a為節(jié)點(diǎn)加速度；v為節(jié)點(diǎn)速度；Δu為節(jié)點(diǎn)位移增量；u為節(jié)點(diǎn)位移全量；m為節(jié)點(diǎn)質(zhì)量；Δt為計(jì)算時(shí)步?；谑剑?）、（2）的交替計(jì)算，即可實(shí)現(xiàn)有限元的顯式求解過程。

本文采用增量法進(jìn)行單元應(yīng)力及節(jié)點(diǎn)變形力的計(jì)算：

式中：[ B ]i、{Δ ε }i、{Δ σ}i、 wi、 Ji分別為高斯點(diǎn)i 的應(yīng)變矩陣、增量應(yīng)變向量、增量應(yīng)力向量、積分系數(shù)及雅克比行列式的值；{σn}i及{σo}i為高斯點(diǎn)i 當(dāng)前時(shí)刻及上一時(shí)刻的應(yīng)力全量；[ D ]、{Δ u }e、{ Fd}e分別為單元的彈性矩陣、節(jié)點(diǎn)增量位移向量及節(jié)點(diǎn)力向量；N為高斯點(diǎn)個(gè)數(shù)。

基于式（3），通過實(shí)時(shí)更新應(yīng)變矩陣[ B]i及節(jié)點(diǎn)坐標(biāo)[19]，可以實(shí)現(xiàn)有限元大位移（包含轉(zhuǎn)動(dòng)及平動(dòng)）、大變形問題的計(jì)算。

2.3 臨界條件

當(dāng)某單元滿足臨界條件時(shí)，該單元即被刪除，并創(chuàng)建顆粒簇。本文采用Mohr-Coulomb準(zhǔn)則及最大拉應(yīng)力準(zhǔn)則作為臨界條件，具體為

式中：σ1及σ3為最小及最大主應(yīng)力；c、φ、T為黏聚力、內(nèi)摩擦角及抗拉強(qiáng)度；Nφ、αp、σp為常數(shù)。

如果fs≤ 0且h≤ 0，則發(fā)生剪切破壞；如果ft≥ 0且h ＞ 0，則發(fā)生拉伸破壞。

2.4 顆粒生成

當(dāng)有限元單元滿足臨界條件被刪除后，需在有限元單元所在區(qū)域產(chǎn)生一定數(shù)量、隨機(jī)分布且微嵌入的顆粒系統(tǒng)。首先對(duì)顆粒半徑進(jìn)行隨機(jī)，隨機(jī)模式為均勻分布模式，設(shè)平均半徑為，則顆粒半徑的分布范圍為

式中：A為二維顆粒的面積；V為三維顆粒的體積；N為產(chǎn)生的顆粒總數(shù)。

接著對(duì)顆粒體心的坐標(biāo)進(jìn)行隨機(jī)，隨機(jī)模式亦為均勻分布模式。當(dāng)顆粒i 的體心坐標(biāo)位于有限元單元內(nèi)部，且顆粒i 到該單元所有面（棱）j 的距離及到所有已經(jīng)創(chuàng)建的顆粒k 的距離均滿足式（7）時(shí)，即創(chuàng)建顆粒i。

式中：Dij為顆粒i 到單元第j 個(gè)面（棱）的距離；dik為顆粒i 到第k 個(gè)顆粒的距離；Ri為顆粒i 的半徑；Rk為顆粒k 的半徑；tol為嵌入量（正值）。

顆粒i 創(chuàng)建后，其質(zhì)量、材料參數(shù)、速度、位移、接觸剛度、接觸力等信息均需從有限元單元中繼承。顆粒質(zhì)量由有限元單元的密度乘以顆粒體積或面積獲得；顆粒的黏聚力、內(nèi)摩擦角及抗拉強(qiáng)度值從有限元單元中完全繼承。

顆粒的速度、位移可通過單元節(jié)點(diǎn)插值獲得：

式中：vP、uP為顆粒速度、位移；Wj為單元第j個(gè)點(diǎn)的插值系數(shù)；為單元第j 個(gè)點(diǎn)的速度、位移；Ne為有限元單元的節(jié)點(diǎn)數(shù)。

假設(shè)顆粒間的接觸為有限面積的接觸（見圖2），則顆粒間的接觸剛度可根據(jù)顆粒的接觸面積及有限元單元的彈性模量、剪切模量獲得：

式中：Kn、Ks為顆粒間的法向及切向剛度；R1、R2為兩個(gè)接觸顆粒的半徑；E、G為有限元單元的彈性模量及剪切模量；Ac為顆粒間的接觸面積（可通過式（10）計(jì)算）。

圖2 顆粒間的有限接觸模型Fig.2 Finite contact model between particles

顆粒間的初始接觸力可通過插值獲得的位移進(jìn)行計(jì)算：

式中：Fn、Fs為顆粒間的法向及切向接觸力；Δun、Δ us為兩個(gè)接觸顆粒間的法向及切向位移差。

2.5 顆粒離散元的計(jì)算

本文在顆粒離散元計(jì)算時(shí)同時(shí)考慮力與力矩的傳遞作用，兩者的計(jì)算均基于增量法，且采用顯式求解方式。顆粒離散元間接觸力的計(jì)算式為

式中：Δdun、Δdus分別為兩個(gè)接觸顆粒間的法向及切向位移增量差。

接觸力矩的計(jì)算式為

式中：Mn、Ms為扭矩及彎矩；I 和J為接觸面的慣性矩及極慣性矩；Δdθn及Δdθs為顆粒間的扭轉(zhuǎn)及彎曲轉(zhuǎn)角增量差。

基于式（12）獲得當(dāng)前時(shí)步的試探接觸力后，需根據(jù)Mohr-Coulomb準(zhǔn)則及最大拉應(yīng)力準(zhǔn)則對(duì)接觸力進(jìn)行修正（法向力以壓為正）：

當(dāng)滿足不等式（16）中的任何一個(gè)式子時(shí)，顆粒間的接觸將不再傳遞力矩（即不進(jìn)行式（13）的計(jì)算）。其中

為了更加真實(shí)地模擬細(xì)觀顆粒的運(yùn)動(dòng)破壞過程，本文在計(jì)算時(shí)考慮了顆粒的轉(zhuǎn)動(dòng)效應(yīng)（計(jì)算示意圖如圖3所示。

圖3 顆粒的轉(zhuǎn)動(dòng)計(jì)算模型Fig.3 Particle-rolling computational model

計(jì)算時(shí)，首先利用式（17）獲得全局增量位移向量Δd，接著將Δd 轉(zhuǎn)換至接觸局部坐標(biāo)系，利用式（12）、（15）獲得局部坐標(biāo)系下的接觸力，進(jìn)而將接觸力轉(zhuǎn)化至整體坐標(biāo)系，并利用式（18）計(jì)算施加至顆粒1、2上的轉(zhuǎn)矩。

式中：w1及 w2為顆粒1、2的轉(zhuǎn)動(dòng)角速度向量；r1、r2為顆粒1、2到接觸點(diǎn)的相對(duì)位置向量（由顆粒質(zhì)心指向接觸點(diǎn)）；v1、v2為顆粒1、2質(zhì)心處的平動(dòng)速度向量。

式中：M1、M2為施加至顆粒1、2上的轉(zhuǎn)矩；F(G)為全局坐標(biāo)下的接觸力。

2.6 有限元與離散元的耦合

本文采用接觸模型進(jìn)行有限元與顆粒離散元的耦合。接觸模型包含接觸檢測(cè)及接觸力的計(jì)算兩個(gè)部分，接觸力的計(jì)算仍然采用增量法進(jìn)行，與顆粒-顆粒間的接觸力計(jì)算基本一致，在此處不再贅述。

接觸檢測(cè)包含初步及精確檢測(cè)兩個(gè)步驟。初步檢測(cè)采用子空間法實(shí)現(xiàn)，通過顆粒、單元與格子的映射關(guān)系縮小搜索范圍。精確檢測(cè)根據(jù)數(shù)值模型的維度可分為兩類：對(duì)于二維問題，實(shí)際上是顆粒與有限元邊界棱之間的接觸，采用點(diǎn)-棱接觸模型進(jìn)行描述；對(duì)于三維情況，實(shí)際上是有限元與邊界面之間的接觸，采用點(diǎn)-面接觸模型進(jìn)行描述。

點(diǎn)-棱接觸（見圖4）的創(chuàng)建需同時(shí)滿足兩個(gè)條件；一是顆粒質(zhì)心到邊界棱的距離（式（19））小于等于顆粒半徑(d ≤ R)；二是顆粒質(zhì)心在邊界棱上的投影點(diǎn)位于棱的內(nèi)部( dik≤ dij，djk≤ dij)。其中，Vpi為i 指向顆粒p 的相對(duì)位置向量，n為棱的單位外法向量。一旦顆粒與某邊界棱建立了接觸關(guān)系，法向彈簧及切向彈簧將自動(dòng)創(chuàng)建，接觸點(diǎn)k 的插值系數(shù)也將通過式（20）自動(dòng)求解。

圖4 點(diǎn)-棱接觸模型Fig.4 Point-edge contact model

點(diǎn)-面接觸（見圖5）的創(chuàng)建與點(diǎn)-棱接觸的創(chuàng)建過程基本一致。首先需判斷顆粒到某邊界面的距離是否小于顆粒半徑；接著將顆粒質(zhì)心點(diǎn)投影至邊界面上，判斷投影點(diǎn)是否位于某邊界面的內(nèi)部。如果上述兩個(gè)條件均滿足，則建立點(diǎn)-面接觸，進(jìn)而利用形函數(shù)計(jì)算接觸點(diǎn)的插值系數(shù)。

圖5 點(diǎn)-面接觸模型Fig.5 Point-face contact model

考慮到三維有限元的邊界面只有三角形及四邊形兩種形式，因此，可利用式（21）判斷某投影點(diǎn)是否位于三角形內(nèi)，利用式（22）判斷某投影點(diǎn)是否位于四邊形內(nèi)。

式中：A為面積；?為容差。

3 算 例

3.1 顆粒球與有限元厚板的碰撞分析

利用有限元轉(zhuǎn)化為顆粒離散元進(jìn)行巖體漸進(jìn)破壞過程的分析，有限元與顆粒之間的接觸計(jì)算不可避免。本算例通過顆粒球下落與有限元厚板的碰撞分析，驗(yàn)證顆粒與有限元接觸模型的正確性。數(shù)值模型如圖6所示。有限元厚板長(zhǎng)為10 m，寬為10 m，厚為0.5 m，由2 968個(gè)四面體單元組成；顆粒球直徑為6 m，球心位于厚板正上方5 m處，由1 277個(gè)顆粒組成，顆粒平均半徑為0.15 m。顆粒在重力作用下向Y 軸負(fù)方向運(yùn)動(dòng)，有限元板不考慮重力，并對(duì)板X=0 m及X=10 m處的節(jié)點(diǎn)進(jìn)行全約束。顆粒的彈性模量為0.1 GPa，泊松比為0.22，密度為2.5 g/cm3，內(nèi)摩擦角為20°，黏聚力及抗拉強(qiáng)度為0 MPa；有限元板的彈性模量為3 MPa，泊松比為0.25，密度為2.5 g/cm3。顆粒下落及與有限元厚板的接觸碰撞過程如圖7所示。由圖可得，隨著顆粒球與有限元板的接觸，顆粒球逐漸散開，有限元厚板出現(xiàn)較大變形；整個(gè)計(jì)算過程中顆粒與有限元之間未出現(xiàn)嵌入情況，且顆粒球與有限元厚板的運(yùn)動(dòng)與物理過程相符，證明了顆粒與有限元單元間接觸模型的正確性。

圖6 顆粒球與厚板碰撞模型Fig.6 Impact model for a particle ball on a slab

圖7 碰撞過程描述Fig.7 Description of impacting process

3.2 巖石單軸壓縮試驗(yàn)?zāi)M

建立10 cm×20 cm的巖石試樣，用468個(gè)三角形單元進(jìn)行剖分。模型底部全約束，模型頂部施加準(zhǔn)靜態(tài)豎直向下的速度載荷。采用傳統(tǒng)有限元及有限元轉(zhuǎn)化為顆粒離散元等兩種思路進(jìn)行巖石漸進(jìn)破壞過程的模擬。巖石彈性模量為30 GPa，泊松比為0.25，密度為2.5 g/cm3，抗拉強(qiáng)度及黏聚力均為3 MPa，內(nèi)摩擦角為30°。利用傳統(tǒng)有限元進(jìn)行模擬時(shí)，采用Mohr-Coulomb理想彈塑性模型及關(guān)聯(lián)流動(dòng)法則。

單軸壓縮過程中，上述兩種方法給出的巖石應(yīng)力-應(yīng)變的關(guān)系如圖8所示。由圖可得，在初期階段，兩種方法獲得的應(yīng)力-應(yīng)變曲線基本一致，軸向應(yīng)力均隨著軸向應(yīng)變的增加而線性增加。到達(dá)臨界狀態(tài)后，隨著軸向應(yīng)變的增加，傳統(tǒng)有限元方法給出的軸向應(yīng)力將保持不變，巖石出現(xiàn)理想彈塑性壓剪變形；而本文所述方法獲得的軸向應(yīng)力將迅速減小為0，巖石發(fā)生脆性壓剪破壞。由此，與傳統(tǒng)的有限元方法相比，本文所述方法可以更為真實(shí)地反映巖石單軸壓縮過程中的脆性破壞過程。

此外，本文所述方法獲得的峰值應(yīng)力約為10.15 MPa，與式（23）估算的臨界值（10.39 MPa）基本一致，表明了本文所述算法的準(zhǔn)確性。

圖8 巖石的應(yīng)力-應(yīng)變曲線Fig.8 Strain-stress curves of rock

利用有限元轉(zhuǎn)化為離散元方法計(jì)算獲得的巖石失穩(wěn)后的破裂形態(tài)如圖9所示。由圖可知，巖石在準(zhǔn)靜態(tài)豎向載荷作用下，發(fā)生了明顯的X型剪切破壞，破壞單元被一系列顆粒所填充，并由這些填充的顆粒完成了后續(xù)的破裂演化過程。

利用傳統(tǒng)有限元方法獲得的巖石失穩(wěn)后的水平方向位移及塑性剪應(yīng)變?cè)茍D如圖10所示。由圖可得，水平位移及塑性剪應(yīng)變均呈現(xiàn)出明顯的X型破壞特征，但材料并沒有發(fā)生斷裂，仍然處于連續(xù)變形狀態(tài)。

圖9 巖石的破裂形態(tài)Fig.9 Rock failure pattern

圖10 有限元計(jì)算獲得的水平位移及塑性剪應(yīng)變?cè)茍DFig.10 Contours of X-displacement and shear plastic strain based on the FEM method

3.3 巖石切割過程模擬

巖石尺寸為0.2 m×0.4 m，由112個(gè)有限元單元進(jìn)行離散；刀頭為三角形，三邊的長(zhǎng)度分別為6.8、10.8、11.1 cm。巖石底部全約束，刀頭以10 m/s的切割速度向巖石推進(jìn)。巖石的彈性模量為30 GPa，泊松比為0.25，密度為2.5 g/cm3，黏聚力為5 MPa，內(nèi)摩擦角為40°，抗拉強(qiáng)度為2 MPa。切割過程中巖石的漸進(jìn)破壞過程如圖11所示。巖石頂部左右兩側(cè)節(jié)點(diǎn)P1及P2（具體位置見圖11）的水平位移隨時(shí)間的變化曲線如圖12所示。

圖11 巖石切割破裂過程Fig.11 Rock failing process during cutting

圖12 節(jié)點(diǎn)P1及P2的水平位移時(shí)程曲線Fig.12 X-displacement vs.time at points P1and P2

由圖11可知，隨著刀片的推進(jìn)，刀片附近的巖石逐漸破碎，轉(zhuǎn)化為顆粒流；巖石中宏觀裂縫的擴(kuò)展過程通過顆粒之間的細(xì)觀接觸力演化自然獲得。由圖12可知，隨著時(shí)間的推移，巖石頂部左、右兩側(cè)節(jié)點(diǎn)P1、P2的水平位移逐漸增大；當(dāng)4.8 ms以后，位移的增加速率明顯增大，表明此時(shí)巖石已經(jīng)出現(xiàn)了失穩(wěn)破壞。圖11、12展示的巖石破裂過程及位移時(shí)程曲線與實(shí)際較為相符，證明了本文所述方法的合理性。

4 結(jié) 論

（1）本文所述有限元與顆粒離散元的耦合轉(zhuǎn)化模型，結(jié)合了有限元與顆粒離散元各自的優(yōu)勢(shì)，較為適合模擬地質(zhì)體在各類靜、動(dòng)載荷作用下的漸進(jìn)破壞過程。

（2）在有限元與顆粒離散元的耦合轉(zhuǎn)化模型中，地質(zhì)體初期的宏觀應(yīng)力-應(yīng)變由有限元計(jì)算，地質(zhì)體后期的漸進(jìn)破壞過程由顆粒離散元模擬；并利用Mohr-Coulomb準(zhǔn)則及最大拉應(yīng)力準(zhǔn)則進(jìn)行轉(zhuǎn)換點(diǎn)的判斷。

（3）顆粒球與有限元板的碰撞分析、巖石單軸壓縮過程模擬、巖石切割過程模擬等案例，證明了有限元與顆粒離散元接觸模型、轉(zhuǎn)化模型的合理性及正確性。

（4）有限元與顆粒離散元耦合轉(zhuǎn)化模型的后續(xù)研究重點(diǎn)將包括如下幾個(gè)方面：轉(zhuǎn)化點(diǎn)所用宏觀準(zhǔn)則的優(yōu)選、隨機(jī)顆粒系統(tǒng)生成方案的優(yōu)選、顆粒不同信息繼承方式的對(duì)比、本文所述方法的進(jìn)一步校核（與試驗(yàn)及其他數(shù)值模擬方法的對(duì)比），利用本文方法進(jìn)行三維巖體漸進(jìn)破壞過程的模擬等。

[1]李世海，劉天蘋，劉曉宇.論滑坡穩(wěn)定性分析方法[J].巖石力學(xué)與工程學(xué)報(bào)，2009，28(增刊2):3309－3324.LI Shi-hai，LIU Tian-ping，LIU Xiao-yu.Analysis method for landslide stability[J].Chinese Journal of Rock Mechanics and Engineering，2009，28(Supp.2):3309－3324.

[2]MIEHE C，GüRSES E.A robust algorithm for configurational-force-driven brittle crack propagation with R-adaptive mesh alignment[J].International Journal for Numerical Methods in Engineering，2007，72(2):127－155.

[3]KHOEI A R，AZADI H，MOSLEMI H.Modeling of crack propagation via an automatic adaptive mesh refinement based on modified superconvergent patch recovery technique[J].Engineering Fracture Mechanics，2008，75(10):2921－2945.

[4]DOLBOW J，BELYTSCHKO T.A finite element method for crack growth without remeshing[J].International Journal for Numerical Methods in Engineering，1999，46:131－150.

[5]ASADPOURE A，MOHAMMADI S.Developing new enrichment functions for crack simulation in orthotropic media by the extended finite element method[J].International Journal for Numerical Methods in Engineering，2007，69(10):2150－2172.

[6]GORDELIY E，PEIRCE A.Implicit level set schemes for modeling hydraulic fractures using the XFEM[J].Computer Methods in Applied Mechanics and Engineering，2013，266:125－143.

[7]JIANG Y，LI B，YAMASHITA Y.Simulation of cracking near a large underground cavern in a discontinuous rock mass using the expanded distinct element method[J].International Journal of Rock Mechanics and Mining Sciences，2009，46(1):97－106.

[8]AMARASIRI A L，KODIKARA J K.Numerical modelling of desiccation cracking using the cohesive crack method[J].International Journal of Geomechanics，2011，13(3):213－221.

[9]BARLA M，PIOVANO G，GRASSELLI G.Rock slide simulation with the combined finite-discrete element method[J].International Journal of Geomechanics，2011，12(6):711－721.

[10]LI SH，ZHAO MH，WANG YN，et al.A new numerical method for DEM-block and particle model[J].International Journal of Rock Mechanics and Mining Sciences，2004，41(3):436－436.

[11]COTTRELL M G，YU J，WEI Z J，et al.The numerical modelling of ceramics subject to impact using adaptive discrete element techniques[J].Engineering Computations，2003，20(1):82－106.

[12]王杰，李世海，周東，等.模擬巖石破裂過程的塊體單元離散彈簧模型[J].巖土力學(xué)，2013，34(8):2355－2362.WANG Jie，LI Shi-hai，ZHOU Dong，et al.A blockdiscrete-spring model to simulate failure process of rock[J].Rock and Soil Mechanics，2013，34(8):2355－2362.

[13]VESGA L F，VALLEJO L E，LOBO GUERRERO S.DEM analysis of the crack propagation in brittle clays under uniaxial compression tests[J].International Journal for Numerical and Analytical Methods in Geomechanics，2008，32(11):1405－1415.

[14]YANG D，SHENG Y，YE J，et al.Dynamic simulation of crack initiation and propagation in cross-ply laminates by DEM[J].Composites Science and Technology，2011，71(11):1410－1418.

[15]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.

[16]YANG B，JIAO Y，LEI S.A study on the effects of microparameters on macroproperties for specimens created by bonded particles[J].Engineering Computations，2006，23(6):607－631.

[17]O?ATE E，ROJEK J.Combination of discrete element and finite element methods for dynamic analysis of geomechanics problems[J].Computer Methods in Applied Mechanics and Engineering，2004，193(27):3087－3128.

[18]ROJEK J，O?ATE E.Multiscale analysis using a coupled discrete/finite element model[J].Interaction and Multiscale Mechanics，2007，1(1):1－31.

[19]FENG C，LI SH，LIU XY，et al.A semi-spring and semi-edge combined contact model in CDEM and its application to analysis of Jiweishan landslide[J].Journal of Rock Mechanics and Geotechnical Engineering，2014，6(1):26－35.