降雨作用下的土質邊坡變形破壞顆粒流仿真模擬

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

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

滑坡是指斜坡巖土體在各種內外營力作用下沿著貫通破裂面發生的以水平位移為主的滑動現象，具有突發性和極強的破壞性。同時，由于坡體本身具有不同的結構和復雜的物質組成情況，使滑坡體的變形破壞模式多種多樣。特別是中國幅員遼闊，地形復雜，經常發生滑坡災害，例如四川宣漢天臺鄉滑坡[1]、重慶武隆縣雞尾山滑坡[2]、四川茂縣新磨村滑坡[3]等均造成人員及財產損失。

自然界多數滑坡的發生都與降雨密切相關，因此，眾多學者通過不同維度對降雨型滑坡進行了不斷的探索與總結[4-6]。對于滑坡的研究，一般是通過物理試驗、數值模擬并結合理論分析來進行，其中物理試驗通常會受到邊界效應、尺寸效應等因素的影響，且效率低、工作量大；而數值模擬憑借高效率、低成本和可重復性好的優點應用越來越廣泛。前人對數值模擬在滑坡中的應用進行了不斷地探索，但大多都是基于有限元的基礎上進行[7-8]，而滑坡體具有不均勻性、非連續性和離散性的特點，通過有限元方法在求解位移不連續及大變形問題上存在不足。Cundall[9]提出離散元分析方法為滑坡數值模擬提供了新途徑，并由此產生了顆料流(PFC)、離散單元法程序(UDEC)、三維離散單元法程序(3DEC)等離散元軟件。PFC2D相對于其他離散元軟件有效率高和可以有效模擬大變形的優點，因此得到了學者廣泛關注。畢鈺璋等[10]應用PFC2D離散元軟件對牛眠溝高速遠程滑坡全過程進行模擬，得到不同摩擦系數下的滑坡變形及堆積特征；陳達等[11]運用PFC2D以陜西省洛南縣劉澗滑坡為研究對象,得到不同位置特征點的加速度隨時間的變化情況；周健等[12]運用顆粒流PFC2D分別對黏性邊坡和砂性邊坡進行數值模擬，結果表明土的性質對邊坡的破壞形式有很大的影響，隨著顆粒黏性的減小邊坡的破壞形式從脆性向塑性轉化；趙洲等[13]運用PFC2D對楊家灣堆積層滑坡的破壞過程進行數值模擬，得到堆積層滑坡的破壞模式是一種典型的牽引式漸進破壞過程。目前，學者們運用顆粒流對滑坡進行的研究主要集中在滑坡的破壞過程及各種因素對堆積形式的影響，并沒有對處于時效變形階段的坡體在降雨或其他工況下的發展趨勢進行研究。

以四川省古藺縣的竹林溝滑坡為例，通過現場勘查、理論分析和數值模擬對滑體進行綜合研究。運用顆粒流離散元軟件PFC2D對滑動體在持續降雨作用下的變形破壞過程進行仿真模擬，并通過在模型上設置監測點，得到坡內不同位置的變形特征，由此分析預測滑坡的發展趨勢，為滑坡防治提供科學依據。

1 竹林溝滑坡狀況

1.1 竹林溝滑坡地質概況

圖1 竹林溝滑坡縱剖面

古藺縣位于亞熱帶季風氣候區，云雨天氣較多，月平均降雨量如圖2所示。據已有資料[14]顯示，在2000年一次降雨中，滑坡后緣一戶居民房出現裂縫，裂縫走向大致與主滑方向垂直，并于2007年一次連續兩天降雨中(降雨量約200 mm)墻體倒塌，同時另外十幾戶居民房出現不同程度的開裂(圖3)。2008年汶川地震后，當地政府組織人員對滑坡進行過安全評估，結果顯示滑體在地震情況下并未出現大的變形。2014年當地政府對滑坡進行治理，在變形較大部位布設抗滑樁(樁埋深7 m)，但據當地村民介紹，坡體在連續降雨或大暴雨情況下仍在變形，持續至今。由此可見降雨是坡體產生變形的主要原因。滑坡已經過多年變形，累積變形較大，存在重大隱患，對其在降雨作用下的穩定性研究顯得尤為必要。

圖2 古藺縣月平均降雨量

圖3 現場房屋受損情況

1.2 滑坡滑動原因分析

竹林溝滑坡的成因與眾多因素有關，滑坡全貌如圖 4 所示。根據已有資料[14]及現場勘查分析，竹林溝滑坡的成因主要有以下3個方面。

圖4 竹林溝滑坡全貌

(1)當地居民在坡體上修建房屋、開墾荒地和開展其他工程活動，造成坡體表面土質疏松裸露，為降雨條件下的雨水下滲創造了有利條件。雨水下滲的過程中使得坡體飽水加載，導致坡體下滑力增加而抗剪強度降低。

(2)滑坡體整體坡度較大(其中中后緣坡度為15°～30°，中部坡度為15°～40°，中前緣坡度相對較緩，為10°～20°)，自重應力沿滑動方向的分力較大，對坡體穩定性不利。坡腳處由于河谷下切，使得前緣臨空，缺乏有效的支擋結構。

(3)由于人類活動頻繁，生活排水、稻田灌溉及工程活動等均對滑體產生不利影響。而在坡體中后緣由于道路修建進行的大面積的開挖，破壞了坡體的整體性，同樣也加劇了坡體的破壞進程。

2 微觀參數選取過程

2.1 顆粒流離散元理論簡介

(1)

(2)

式中：kn和ks分別為接觸點法向和切向剛度；un和us分別為接觸點法向和切向的微小位移。

(3)

式(3)中：c為顆粒接觸點的黏聚力；φ為顆粒接觸點的摩擦角。

顆粒在整個計算過程中，循環運用牛頓運動定律(Newton’s laws of motion)和胡克定律(Hooke’s law)，根據牛頓運動定律不斷變換顆粒之間和顆粒與墻體之間的接觸關系，根據胡克定律確定顆粒新的接觸關系中的接觸力(圖5)。當顆粒之間或顆粒與墻體之間的作用力大于黏結強度時，顆粒就會發生相對錯動。

圖5 顆粒計算循環過程

2.2 模型微觀參數的選取

離散元 PFC2D采用的是微觀的力學參數，通過建立正確的微觀和宏觀之間的力學關系，把滑坡體宏觀現象用微觀參數表達出來。微觀參數的選取通常是在 PFC2D中進行反復的虛擬雙軸壓縮試驗，得到應力-應變曲線與室內壓縮試驗曲線進行匹配，最終得出滑坡體的微觀力學參數。由于坡體滑動的主要誘發因素是降雨，因此直接對坡體土樣在飽和狀態下的強度參數進行雙軸試驗標定以模擬持續降雨工況，坡體的宏觀力學參數如表 1 所示。

表1 土的宏觀力學參數

通過 PFC2D軟件建立雙軸壓縮試驗的模型，隨后在模型內部填充顆粒。考慮到實際顆粒大小具有不均勻性，同時保證計算效率，最終選擇粒徑介于0.36～0.54 m的顆粒進行填充，孔隙率為0.155。為了避免尺寸效應對試驗造成影響，模型的短邊尺寸還應大于40倍的顆粒平均半徑[16]。各粒徑顆粒數量從Rmin到Rmax服從正態分布。

模型建立完成對其施加圍壓，并通過頂、底板施加軸向壓力進行試驗。經過大量虛擬雙軸壓縮試驗得到多條應力-應變曲線，與被模擬材料飽和狀態下的室內試驗所得結果匹配，最終得出匹配程度最高的應力-應變曲線如圖6所示。對應微觀參數的取值如表2所示。

表2 雙軸數值模型細觀參數

圖6 雙軸壓縮應力應變曲線

圖7為模型賦予參數后在不同圍壓下(50、100、150 kPa)模擬雙軸壓縮實驗所得出的莫爾破壞包絡線，最終反算出的宏觀參數內摩擦角φ為16.3°，黏聚力c為21.7 kPa，得到的結果與被模擬滑坡土體在持續降雨工況下的c、φ相近，這也論證了所選微觀參數數值的可行性。

圖7 莫爾破壞包絡線

3 數值模擬結果分析

3.1 模型的建立

PFC 提供兩種建模方法，ball-ball 法和 ball-wall 法。ball-ball 法是滑坡體和基巖均采用顆粒進行填充，因此所用顆粒較多，計算過程較長，適用于滑動面未知的情況。ball-wall 法是滑坡體采用顆粒填充而基巖采用墻體進行模擬，所用顆粒較少，運算效率高，適用于滑動面已知的情況。根據文獻[14]和現場勘查已確定滑坡潛在滑動面的位置，模型采用 ball-wall 法進行建模。具體建模過程如下。

(1)在 AutoCAD 中繪制滑體和滑床的邊界線，尺寸與被模擬滑坡真實尺寸一致，將圖層存儲為 DXF 文件格式。

(2)在 PFC 中使用 geometry import 命令和wall import geometry 命令將 DXF 文件導入后根據滑體和滑床邊界線生成邊界墻。

(3)在邊界墻形成的空間內填入足夠多的 Brick，并根據表 2 賦予顆粒接觸剛度、摩擦系數等微觀參數。

(4)刪除超出邊界線以外的顆粒，進行多次迭代計算使得模型達到平衡，并在坡體上設置需要的監測區域。

最終所建模型如圖 8 所示，模型運用 PFC 軟件自帶的塊體填充法(brick)在墻體內部填充足夠多的顆粒[17]，該方法規避了膨脹法由于膨脹系數難以準確選擇，而導致顆粒接觸處應力過大，顆粒飛出的現象，且操作簡單。

圖8 滑坡ball-wall模型及測點分布

模型建立完成后，將模型的初始位移和速度置零并去除邊界約束墻，對坡體施加重力使坡體在重力作用下開始緩慢變形。

3.2 滑坡破壞過程分析

滑坡體不同時步的位移云圖如圖9所示。模型坡體邊界線只作為坡體變形的參考。

圖9 位移云圖

由圖9(a)可知，4 000 時步時坡體并未出現較大變形,僅在中部一區域內產生小變形，說明滑坡的變形是從中部開始的。坡體上、下兩部分覆蓋層厚度及滑床坡度均存在差異，導致 12 000 時步時滑體表現出二級滑動特征，一級滑動面高程約 570 m，二級滑動面高程約 680 m。一級滑動發生后，坡體中部產生拉張裂縫，為二級滑動創造了臨空條件。20 000 時步時坡腳處有部分土體在上部擠壓作用下剪出，坡體中部和后緣拉張裂縫增大，后緣出現“下坐”現象圖9(c)，這是因為后緣裂縫增大到一定程度，使得裂縫后土體失去有效支擋從而向前傾倒，伴隨著裂縫的縮小和土體的“下坐”現象，這種現象也表明坡體進入累進性破壞階段[18]，進入該階段后坡體已不能再承受過多變形，其之后的變化發展也具有不可控性，當潛在剪切破壞面的鎖固段被貫通就會發生滑坡。圖9(d)為滑坡在 100 000 時步的發展變化情況，此時滑坡剪切面已經被貫通，坡體發生整體滑動，坡腳土體滑出數十米。根據滑坡破壞特征可知，該滑坡的破壞模式是典型的蠕滑-拉裂。

3.3 滑坡運動特征分析

為了能夠精確地研究其滑動過程中各部位的發展變化及力學響應情況，在坡體中分別布置 7 個控制測量點及監測區域如圖 8 所示。控制點用來監測坡體位移，監測區域用來測量滑動過程的應力變化情況，控制點從右到左序號分別為 1～7，若以坡頂點在水平面上投影為坐標原點，監測點坐標分別為(437.28，61.34)、(362.75，80.83)、(273.95，111.06)、(212.24，137.82)、(141.65，162.89)、(86.67，186.81)、(39.58，198.86)。模型將控制點布置在監測區域中，通過分析同一位置位移及應力的發展趨勢對坡體變形破壞特征進行研究。

由于滑坡一旦進入累進性破壞，其穩定性無論在工程措施還是時間上都極難甚至不可能得到控制[19]，因此滑坡防治的最佳時間是在時效變形階段，所以對時效變形階段進行研究更具意義。由此，僅對滑坡進入累進性破壞前(20 000 時步前)的位移及應力的變化過程進行研究，以向滑坡治理提供科學依據，對滑坡堆積過程不做詳細分析。

圖10為各監測點的位移時程曲線。由圖10可知，在 5 000 時步之前滑坡各測點位移曲線斜率均較小，且位移的增長基本保持一致，這表明在滑坡浸水飽和后穩定性降低，滑坡發生整體的微小位移。在 5 000 時步后 7 號測點位移曲線的增長速率明顯減緩，且之后該點位移量與其余監測點位移差距逐漸增大，這是因為坡體在滑動過程中存在不協調變形[20]，滑坡前部滑動較快，而坡體后緣坡度緩、滑動慢，促使后緣產生拉張裂隙，隨著裂隙不斷增大，向坡內延伸至基巖層面后，滑坡后緣完全失去了前部的牽引作用，僅在自重作用下滑動，因此 7 號測點區域土體與前部的位移差越來越大。12 000 時步之后，坡體各區域土體位移變化量差異明顯增大，表現為從后緣到前緣位移量漸增的特征，且1、2、3號成一組變化，4、5、6號成一組變化，一方面是由于此時坡體出現二級滑動趨勢，另一方面是坡體中-后段變形受前部土體抗滑力約束，而中-前段土體發生移動后產生臨空面，為坡體后緣土體的失穩提供變形空間，即表現為從下到上的牽引式破壞。在20 000時步時7號測點位移量達到 1.05 m，1～6 號測量圓位移量達到1.65 m左右。

圖10 滑坡測量點位移曲線

圖11為各測量圓 20 000 時步內的應力曲線圖，變形初期，坡體中部位移較前部大，致使 1、2 號測量圓區域應力水平不斷積累，在 6 000 時步左右，坡體整體出現較大位移，1、2 號測量圓區域應力得到釋放，應力值陡降。其中，1 號測量圓位于坡腳處，坡度較緩，因此X方向應力明顯大于Y方向，約為Y方向的1.5倍；2號測量圓所在位置坡度較1號陡，因此X方向應力略大于Y方向，1、2號測量圓的峰值應力均產生在5 500時步。由于滑床坡度整體表現為前緩后陡，因此3、4、5號測量圓X、Y方向的應力差依次縮小，應力大小表現為σ3>σ4>σ5。6號測量圓位于滑坡后緣受拉區附近，受到中部坡體產生的拉應力作用，表現為X方向應力大于Y方向應力。7號測量圓滑床坡度緩，因此測量圓Y方向的應力比X方向大。

圖11 各測量圓內應力變化曲線

4 結論

根據現場勘查資料，結合顆粒流 PFC2D技術對竹林溝土質滑坡在持續降雨工況下的發展破壞進行模擬研究，得到結論如下。

(1)坡體在持續降雨工況下的破壞是從坡腳開始的，隨后發生從下到上的牽引式漸進破壞。中-前段坡體破壞為中-后段坡體提供了臨空條件，因此坡體破壞表現出二級滑動特征，滑坡整體上表現為蠕滑-拉裂破壞模式。

(2)滑坡各區域變形存在明顯的不協調，后緣張拉裂縫貫通后，失去前部的牽引作用，僅在自重應力下產生位移。二級滑動也使得中-前部與中-后部呈兩組不同發展趨勢。

(3)坡體前緣應力和位移均大于后緣，一級滑動發生后坡腳應力得到釋放，前緣應力水平有所下降，二級滑動受到前部土體阻擋，應力水平下降較少。