金愛兵，王 凱，張秀鳳，孟新秋，楊振偉

（1.北京科技大學(xué) 土木與環(huán)境工程學(xué)院，北京 100083；2.溫州市城建設(shè)計院，浙江 溫州 325000）

1 引 言

土體作為工程材料的一種，在外荷載作用下的力學(xué)響應(yīng)不能完全由經(jīng)典力學(xué)描述。特別是軟土，由于其內(nèi)部結(jié)構(gòu)特殊，其應(yīng)力-應(yīng)變關(guān)系呈現(xiàn)出一定的非線性流變特征。陳宗基[1]通過將Maxwell 體與Hook 彈性體并聯(lián)，得到一種新模型來描述土體變形規(guī)律；孫鈞[2]對上海地區(qū)3 種典型軟土進(jìn)行了大量的單軸剪切與三軸壓縮流變試驗，并對試驗數(shù)據(jù)進(jìn)行詳細(xì)分析后，提出了上海軟土蠕變特性的方程式。粉質(zhì)黏土作為一種常見軟土，也常常表現(xiàn)出一定的流變特性。王坤等[3]對深圳地區(qū)深層粉質(zhì)黏土進(jìn)行了低圍壓作用下的三軸流變試驗，并得到了粉質(zhì)黏土的三軸流變規(guī)律。而對于服役時間較長的露天礦山邊坡，粉質(zhì)黏土的流變效應(yīng)對邊坡工程的安全影響更為顯著。

目前，對巖土流變特性的研究主要集中在室內(nèi)試驗[4－6]和流變模型開發(fā)[7－9]等方面。室內(nèi)試驗研究主要有單軸試驗[6]和不同圍壓下的三軸試驗[10－14]。熊良宵等[15]將修正后的Burgers 模型應(yīng)用于蠕變試驗的擬合分析；李娜等[16]用改進(jìn)的西原模型對循環(huán)加、卸載試驗結(jié)果進(jìn)行擬合，所得效果良好。但室內(nèi)試驗研究往往存在操作復(fù)雜、花費(fèi)時間長、試驗費(fèi)用高等缺點。

作為巖土工程穩(wěn)定性研究的常用工具，數(shù)值模擬方法已經(jīng)被許多科研工作者應(yīng)用于巖土流變特性的研究。張志沛等[17]利用FLAC3D中的Burgers 模型對不同形狀的試樣進(jìn)行了蠕變特性模擬分析；熊良宵等[18]在FLAC3D中開發(fā)出新的流變模型，并將其應(yīng)用到巖土應(yīng)變軟化的相關(guān)工程計算中；劉珊珊等[19]在FLAC3D中二次開發(fā)出黏彈性廣義Kelvin 模型，并用算例驗證了所開發(fā)模型的正確性；王濤等[20]開發(fā)出適用于PFC2D的廣義Kelvin 接觸本構(gòu)模型，并將其應(yīng)用到了隧洞工程問題中。

本文在對PFC2D中自定義本構(gòu)模型進(jìn)行研究的基礎(chǔ)上，利用Microsoft Visual Studio 平臺，開發(fā)出廣義Kelvin 模型的動態(tài)鏈接庫(DLL)文件。通過算例驗證程序編制的準(zhǔn)確性，并在此基礎(chǔ)上，將其應(yīng)用到露天煤礦含粉質(zhì)黏土邊坡的失穩(wěn)研究中。

2 廣義Kelvin 模型本構(gòu)關(guān)系

廣義Kelvin 模型由一個Maxwell 體和任意個Kelvin 體串聯(lián)組成，其模型如圖1 所示。其中，Maxwell 體由一個彈簧和一個阻尼器串聯(lián)而成，Kelvin 體由一個彈簧和一個阻尼器并聯(lián)而成，Kn和Ks代表彈簧法向和切向剛度，Cn和Cs代表阻尼器法向和切向阻尼，下標(biāo)數(shù)字1,2…,n為原件編號。

圖1 廣義Kelvin 模型Fig.1 Generalized Kelvin model

由圖1 可知，廣義Kelvin 模型的總位移u 等于一個Maxwell 體的位移 um1與若干個Kelvin 體的位移ukr（下標(biāo)r為開爾文體的編號，r=2,3,…,n）之和，如式（1）所示。其中，um1由Maxwell 體中彈簧的位移 uk1和阻尼器的位移 uc1共同組成。

對廣義Kelvin 模型的本構(gòu)關(guān)系推導(dǎo)如下，下列各式中的“±”(或“?”)分別對應(yīng)著法向和切向的情況；Kr和Cr(r=1,2,…,n) 分別為廣義Kelvin 模型中彈簧的剛度和阻尼器的阻尼；f為模型中顆粒間的接觸力；Ar、Br、C和D 是分別與Kr、Cr和Δt 有關(guān)的中間變量。

對式（1）進(jìn)行一階求導(dǎo)，得

接觸力f 可表示為

將式（3）、（5）、（6）代入式（2）得

利用有限差分格式的中心差分法近似對時間t進(jìn)行差分，其中對t 取平均值得

對于Kelvin 體，則有

整理后，接觸力ft+1可表示為

ft+1可由已知量ut+1、ut、和ft求出。

3 廣義Kelvin 模型開發(fā)

PFC 程序中的用戶自定義接觸模型的編寫包括基類描述、模型注冊、成員函數(shù)描述、模型和PFC程序間的信息交換。軟件還提供了自帶本構(gòu)模型的源代碼，用戶在進(jìn)行本構(gòu)模型開發(fā)時，可以利用和修改這些代碼。下面介紹如何修改本構(gòu)模型的頭文件和源文件，從而創(chuàng)建用戶自己的本構(gòu)模型。

在頭文件中進(jìn)行新的本構(gòu)模型派生類的聲明，修改模型的編號、名稱和版本，修改派生類的私有成員，包括模型的基本變量及主要中間變量。為避免沖突，新自定義模型的編號需不小于100。

在源文件中修改以下各個成員函數(shù)：

（1）修改構(gòu)造函數(shù)Kelvin::Kelvin()，對其私有成員進(jìn)行初始化。

（2）修改Kelvin::Name()函數(shù)，返回用戶自定義模型的名稱。

（3）根據(jù)派生類私有成員的名稱字符串，依次修改Kelvin::PropNames()、Kelvin::ReturnProp()和Kelvin::AcceptProp()3 個函數(shù)。PFC 的PROPERTY命令可用來調(diào)用Kelvin::AcceptProp()函數(shù)，對模型的參數(shù)進(jìn)行賦值。

（4）修改Kelvin::FDlaw()函數(shù)，此函數(shù)在計算離散元的力與位移關(guān)系時被調(diào)用。用戶根據(jù)推導(dǎo)的模型中力與位移數(shù)學(xué)關(guān)系，編寫C++語句，由應(yīng)變增量計算得到應(yīng)力增量，從而獲得新的應(yīng)力。本構(gòu)模型將在每個循環(huán)中更新顆粒間接觸的力與力矩，以及法向和切向剛度。

（5）修改USERCM1::SaveRestore()函數(shù)，方法同步驟（3），此函數(shù)用來保存與恢復(fù)實型、整型和邏輯型的內(nèi)部數(shù)據(jù)，PFC 程序給出save和restore命令時調(diào)用此函數(shù)。

另外，針對本文自定義的廣義Kelvin本構(gòu)模型，程序中設(shè)置控制變量i來靈活控制Kelvin體的個數(shù)。編寫Kelvin.h和Kelvin.cpp 文件后，編譯動態(tài)鏈接庫文件。至此，所開發(fā)的廣義Kelvin 模型以動態(tài)鏈接庫文件的形式存儲起來，可供PFC2D程序調(diào)用。

4 廣義Kelvin 模型驗證

當(dāng)控制變量i=0 時，即廣義Kelvin 模型中只含有Maxwell 體，無Kelvin 體時，模型即為Maxwell模型；當(dāng)控制變量i=1 時，即廣義Kelvin 模型中只含有Maxwell 體和1 個Kelvin 體時，模型即為Burgers 模型。本文分別采用i=0和1 兩種情況下的Maxwell 模型和Burgers 模型進(jìn)行松弛試驗的驗證。

（1）Maxwell 模型理論值求解

當(dāng)i=0，無Kelvin 體時，式（7）可表示為

（2）Burgers 模型理論值求解

當(dāng)i=1 時，對于單位輸入函數(shù) u( t)=1，Burgers模型的理論值為

兩種模型中Maxwell 體彈簧法向、切向剛度：K1n=K1s=100MPa；阻尼器法向、切向阻尼：C1n=C1s=100MPa·s；Kelvin體彈簧法向、切向剛度：K2n=K2s=100MPa，阻尼器法向、切向阻尼：C2n=C2s=100MPa·s。

通過上述計算，得到兩種接觸模型松弛過程的理論值后，采用PFC2D命令流和Fish 語言編寫應(yīng)力松弛程序，然后在PFC2D軟件中調(diào)用該程序，分別加載i為0和1 時的廣義Kelvin 本構(gòu)模型動態(tài)鏈接庫文件，模擬兩個相互接觸球體（兩球重疊量為0.01 m）的應(yīng)力松弛過程，得出球體間法向接觸力隨時間的變化關(guān)系，并將數(shù)值模擬解與理論值進(jìn)行對比分析。其對比結(jié)果如圖2 所示。

由圖2 可知，兩球接觸力隨時間的增加呈指數(shù)趨勢減小，本構(gòu)模型的數(shù)值模擬解與理論解所對應(yīng)曲線相互吻合，從而驗證了所開發(fā)模型的準(zhǔn)確性。

圖2 法向接觸力模擬計算結(jié)果與理論結(jié)果Fig.2 Analytical and numerical results of normal contact forces

5 工程應(yīng)用

5.1 邊坡失穩(wěn)規(guī)律

含軟弱粉質(zhì)黏土邊坡廣泛分布于山西、陜西等黃土高原地區(qū)，尤其是在山西地區(qū)，煤炭資源豐富，露天開采形成的含粉質(zhì)黏土邊坡普遍存在。在山西某露天煤礦，經(jīng)長期實地調(diào)研，發(fā)現(xiàn)此類邊坡的滑坡形式較為特殊，其失穩(wěn)規(guī)律有以下幾個特征：①邊坡的滑坡區(qū)域含有粉質(zhì)黏土層；②巖土接觸帶傾向坡內(nèi)，傾角較小(6°～9°)，且有地下水滲出；③滑坡體前緣平緩，局部有底鼓隆起現(xiàn)象；④滑坡體中部較為平緩；⑤滑坡體后緣較陡，張裂縫嚴(yán)重，下落量小；⑥ 滑動面形式表現(xiàn)異常，為“圓弧+折線”的組合型，最終輪廓表現(xiàn)為帽檐形。

此類含粉質(zhì)黏土邊坡失穩(wěn)破壞的一般過程與現(xiàn)場情況如圖3 所示。

5.2 模型計算參數(shù)

巖土體的物理力學(xué)參數(shù)由巖石壓縮試驗、土工直剪試驗及粉質(zhì)黏土三軸流變試驗獲得，見表1。參照室內(nèi)試驗，進(jìn)行數(shù)值模擬試驗，粉質(zhì)黏土部分調(diào)用本文自主開發(fā)的廣義Kelvin 模型，其他巖土體均采用PFC2D軟件默認(rèn)線彈性接觸模型；數(shù)值模型尺寸均為100 mm×50 mm 的標(biāo)準(zhǔn)試樣；巖土體顆粒間均采用能夠抗拉、抗剪、承受彎矩的平行黏結(jié)。

模擬試驗過程中，對廣義Kelvin 模型中的控制變量i 及模型的細(xì)觀力學(xué)參數(shù)進(jìn)行反復(fù)調(diào)試。最終結(jié)果表明，當(dāng)控制變量i=2，即廣義Kelvin 模型中含1 個Maxwell 體和2 個Kelvin 體時，粉質(zhì)黏土的模擬試驗與室內(nèi)試驗結(jié)果較吻合（見圖4、5）。此時，在廣義Kelvin 模型中，Maxwell 體彈簧法向、切向剛度：K1n=K1s=100 MPa；阻尼器法向、切向阻尼：C1n=C1s=300 MPa·s；兩個Kelvin 體相同，彈簧法向、切向剛度：K2n=K2s=K3n=K3s=100 MPa；阻尼器法向、切向阻尼：C2n=C2s=C3n=C3s=300 MPa·s 。邊坡模型中巖土體的其他細(xì)觀力學(xué)參數(shù)見表2。

圖3 邊坡失穩(wěn)破壞過程與現(xiàn)場Fig.3 Process and scene of slope failure

表1 巖土體力學(xué)參數(shù)Table 1 Mechanical parameters of rock and soils

圖4 數(shù)值試驗與直剪試驗結(jié)果Fig.4 Results of numerical test and direct shear test

圖5 數(shù)值試驗與流變試驗結(jié)果Fig.5 Results of numerical test and rheological test

5.3 邊坡模型構(gòu)建

根據(jù)參數(shù)調(diào)試過程中獲取的廣義Kelvin模型和細(xì)觀力學(xué)參數(shù)，建立邊坡開挖前的PFC2D數(shù)值計算模型，見圖6。根據(jù)工程實際情況，邊坡模型長×高尺寸為：217.51 m×82.66 m。模型內(nèi)部共生成36 261 個巖土體顆粒，模型的外部由1 道水平墻體和2 道豎直墻體構(gòu)成邊界限制條件，墻體法向剛度取 kn=1 000 MN/m，切向剛度取 ks=1 000 MN/m，摩擦系數(shù)取1.0。模型中不同顏色顆粒分別代表不同地層，邊坡自上而下依次為粉土、粉質(zhì)黏土（天然）、粉質(zhì)黏土（飽和）、接觸帶及泥巖。

表2 顆粒細(xì)觀力學(xué)參數(shù)Table 2 Meso mechanical parameters of particles

圖6 邊坡數(shù)值模型Fig.6 Numerical model of slope

對于粉質(zhì)黏土，考慮到吸水-失水循環(huán)對其強(qiáng)度的影響，根據(jù)現(xiàn)場實際情況，在邊坡模型的不同地層位置和計算的不同時段，對粉質(zhì)黏土分別采用天然強(qiáng)度與飽和強(qiáng)度。無降水時下層粉質(zhì)黏土采用天然強(qiáng)度，降水時下層粉質(zhì)黏土采用飽和強(qiáng)度。后續(xù)數(shù)值模擬中，邊坡開挖后通過在計算的不用時步對該層粉質(zhì)黏土賦予相應(yīng)參數(shù)來實現(xiàn)這一過程，根據(jù)當(dāng)?shù)赜昙镜慕邓?guī)律：模擬中每40 萬時步采用天然強(qiáng)度，每10 萬時步采用飽和強(qiáng)度，循環(huán)多次至邊坡最終破壞形成，共計約240 萬時步。

5.4 計算結(jié)果及分析

在PFC2D中編寫命令流，實現(xiàn)對邊坡數(shù)值模型的計算與開挖。對于邊坡中的粉質(zhì)黏土地層，分別采用線彈性接觸模型和本文所開發(fā)的具有流變特性的廣義Kelvin 模型兩種情況進(jìn)行模擬計算，兩算例中吸水-失水循環(huán)模擬過程相同。最終將其計算結(jié)果進(jìn)行對比分析。

計算中采用顆粒間平行黏結(jié)的拉裂或剪斷來表示邊坡巖土體結(jié)構(gòu)的破壞[21]，為便于觀察，將模型中黏結(jié)破壞的顆粒設(shè)置為深藍(lán)色。模型中粉質(zhì)黏土地層采用線彈性接觸模型的計算結(jié)果如圖7 所示，邊坡最終破壞模式呈常見的圓弧形土體滑坡模式，但與現(xiàn)場邊坡實際破壞情況不符。

圖7 線彈性接觸模型邊坡破壞結(jié)果Fig.7 Slope failure predicted by a linear elastic contact model

模型中粉質(zhì)黏土地層采用用戶自定義廣義Kelvin 模型的計算結(jié)果，如圖8 所示，邊坡最終破壞模式比較符合現(xiàn)場的實際情況。計算過程中，邊坡的失穩(wěn)破壞可分為5 個階段。

第1階段：邊坡施工開挖，如圖8(a)所示，各臺階開挖臨空面小范圍內(nèi)巖土體受到擾動，深藍(lán)色顆粒開始出現(xiàn)，并零星散布于邊坡面的淺層。

第2階段：邊坡局部破壞，如圖8(b)所示，邊坡內(nèi)部整體變化微小，因賦存地下水，粉質(zhì)黏土層強(qiáng)度降低，其坡腳處開始有深藍(lán)色顆粒聚集，局部位移增大，破壞初顯。

第3階段：破壞范圍擴(kuò)大，如圖8(c)所示，飽和粉質(zhì)黏土層邊坡坡腳處深藍(lán)色顆粒增多，孔隙率增加，破壞范圍擴(kuò)大，少量黏結(jié)已被破壞的顆粒在重力作用下開始滑移出坡面，邊坡后緣開始呈現(xiàn)出較明顯的圓弧滑面。

第4階段：邊坡出現(xiàn)失穩(wěn)，如圖8(d)所示，坡面周邊大量深藍(lán)色顆粒翻滾下滑，孔隙率大幅增加，邊坡破壞區(qū)域輪廓成“折線+圓弧”的帽檐形。

第5階段：帽檐形滑坡形成，如圖8(e)所示，模型中的顆粒停止下滑，邊坡模型已經(jīng)破壞。

圖8 廣義Kelvin 模型邊坡破壞過程Fig.8 Slope failure process derived from the generalized Kelvin model

通過上述計算，并結(jié)合現(xiàn)場滑坡情況，邊坡失穩(wěn)過程主要表現(xiàn)為：

（1）天然狀態(tài)下，粉質(zhì)黏土層強(qiáng)度較高，所在邊坡具有足夠的自穩(wěn)能力，隨著開采的進(jìn)行，粉質(zhì)黏土層中顆粒黏結(jié)遭到破壞，導(dǎo)致裂隙產(chǎn)生，計算過程中坡體共產(chǎn)生498 條裂隙。

（2）微觀裂隙貫通產(chǎn)生宏觀裂縫，受水影響，飽和粉質(zhì)黏土在接觸帶處強(qiáng)度較低，顆粒黏結(jié)容易遭到破壞，所以該處產(chǎn)生裂隙最多，成為邊坡破壞的主要區(qū)域。

（3）受開采影響，失水-吸水循環(huán)過程中，裂隙進(jìn)一步擴(kuò)張，松散體影響范圍逐漸擴(kuò)大，并在邊坡后緣開始形成圓弧滑面。

（4）因松散體自身承載力不夠，內(nèi)部土體顆粒開始不斷下滑，在粉質(zhì)黏土流變作用下，一段時間后，粉質(zhì)黏土不能承受上部土體的荷載，松散體呈折線形式向上蔓延，聯(lián)通至邊坡后緣圓弧滑面，以“折線+圓弧”的組合滑動面形式破壞，破壞區(qū)最終輪廓呈帽檐形。

含粉質(zhì)黏土邊坡失穩(wěn)模擬結(jié)果表明，該類邊坡失穩(wěn)與普通土質(zhì)邊坡的失穩(wěn)形態(tài)以及失穩(wěn)過程并不一致，粉質(zhì)黏土的流變性對該類邊坡失穩(wěn)具有重要影響；而水是在粉質(zhì)黏土流變過程中導(dǎo)致邊坡失穩(wěn)的重要原因，因此，實際工程中可以采取壓坡阻水的方式處治該類邊坡，以保證邊坡穩(wěn)定。

6 結(jié) 論

（1）自定義廣義Kelvin 本構(gòu)模型包含1 個Maxwell 體和若干個Kelvin 體，實際應(yīng)用中可以根據(jù)具體情況選擇合適數(shù)量的Kelvin 體，以便更好地模擬巖土體流變特性。

（2）在PFC2D軟件中調(diào)用基于C++語言開發(fā)的自定義廣義Kelvin 本構(gòu)模型動態(tài)鏈接庫，進(jìn)行算例驗證，得到的數(shù)值模擬解與理論解吻合，模型具有較好的準(zhǔn)確性。

（3）利用PFC 調(diào)用自定義廣義Kelvin 本構(gòu)模型對國內(nèi)某含粉質(zhì)黏土邊坡模擬表明，含1個Maxwell體和2 個Kelvin 體的用戶自定義廣義Kelvin 模型可很好地模擬粉質(zhì)黏土的流變特性，得到的邊坡失穩(wěn)破壞過程與現(xiàn)場基本一致。

（4）由于流變特性的影響，該類粉質(zhì)黏土邊坡破壞形式較為特殊，與常規(guī)土質(zhì)邊坡破壞呈圓弧滑動不一致，其破壞區(qū)最終輪廓成“折線+圓弧”的帽檐形。

[1]陳宗基.固結(jié)及時間效應(yīng)的單維問題[J].土木工程學(xué)報,1958,5(1):1－10.CHEN Zong-ji.Unidimensional problems of consolidation and time effect[J].China Civil Engineering Journal,1958,5(1):1－10.

[2]孫鈞.巖土材料流變及其工程應(yīng)用[M].北京:中國建筑工業(yè)出版社,1999.

[3]王坤,蔡永昌,魏赟,等.深圳地區(qū)深層粉質(zhì)黏土的流變特性試驗研究[J].巖土工程界,2009,12(12):86－89.WANG Kun,CAI Yong-chang,WEI Yun,et al.Experimental study on the rheological properties of deep silty clay in Shenzhen[J].Geotechnical Engineering World,2009,12(12):86－89.

[4]ZHANG ZHILIANG,XU WEIYA,WANG WEI,et al.Triaxial creep tests of rock from the compressive zone of dam foundation in Xiangjiaba Hydropower Station[J].International Journal of Rock Mechanics and Mining Sciences,2012,50(1):133－139.

[5]胡華,鄭曉栩.動載作用頻率對海相沉積軟土動態(tài)流變特性影響試驗研究[J].巖土力學(xué),2013,34(增刊1):9－13.HU Hua,ZHENG Xiao-xu.Experimental research on dynamic rheological characteristics of marine deposit soft soil under different frequencies of dynamic loading[J].Rock and Soil Mechanics,2013,34(Supp.1):9－13.

[6]LI YONGSHENG,XIA CAICHU.Time-dependent tests on intact rocks in uniaxial compression[J].International Journal of Rock Mechanics and Mining Sciences,2000,37(3):467－475.

[7]ZHANG HUABIN,WANG ZHIYIN,ZHENG YALI,et al.Study on triaxial creep experiment and constitutive relation of different rock salts[J].Safety Science,2012,50(4):801－805.

[8]LIU LANG,WANG GEMIN,CHEN JIANHONG,et al.Creep experiment and rheological model of deep saturated rock[J].Transactions of Nonferrous Metals Society of China,2013,23(2):478－483.

[9]SHAO J F,ZHU Q Z,SU K.Modeling of creep in rock materials in terms of material degradation[J].Computers and Geotechnics,2003,30(7):549－555.

[10]ZHANG YU,XU WEIYA,GU JINJIAN,et al.Triaxial creep tests of weak sandstone from fracture zone of high dam foundation[J].Journal of Central South University,2013,20(9):2528－2536.

[11]YANG SHENGQI,JIANG YUZHOU.Triaxial mechanical creep behavior of sandstone[J].Mining Science and Technology,2010,20(3):339－349.

[12]ZHAO YANLIN,CAO PING,WANG WEIJUN,et al.Viscoelasto-plastic rheological experiment under circular increment step load and unload and nonlinear creep model of soft rocks[J].Journal of Central South University of Technology,2009,16(3):484－491.

[13]李亞麗,于懷昌,劉漢東.三軸壓縮下粉砂質(zhì)泥巖蠕變本構(gòu)模型研究[J].巖土力學(xué),2012,33(7):2035－2040.LI Ya-li,YU Huai-chang,LIU Han-dong.Study of creep constitutive model of silty mudstone under triaxial compression[J].Rock and Soil Mechanics,2012,33(7):2035－2040.

[14]徐進(jìn),張家生,趙同順,等.軟弱路基土體三軸蠕變試驗及蠕變模型研究[J].中南大學(xué)學(xué)報(自然科學(xué)版),2011,42(10):3136－3142.XU Jin,ZHANG Jia-sheng,ZHAO Tong-shun,et al.Creep triaxial tests and constitutive model of embankment foundation soil[J].Journal of Central South University(Science and Technology),2011,42(10):3136－3142.

[15]熊良宵,楊林德,張堯.巖石的非定常Burgers 模型[J].中南大學(xué)學(xué)報(自然科學(xué)版),2010,41(2):679－684.XIONG Liang-xiao,YANG Lin-de,ZHANG Yao.Non-stationary Burgers model for rock[J].Journal of Central South University(Science and Technology),2010,41(2):679－684.

[16]李娜,曹平,衣永亮,等.分級加卸載下深部巖石流變實驗及模型[J].中南大學(xué)學(xué)報(自然科學(xué)版),2011,42(11):3465－3471.LI Na,CAO Ping,YI Yong-liang,et al.Creep properties experiment and model of deep rock with step loading and unloading[J].Journal of Central South University(Science and Technology),2011,42(11):3465－3471.

[17]張志沛,王芝銀,彭惠.陜南泥巖三軸壓縮蠕變試驗及其數(shù)值模擬研究[J].水文地質(zhì)工程地質(zhì),2011,38(1):53－58.ZHANG Zhi-pei,WANG Zhi-yin,PENG Hui.Study on triaxial compression creep test and numerical simulation about mud rock in southern Shaanxi[J].Hydrogeology &Engineering Geology,2011,38(1):53－58.

[18]熊良宵,楊林德,張堯.硬巖的復(fù)合黏彈塑性流變模型[J].中南大學(xué)學(xué)報(自然科學(xué)版),2010 41(4):1540－1548.XIONG Liang-xiao,YANG Lin-de,ZHANG Yao.Composite viscoelasto-plastic rheological model for hard rock[J].Journal of Central South University(Science and Technology),2010,41(4):1540－1548.

[19]劉姍姍,趙同彬.黏彈性廣義Kelvin 模型的FLAC3D二次開發(fā)[J].山東科技大學(xué)學(xué)報(自然科學(xué)版),2010,29(4):20－23.LIU Shan-shan,ZHAO Tong-bin.Secondary development on generalized viscoelastic Kelvin model with FLAC3D[J].Journal of Shandong University of Science and Technology(Natural Science),2010,29(4):20－23.

[20]王濤,呂慶,李楊,等.顆粒離散元方法中接觸模型的開發(fā)[J].巖石力學(xué)與工程學(xué)報,2009,28(增刊2):4040－4045.WANG Tao,Lü Qing,LI Yang,et al.Development of contact model in particle discrete element method[J].Chinese Journal of Rock Mechanics and Engineering,2009,28(Supp.2):4040－4045.

[21]張曉平,吳順川,王思敬.類土質(zhì)路塹邊坡動態(tài)監(jiān)測及數(shù)值模擬分析[J].巖石力學(xué)與工程學(xué)報,2008,27(增刊2):3431－3440.ZHANG Xiao-ping,WU Shun-chuan,WANG Si-jing.Dynamic monitoring and numerical analysis of soil-like cut slope[J].Chinese Journal of Rock Mechanics andEngineering,2008,27(Supp.2):3431－3440.