基于顆粒流程序的廣義Kelvin 模型及其應用

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

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

1 引 言

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

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

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

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

2 廣義Kelvin 模型本構關系

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

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

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

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

對式（1）進行一階求導，得

接觸力f 可表示為

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

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

對于Kelvin 體，則有

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

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

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

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

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

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

（1）修改構造函數Kelvin::Kelvin()，對其私有成員進行初始化。

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

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

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

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

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

4 廣義Kelvin 模型驗證

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

（1）Maxwell 模型理論值求解

當i=0，無Kelvin 體時，式（7）可表示為

（2）Burgers 模型理論值求解

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

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

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

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

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

5 工程應用

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

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

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

5.2 模型計算參數

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

模擬試驗過程中，對廣義Kelvin 模型中的控制變量i 及模型的細觀力學參數進行反復調試。最終結果表明，當控制變量i=2，即廣義Kelvin 模型中含1 個Maxwell 體和2 個Kelvin 體時，粉質黏土的模擬試驗與室內試驗結果較吻合（見圖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 。邊坡模型中巖土體的其他細觀力學參數見表2。

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

表1 巖土體力學參數Table 1 Mechanical parameters of rock and soils

圖4 數值試驗與直剪試驗結果Fig.4 Results of numerical test and direct shear test

圖5 數值試驗與流變試驗結果Fig.5 Results of numerical test and rheological test

5.3 邊坡模型構建

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

表2 顆粒細觀力學參數Table 2 Meso mechanical parameters of particles

圖6 邊坡數值模型Fig.6 Numerical model of slope

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

5.4 計算結果及分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

6 結 論

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

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

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

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

[1]陳宗基.固結及時間效應的單維問題[J].土木工程學報,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]孫鈞.巖土材料流變及其工程應用[M].北京:中國建筑工業(yè)出版社,1999.

[3]王坤,蔡永昌,魏赟,等.深圳地區(qū)深層粉質黏土的流變特性試驗研究[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].巖土力學,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]李亞麗,于懷昌,劉漢東.三軸壓縮下粉砂質泥巖蠕變本構模型研究[J].巖土力學,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].中南大學學報(自然科學版),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].中南大學學報(自然科學版),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].中南大學學報(自然科學版),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]張志沛,王芝銀,彭惠.陜南泥巖三軸壓縮蠕變試驗及其數值模擬研究[J].水文地質工程地質,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]熊良宵,楊林德,張堯.硬巖的復合黏彈塑性流變模型[J].中南大學學報(自然科學版),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].山東科技大學學報(自然科學版),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].巖石力學與工程學報,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]張曉平,吳順川,王思敬.類土質路塹邊坡動態(tài)監(jiān)測及數值模擬分析[J].巖石力學與工程學報,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.