基于顆粒流的土質邊坡失穩機理

李龍起， 王 滔， 趙瑞志， 趙皓璆， 王夢云

(成都理工大學， 地質災害防治與地質環境保護國家重點實驗室， 成都 610059)

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

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

1 工點地質概況

1.1 竹林溝滑坡

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

圖1 滑坡特征Fig.1 Landslide characteristics

圖2 區域地質構造Fig.2 Regional geological structure

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

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

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

2 數值模擬計算

2.1 滑坡參數標定

PFC2D進行數值計算時采用的是巖土的細觀參數，因此必須將宏觀與細觀參數建立聯系，進而對參數進行標定[10]，巖土雙軸壓縮試驗即可實現該過程。考慮到竹林溝滑坡的主要誘發因素是降雨，故將對土體在暴雨工況下的參數進行標定。土體的宏觀力學參數如表1所示[7]。

表1 土的宏觀力學參數[7]Table 1 Macro-mechanical parameters of soil[7]

建立雙軸壓縮試驗時，考慮到計算機的能力以及精度，在數值模擬時將計算顆粒進行適當放大[11]。數值試驗采用不同粒徑的顆粒組成，為了滿足土體的不均勻性和各向異性，顆粒半徑R服從Gauss分布，分別輸入顆粒的各種參數[12]，通過對顆粒的不同參數進行反復調試，最終確定顆粒半徑介于0.06～0.08 m。為避免尺寸效應[13]，即模型短邊應大于40倍的粒徑，最終確定雙軸模型尺寸為5 m×10 m。雙軸壓縮試驗模型如圖4所示。

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

σ1為軸向應力；σ3為側向圍壓圖4 雙軸試驗模型Fig.4 Biaxial test model

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

圖6 摩爾應力破壞包線Fig.6 Molar stress failure envelope

表2 細觀參數Table 2 Mesoscopic parameters

2.2 滑坡數值模型建立

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

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

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

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

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

(4)對生成的土體顆粒賦予表1的相關參數，最后刪除相應滑體的邊界墻體，因墻體刪除后顆粒之間的應力會使顆粒大量分離，因此必須再進行一次迭代平衡，使顆粒的最大不平衡接觸力接近0，表明整個模型體系已經達到了平衡狀態，如圖8所示。

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

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

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

A～E為測量圓編號圖9 數值邊坡模型Fig.9 Numerical slope model

3 滑坡數值模擬結果分析

3.1 破壞特征

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

根據位移云圖(圖11)也可以看出，當運行至30 000時步時，滑坡前緣部分坡體在中部坡體的擠壓作用下向前剪出?；轮胁恳约昂缶墢埨芽p逐漸增大，后緣土體有向前滑塌的趨勢，這也說明滑坡已經處于漸進破壞發展階段[7]。由于張拉裂隙增大，坡體不同位置顆粒的位移也有所不同，整體呈現前緣大于后緣的趨勢。運行至80 000時步，此時滑動面已經貫通，坡體已經整體向前滑動。運行至300 000時步時，坡體滑出滑床，速度逐漸降低，位移不在增加，最后滑體堆積在地面。

3.2 運動特征

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

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

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

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

圖13 速度曲線Fig.13 Velocity curve

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

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

通過觀察圖15顆粒的位移-時間曲線可知，處于滑坡后緣的1號顆粒的位移最大，達到了600 m，處于滑坡前緣的4號顆粒的位移最小，最小位移為310 m，所監測的顆粒達到最大位移的時間均在25 s左右?；抡w的位移趨勢是顆粒所處的高度越高，則最終達到的位移也就越大。

圖15 位移曲線Fig.15 Displacement curve

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

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

4 結論

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

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

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

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