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

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

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

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

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

1 竹林溝滑坡狀況

1.1 竹林溝滑坡地質(zhì)概況

圖1 竹林溝滑坡縱剖面

古藺縣位于亞熱帶季風(fēng)氣候區(qū)，云雨天氣較多，月平均降雨量如圖2所示。據(jù)已有資料[14]顯示，在2000年一次降雨中，滑坡后緣一戶居民房出現(xiàn)裂縫，裂縫走向大致與主滑方向垂直，并于2007年一次連續(xù)兩天降雨中(降雨量約200 mm)墻體倒塌，同時(shí)另外十幾戶居民房出現(xiàn)不同程度的開裂(圖3)。2008年汶川地震后，當(dāng)?shù)卣M織人員對滑坡進(jìn)行過安全評估，結(jié)果顯示滑體在地震情況下并未出現(xiàn)大的變形。2014年當(dāng)?shù)卣畬逻M(jìn)行治理，在變形較大部位布設(shè)抗滑樁(樁埋深7 m)，但據(jù)當(dāng)?shù)卮迕窠榻B，坡體在連續(xù)降雨或大暴雨情況下仍在變形，持續(xù)至今。由此可見降雨是坡體產(chǎn)生變形的主要原因?；乱呀?jīng)過多年變形，累積變形較大，存在重大隱患，對其在降雨作用下的穩(wěn)定性研究顯得尤為必要。

圖2 古藺縣月平均降雨量

圖3 現(xiàn)場房屋受損情況

1.2 滑坡滑動(dòng)原因分析

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

圖4 竹林溝滑坡全貌

(1)當(dāng)?shù)鼐用裨谄麦w上修建房屋、開墾荒地和開展其他工程活動(dòng)，造成坡體表面土質(zhì)疏松裸露，為降雨條件下的雨水下滲創(chuàng)造了有利條件。雨水下滲的過程中使得坡體飽水加載，導(dǎo)致坡體下滑力增加而抗剪強(qiáng)度降低。

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

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

2 微觀參數(shù)選取過程

2.1 顆粒流離散元理論簡介

(1)

(2)

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

(3)

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

顆粒在整個(gè)計(jì)算過程中，循環(huán)運(yùn)用牛頓運(yùn)動(dòng)定律(Newton’s laws of motion)和胡克定律(Hooke’s law)，根據(jù)牛頓運(yùn)動(dòng)定律不斷變換顆粒之間和顆粒與墻體之間的接觸關(guān)系，根據(jù)胡克定律確定顆粒新的接觸關(guān)系中的接觸力(圖5)。當(dāng)顆粒之間或顆粒與墻體之間的作用力大于黏結(jié)強(qiáng)度時(shí)，顆粒就會(huì)發(fā)生相對錯(cuò)動(dòng)。

圖5 顆粒計(jì)算循環(huán)過程

2.2 模型微觀參數(shù)的選取

離散元 PFC2D采用的是微觀的力學(xué)參數(shù)，通過建立正確的微觀和宏觀之間的力學(xué)關(guān)系，把滑坡體宏觀現(xiàn)象用微觀參數(shù)表達(dá)出來。微觀參數(shù)的選取通常是在 PFC2D中進(jìn)行反復(fù)的虛擬雙軸壓縮試驗(yàn)，得到應(yīng)力-應(yīng)變曲線與室內(nèi)壓縮試驗(yàn)曲線進(jìn)行匹配，最終得出滑坡體的微觀力學(xué)參數(shù)。由于坡體滑動(dòng)的主要誘發(fā)因素是降雨，因此直接對坡體土樣在飽和狀態(tài)下的強(qiáng)度參數(shù)進(jìn)行雙軸試驗(yàn)標(biāo)定以模擬持續(xù)降雨工況，坡體的宏觀力學(xué)參數(shù)如表 1 所示。

表1 土的宏觀力學(xué)參數(shù)

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

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

表2 雙軸數(shù)值模型細(xì)觀參數(shù)

圖6 雙軸壓縮應(yīng)力應(yīng)變曲線

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

圖7 莫爾破壞包絡(luò)線

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

3.1 模型的建立

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

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

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

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

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

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

圖8 滑坡ball-wall模型及測點(diǎn)分布

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

3.2 滑坡破壞過程分析

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

圖9 位移云圖

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

3.3 滑坡運(yùn)動(dòng)特征分析

為了能夠精確地研究其滑動(dòng)過程中各部位的發(fā)展變化及力學(xué)響應(yīng)情況，在坡體中分別布置 7 個(gè)控制測量點(diǎn)及監(jiān)測區(qū)域如圖 8 所示。控制點(diǎn)用來監(jiān)測坡體位移，監(jiān)測區(qū)域用來測量滑動(dòng)過程的應(yīng)力變化情況，控制點(diǎn)從右到左序號(hào)分別為 1～7，若以坡頂點(diǎn)在水平面上投影為坐標(biāo)原點(diǎn)，監(jiān)測點(diǎn)坐標(biāo)分別為(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)。模型將控制點(diǎn)布置在監(jiān)測區(qū)域中，通過分析同一位置位移及應(yīng)力的發(fā)展趨勢對坡體變形破壞特征進(jìn)行研究。

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

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

圖10 滑坡測量點(diǎn)位移曲線

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

圖11 各測量圓內(nèi)應(yīng)力變化曲線

4 結(jié)論

根據(jù)現(xiàn)場勘查資料，結(jié)合顆粒流 PFC2D技術(shù)對竹林溝土質(zhì)滑坡在持續(xù)降雨工況下的發(fā)展破壞進(jìn)行模擬研究，得到結(jié)論如下。

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

(2)滑坡各區(qū)域變形存在明顯的不協(xié)調(diào)，后緣張拉裂縫貫通后，失去前部的牽引作用，僅在自重應(yīng)力下產(chǎn)生位移。二級滑動(dòng)也使得中-前部與中-后部呈兩組不同發(fā)展趨勢。

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