基于顆粒流的土質(zhì)邊坡失穩(wěn)機(jī)理

李龍起， 王 滔， 趙瑞志， 趙皓璆， 王夢(mèng)云

(成都理工大學(xué)， 地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國家重點(diǎn)實(shí)驗(yàn)室， 成都 610059)

滑坡是指在降雨、地震、地下水、人類工程活動(dòng)等作用下，巖土體沿一定軟弱面向下滑動(dòng)的自然現(xiàn)象[1]。目前在滑坡的研究上，一般采用物理模型試驗(yàn)、數(shù)值模擬以及理論分析等，室內(nèi)物理模型試驗(yàn)易受到尺寸效應(yīng)或邊界效應(yīng)等影響，理論分析只是采用近似或者定性的分析，因此現(xiàn)階段中外大量專家學(xué)者多采用數(shù)值模擬手段對(duì)滑坡破壞過程進(jìn)行研究。由于土質(zhì)邊坡是復(fù)雜的非連續(xù)介質(zhì)，相比于運(yùn)用于連續(xù)介質(zhì)的FLAC模擬軟件以及應(yīng)用于UDEC軟件，PFC2D數(shù)值模擬技術(shù)在模擬土質(zhì)邊坡變形過程中不同部位的應(yīng)力狀態(tài)、速度、位移等出現(xiàn)的變化存在明顯優(yōu)勢(shì)。宋浩燃等[2]通過PFC3D對(duì)某棄渣邊坡的變形以及滑動(dòng)破壞過程進(jìn)行了模擬，取得了較好的數(shù)值結(jié)果，結(jié)果表明棄渣邊坡的變形集中在邊坡表層，而破壞的規(guī)模受邊坡角度的影響較大。王宇等[3]將平行黏結(jié)模型引入白河滑坡，對(duì)白河滑坡的破壞進(jìn)行了模擬，明確了白河滑坡的時(shí)空演化規(guī)律，其也為后續(xù)工程決策提供了理論依據(jù)。曹文等[4]以紅巖石滑坡為例，使用了Brick分別建立二維和三維模型，提出了一個(gè)比較健全的建模方法，對(duì)后續(xù)的滑坡模擬建模提供了有效的方法參考。Wei等[5]基于支持向量機(jī)對(duì)馬邊滑坡的微觀參數(shù)進(jìn)行了反演，通過高精度數(shù)字高程模型(digital elevation model，DEM)數(shù)據(jù)進(jìn)行了PFC3D建模，并對(duì)滑坡的速度、位移以及能量進(jìn)行了詳細(xì)的分析，數(shù)值效果與實(shí)際特征較為吻合。趙洲等[6]等應(yīng)用顆粒流離散元對(duì)新府山滑坡的破壞過程進(jìn)行了研究，通過PFC2D雙軸壓縮試驗(yàn)來標(biāo)定相應(yīng)的滑坡模型的細(xì)觀參數(shù)，從變形特征和破壞特征分析新府山滑坡的破壞過程以及破壞機(jī)理。

基于上述研究，根據(jù)前期的地質(zhì)調(diào)研對(duì)竹林溝滑坡進(jìn)行了理論分析，隨后對(duì)該滑坡進(jìn)行數(shù)值仿真模擬，模擬結(jié)果與現(xiàn)場(chǎng)實(shí)際情況較為吻合。通過對(duì)滑坡各位置處的速度、位移等參數(shù)的分析，對(duì)滑坡表生改造階段-時(shí)效變形階段-破壞發(fā)展階段[7]進(jìn)行了深入分析，對(duì)該滑坡發(fā)展趨勢(shì)有一定預(yù)測(cè)作用及防災(zāi)減災(zāi)提供理論依據(jù)，也為同類型降雨作用下土質(zhì)邊坡穩(wěn)定性變化提供一定的理論支持。

1 工點(diǎn)地質(zhì)概況

1.1 竹林溝滑坡

竹林溝滑坡位于四川省瀘州市古藺縣土城鄉(xiāng)天井村，106°14′9.56″E，北緯28°4′28.82″N,交通較便利，有公路通過?；赂叱淘?90～750 m，滑坡平面呈長(zhǎng)椅狀，坡體長(zhǎng)×寬=520 m×150 m，覆蓋層厚約12 m，方量約9×105m3，主滑方向280°，滑坡坡度在20°左右，滑坡現(xiàn)場(chǎng)如圖1所示。研究區(qū)地表以殘坡積粉質(zhì)黏土(Q4del)為主，基巖為三疊系飛仙關(guān)組(T1f)為主，巖性為紫紅色砂巖夾泥巖，區(qū)域地質(zhì)構(gòu)造圖如圖2所示。根據(jù)現(xiàn)場(chǎng)勘測(cè)及歷史資料[8]可知，2000年時(shí)，在一次降雨中少數(shù)住戶的地面及墻體出現(xiàn)了細(xì)微裂縫，2007年時(shí)，在連續(xù)降雨的作用下，又有多戶居民住房產(chǎn)生多條張拉裂縫。竹林溝滑坡地質(zhì)剖面圖如圖3所示。

圖1 滑坡特征Fig.1 Landslide characteristics

圖2 區(qū)域地質(zhì)構(gòu)造Fig.2 Regional geological structure

圖3 滑坡剖面圖Fig.3 Landslide profile

1.2 滑坡影響因素以及成因分析

竹林溝滑坡是在其內(nèi)外因素共同作用才形成的。坡體縱向成兩級(jí)臺(tái)地，前緣坡度約25°，后緣略緩，坡度約20°，前緣與后緣的高程差較大。另外，由于前緣河谷下切(表生改造階段)，促使坡腳產(chǎn)生臨空面，為滑坡失穩(wěn)創(chuàng)造了有利的條件。在長(zhǎng)期暴雨工況下，坡內(nèi)土體的含水率增加、重度增大，導(dǎo)致土體力學(xué)性能劣化，使坡體發(fā)生剪切蠕變[9]，并在后緣及中部出現(xiàn)了較多的張拉裂縫，為坡表水滲入坡體內(nèi)部創(chuàng)造了有利條件，進(jìn)一步減弱土體的抗剪強(qiáng)度，最終坡體發(fā)生破壞，破壞模式為蠕滑—拉裂。

2 數(shù)值模擬計(jì)算

2.1 滑坡參數(shù)標(biāo)定

PFC2D進(jìn)行數(shù)值計(jì)算時(shí)采用的是巖土的細(xì)觀參數(shù)，因此必須將宏觀與細(xì)觀參數(shù)建立聯(lián)系，進(jìn)而對(duì)參數(shù)進(jìn)行標(biāo)定[10]，巖土雙軸壓縮試驗(yàn)即可實(shí)現(xiàn)該過程。考慮到竹林溝滑坡的主要誘發(fā)因素是降雨，故將對(duì)土體在暴雨工況下的參數(shù)進(jìn)行標(biāo)定。土體的宏觀力學(xué)參數(shù)如表1所示[7]。

表1 土的宏觀力學(xué)參數(shù)[7]Table 1 Macro-mechanical parameters of soil[7]

建立雙軸壓縮試驗(yàn)時(shí)，考慮到計(jì)算機(jī)的能力以及精度，在數(shù)值模擬時(shí)將計(jì)算顆粒進(jìn)行適當(dāng)放大[11]。數(shù)值試驗(yàn)采用不同粒徑的顆粒組成，為了滿足土體的不均勻性和各向異性，顆粒半徑R服從Gauss分布，分別輸入顆粒的各種參數(shù)[12]，通過對(duì)顆粒的不同參數(shù)進(jìn)行反復(fù)調(diào)試，最終確定顆粒半徑介于0.06～0.08 m。為避免尺寸效應(yīng)[13]，即模型短邊應(yīng)大于40倍的粒徑，最終確定雙軸模型尺寸為5 m×10 m。雙軸壓縮試驗(yàn)?zāi)Ｐ腿鐖D4所示。

擬通過建立顆粒流離散元二維雙軸壓縮模型來模擬室內(nèi)剪切條件下粉質(zhì)黏土的力學(xué)試驗(yàn)，使用PFC2D中的圓形顆粒來模擬試驗(yàn)中的粉質(zhì)黏土，使用剛性墻體來模擬加載設(shè)備，并將PFC中的接觸黏結(jié)模型引入，利用PFC的伺服控制程序[14]控制雙軸壓縮過程中的圍壓，通過大量的試算，最終確定使用表2中的細(xì)觀參數(shù)進(jìn)行模擬，并從PFC軟件中提取出了圍壓分別為100、200、300 kPa的應(yīng)力-應(yīng)變曲線，如圖5所示，隨后根據(jù)應(yīng)力-應(yīng)變曲線繪制摩爾-庫倫應(yīng)力圓，如圖6所示，最終得出土體的黏聚力c=22.3 kPa，內(nèi)摩擦角φ=13.2°，經(jīng)驗(yàn)證該數(shù)據(jù)與室內(nèi)試驗(yàn)結(jié)果相近，故可以用于后期的滑坡運(yùn)動(dòng)過程模擬研究中。

σ1為軸向應(yīng)力；σ3為側(cè)向圍壓圖4 雙軸試驗(yàn)?zāi)Ｐ虵ig.4 Biaxial test model

圖5 不同圍壓下應(yīng)力-應(yīng)變曲線Fig.5 Stress-strain curves under different confining pressures

圖6 摩爾應(yīng)力破壞包線Fig.6 Molar stress failure envelope

表2 細(xì)觀參數(shù)Table 2 Mesoscopic parameters

2.2 滑坡數(shù)值模型建立

PFC2D滑坡模型建立的方法通常有兩種，一種是用墻單元來模擬滑床而滑體使用顆粒單元填充，這種建模方法稱為Ball-Wall法；另一種則是ball-ball法，即用顆粒單元來填充滑體和滑床。根據(jù)現(xiàn)場(chǎng)勘查和相關(guān)資料表明，已經(jīng)確定竹林溝滑坡的潛在滑動(dòng)面；其次，鑒于Ball-Wall法所需顆粒較少，計(jì)算效率會(huì)大大提高。因此將采用Ball-Wall法建立模型建模方法如下。

(1)首先根據(jù)坡體的地質(zhì)剖面圖借助AUTOCAD繪制相同尺寸的邊坡幾何模型。

(2)在PFC中使用geometry import以及wall import等相關(guān)命令將前面生成的幾何模型導(dǎo)入到PFC中并生成相應(yīng)的墻體邊界。

(3)在墻體中使用膨脹法[15]生成相對(duì)應(yīng)孔隙率的顆粒，并經(jīng)過多次試算，使其重力作用下平衡?；鲁跏寄Ｐ腿鐖D7所示。

圖7 滑坡初始模型Fig.7 Landslide initial model

(4)對(duì)生成的土體顆粒賦予表1的相關(guān)參數(shù)，最后刪除相應(yīng)滑體的邊界墻體，因墻體刪除后顆粒之間的應(yīng)力會(huì)使顆粒大量分離，因此必須再進(jìn)行一次迭代平衡，使顆粒的最大不平衡接觸力接近0，表明整個(gè)模型體系已經(jīng)達(dá)到了平衡狀態(tài)，如圖8所示。

圖8 不平衡力圖Fig.8 Unbalanced force diagram

(5)對(duì)滑坡模型進(jìn)行速度和位移的清零，使用set gravity 命令對(duì)模型施加重力，使模型在重力作用下滑動(dòng)，坡體輪廓線只作為坡體滑動(dòng)的參考。

建立的數(shù)值計(jì)算模型如圖9所示，長(zhǎng)520 m，高200 m，坡度約20°，顆粒最小半徑為0.06 m，最大半徑為0.08 m。

A～E為測(cè)量圓編號(hào)圖9 數(shù)值邊坡模型Fig.9 Numerical slope model

3 滑坡數(shù)值模擬結(jié)果分析

3.1 破壞特征

根據(jù)圖10關(guān)鍵位置處土體運(yùn)動(dòng)速度云圖可知，運(yùn)行至2 000時(shí)步時(shí)，坡體中后部的顆粒首先達(dá)到較大速度，主要是由于中部及上部的坡體蠕滑變形累積所造成，坡后緣及中部出現(xiàn)少量的張拉裂隙，并且裂隙有不斷擴(kuò)張的趨勢(shì)，表明坡體處于蠕滑變形階段。運(yùn)行至10 000時(shí)步時(shí)，張拉裂隙向坡內(nèi)延伸貫通，坡腳發(fā)生剪切破壞，整個(gè)坡體在重力作用下向下運(yùn)動(dòng)；坡體不同部位形變速度產(chǎn)生差異，滑坡中部和前緣坡腳的顆粒速度較大，坡體后緣顆粒速度較小，此時(shí)坡體已經(jīng)處于時(shí)效變形階段。運(yùn)行至20 000時(shí)步時(shí)，坡腳顆粒已經(jīng)快速剪出，帶動(dòng)中部及后緣坡體向前運(yùn)動(dòng)。由于自后緣到前緣的坡度逐漸增大，因此后緣坡體速度增加較快。

根據(jù)位移云圖(圖11)也可以看出，當(dāng)運(yùn)行至30 000時(shí)步時(shí)，滑坡前緣部分坡體在中部坡體的擠壓作用下向前剪出。滑坡中部以及后緣張拉裂縫逐漸增大，后緣土體有向前滑塌的趨勢(shì)，這也說明滑坡已經(jīng)處于漸進(jìn)破壞發(fā)展階段[7]。由于張拉裂隙增大，坡體不同位置顆粒的位移也有所不同，整體呈現(xiàn)前緣大于后緣的趨勢(shì)。運(yùn)行至80 000時(shí)步，此時(shí)滑動(dòng)面已經(jīng)貫通，坡體已經(jīng)整體向前滑動(dòng)。運(yùn)行至300 000時(shí)步時(shí)，坡體滑出滑床，速度逐漸降低，位移不在增加，最后滑體堆積在地面。

3.2 運(yùn)動(dòng)特征

在PFC模型計(jì)算過程中，無法對(duì)坡體進(jìn)行整體的變形、運(yùn)動(dòng)特征的監(jiān)測(cè)，因此，分別選取滑坡中幾個(gè)比較關(guān)鍵的位置設(shè)置5個(gè)半徑為2 m的測(cè)量圓，如圖9中紅色圓點(diǎn)所示，測(cè)量圓自后緣到前緣編號(hào)依次為A、B、C、D、E。其測(cè)量點(diǎn)主要是用于監(jiān)測(cè)關(guān)鍵位置處土體孔隙率的實(shí)時(shí)變化進(jìn)而分析滑坡的破壞過程。另外，記錄坡體上部，中部，下部5個(gè)顆粒的速度、位移變化，顆粒坐標(biāo)分別為(38.746 7，200.214)，(96.631 3，182.98)，(249.96，120.184)，(344.047，85.154)，(453.484，56.970 8)，編號(hào)以依次為1、2、3、4、5，通過監(jiān)測(cè)顆粒的速度和位移變化來分析整個(gè)坡體的運(yùn)動(dòng)過程。由此可以得到多個(gè)曲線圖，進(jìn)而與模型邊坡的整體失穩(wěn)破壞過程相結(jié)合進(jìn)行分析，探究坡體形變與特殊物理量之間的內(nèi)在聯(lián)系。

通過觀察圖12孔隙率-時(shí)步曲線可知，在模擬開始時(shí)，5個(gè)測(cè)量圓中的孔隙率都較接近初始給定孔隙率，但整個(gè)坡體由于重力的作用，顆粒之間會(huì)相互擠壓，導(dǎo)致孔隙率有略微的下降；隨著滑坡體向下運(yùn)動(dòng)，坡腳剪出較快，因此坡腳的2個(gè)測(cè)量圓孔隙率會(huì)逐漸上升，而中部的測(cè)量圓由于坡體的鎖固作用，土體會(huì)繼續(xù)被擠壓密實(shí)，孔隙率會(huì)逐漸降低，一旦滑動(dòng)剪切面貫通，鎖固段擠壓作用消失，孔隙率則會(huì)上升。坡體后緣由于張拉裂隙的產(chǎn)生，孔隙率會(huì)逐漸上升。其中，坡腳的孔隙率增長(zhǎng)較快，說明坡腳顆?？赡軙?huì)發(fā)生較大的位移。

圖12 孔隙率變化曲線Fig.12 Porosity curve

通過觀察圖13速度-時(shí)步曲線可知：不同部位的顆粒的運(yùn)動(dòng)速度變化趨勢(shì)整體上較相近，均表現(xiàn)為先增后減，但不同部位顆粒達(dá)到峰值的時(shí)刻不同，整體呈現(xiàn)由前緣到后緣依次達(dá)到速度峰值，高程越大，達(dá)到的速度峰值也越大。坡體不同位置的速度峰值介于0～20 m/s，其中位于坡體后緣的顆粒的速度峰值最大，17 s時(shí)到達(dá)了20 m/s?；瑒?dòng)初期時(shí)，由于后緣坡度較緩，1號(hào)顆粒的加速度較其他部位小，同一時(shí)刻速度較小?；潞缶壍?號(hào)顆粒和滑坡前緣的5號(hào)顆粒的斜率較大，表明兩個(gè)顆粒的加速度比較大，可能會(huì)產(chǎn)生較大位移。而隨著時(shí)間的推移，2號(hào)顆粒的速度仍然增長(zhǎng)較快，表明后緣局部形成張拉裂隙。伴隨整個(gè)剪切滑動(dòng)面貫通，坡體中部的2、3、4顆粒的速度均在15 s左右達(dá)到速度的峰值。

圖13 速度曲線Fig.13 Velocity curve

根據(jù)坡內(nèi)關(guān)鍵位置處5個(gè)顆粒運(yùn)動(dòng)速度，繪制出坡體平均速度變化時(shí)程曲線，如圖14所示，其曲線表明：0～10 s平均速度近似線性增大，15 s時(shí)達(dá)到速度峰值14 m/s，達(dá)到峰值后平均速度開始下降，大約達(dá)到21 s左右，滑坡平均速度下降變緩，但顆粒最終仍然有一個(gè)較小的速度，為0.4 m/s，隨著計(jì)算時(shí)間的增加，坡體的速度將趨近于0。

圖14 平均速度曲線Fig.14 Average speed curve

通過觀察圖15顆粒的位移-時(shí)間曲線可知，處于滑坡后緣的1號(hào)顆粒的位移最大，達(dá)到了600 m，處于滑坡前緣的4號(hào)顆粒的位移最小，最小位移為310 m，所監(jiān)測(cè)的顆粒達(dá)到最大位移的時(shí)間均在25 s左右。滑坡整體的位移趨勢(shì)是顆粒所處的高度越高，則最終達(dá)到的位移也就越大。

圖15 位移曲線Fig.15 Displacement curve

通過觀察位移變化云圖，可以認(rèn)為每個(gè)顆粒位移變化可以代表其附近區(qū)域的位移變化特征。根據(jù)滑坡內(nèi)5個(gè)顆粒的位移曲線，繪制相應(yīng)的平均位移變化時(shí)程曲線，如圖16所示，根據(jù)平均位移變化曲線可以看出，0～10 s范圍內(nèi)，平均位移斜率逐漸增大，說明速度增加較快，因此位移也逐漸增大，該階段為加速階段。10～18 s平均位移呈現(xiàn)線性增大的趨勢(shì)，說明滑坡此時(shí)已經(jīng)完全破壞，整體向下滑動(dòng)。18～25 s，平均位移的斜率逐漸減小，坡體的速度開始下降，平均位移仍然增大，該階段為滑坡的減速階段。25 s后平均位移基本不在變化，說明滑坡已經(jīng)停止滑動(dòng)，且最大平均位移為450 m。

圖16 平均位移曲線Fig.16 Average displacement curve

4 結(jié)論

以竹林溝滑坡為依托，結(jié)合數(shù)值模擬軟件PFC2D建立計(jì)算模型，選用Ball-Wall的方法進(jìn)行建模，并將接觸黏結(jié)模型概念引入到模型中，由此研究該滑坡的破壞模式及坡內(nèi)土體的力學(xué)響應(yīng)，探究坡體破壞模式與重要物理量之間的內(nèi)在聯(lián)系，得出如下結(jié)論。

(1)滑坡是在持續(xù)降雨工況下破壞的，整體破壞特征為：前緣牽引-中部剪切-后緣拉裂。表現(xiàn)為典型的蠕滑-拉裂式破壞。

(2)滑坡初期由于顆粒之間的擠壓，坡體的孔隙率均有所下降，隨著計(jì)算時(shí)步的增加，坡體前緣開始剪出，坡體不同位置逐漸出現(xiàn)張拉裂隙，孔隙率會(huì)逐漸增大。

(3)通過對(duì)滑坡的速度進(jìn)行監(jiān)測(cè)，可以得到：該滑坡整個(gè)破壞過程持續(xù)25 s左右，0～10 s為加速階段，滑坡平均速度達(dá)14 m/s，局部最大速度達(dá)到20 m/s，15～25 s為減速階段，最終在25 s左右坡體速度趨近于0。滑坡整體位移約450 m，處于后緣的顆粒位移最大，約600 m，前緣顆粒位移較小，最小約310 m。