基于PFC2D的柿子樹坪滑坡變形特征模擬研究

惠航，劉鋒，馮兵，趙法鎖

(1.中國地震局第二監測中心，西安 710054；2.陜西核工業工程勘察院有限公司，西安 710054;3.長安大學地質工程與測繪學院，西安 710054)

秦巴山區地質構造復雜，地形變化大，季節性降雨明顯且多集中于夏秋季，是滑坡、泥石流等災害高發地帶[1],隨著基礎設施建設的深入,大量的滑坡、泥石流等地質災害事件的頻發,其中由于滑坡造成的人員傷亡人數呈逐漸上升趨勢[2-3]，滑坡的研究與防治工作依舊繁重。

近些年，基于數值模擬方法的滑坡致災機理分析以及危險性評價是滑坡災害研究的重點之一[4-5]。目前多以有限元、有限差分等方法分析為主，而滑坡體具有不均勻性、非連續性和離散性的特點，通過以上方法在求解位移不連續及大變形問題上稍顯不足。20世紀70年代，Cundall[6]提出了適用于非連續介質動力特性的研究的離散單元法，相對于以上兩種方法有效率高和可以有效模擬大變形的優點，為滑坡數值模擬提供了新方法。周健等[7]利用顆粒流的方法研究細觀參數取值差異對于土質邊坡破壞形式的影響,展示了顆粒流應用用于觀察巖土體的破壞過程的優勢。文獻[8-11]研究了滑坡變形破壞機理以及強度定量化，得出滑坡破壞的突發性與堆積的緩慢性，并得出滑坡強度隨坡腳距的增加而呈現指數衰減特征的沖擊力指標。文獻[12-13]從某含軟弱夾層的邊坡入手，用PFC2D模擬滑坡的漸進失穩過程,得出當靠近坡面的軟弱結構面發生滑動時，因為上覆滑體的擠壓和摩擦作用使得滑床巖體破碎,并產生不斷向深處發展的滑移機理。文獻[14-15]對高速遠程類滑坡運動的全過程進行顆粒流模擬，得到不同坡面摩擦因數下的滑坡變形及堆積特征。文獻[16]對巖層傾角與邊坡走向夾角變化時，邊坡的破壞模式以及穩定性進行的模擬研究。文獻[17]基于廣東某核電邊坡現場發生的局部滑坡校核巖土力學參數，利用顆粒流方法探討降雨作用下的滑坡的穩定性。此外文獻[18-20]通過離散元的方法對滑坡的建模及參數的標定做了大量的研究，為后續的持續性探索奠定了基礎。

綜上所述，國內外學者對滑坡的研究多集中在巖土體細觀參數對滑坡破壞方式的影響、滑坡穩定性計算等方面，而對滑坡破壞過程中的力學機理分析、滑坡速度變化以及坡體處于時效變形階段的發展趨勢的分析和研究相對較少。因此，以陜西省安康市鎮坪縣柿子樹坪滑坡為研究對象，對其穩定性進行分析，并利用PFC2D對滑坡的變形破壞過程進行模擬，監測滑坡內不同位置的位移、應力指標的變化情況，分析該滑坡在運動過程中呈現出的階段性破壞變形特征，為該滑坡的預防治理提供參考依據。

1 柿子樹坪滑坡概況

1.1 自然地理位置

柿子樹坪滑坡位于陜西省安康市鎮坪縣城關鎮菜村四組，中心點地理位置：北緯31°52′53.08″，東經109°31′33.17″。東南部與湖北省竹溪縣接壤，南部與重慶市巫溪縣、城口縣毗鄰，西北與陜西省平利縣相鄰。S207、S225省道可直通重慶市、湖北省區縣，嵐(皋)—鎮(坪)公路、平(利)—鎮(坪)公路可通往鄰近縣及安康市，交通較為便利(見圖1)。

圖1 滑坡位置圖Fig.1 Location of the landslides

1.2 滑坡的形態、規模特征

根據詳細調查及勘查情況，該滑坡發育于南江河右岸，中低山地貌，微地貌單元屬山前斜坡地段，主要由兩個滑坡及一個局部小滑塌組成，以中部沖溝為界(圖2)，滑坡前緣坡度為35°～40°，中部、后緣坡度為30°～38°，滑坡體前緣以下為居民集中區，該滑坡體上部出現裂縫(圖3)，中下部表層土體滑動變形，村村通道路處形成0.35 m寬的裂縫，最大可見深度達0.5 m(圖4)。

圖2 滑坡形態規模Fig.2 Scale of landslides

圖4 道路旁裂縫Fig.4 Roadside fractures

1#滑坡位于沖溝西側，平面形態較不規則，滑坡后緣可隱約見陡坎，東側緣中上部局部可見基巖，西側緣下部有較大面積出露強風化基巖，前緣土層較薄，滑坡周界較清晰。

2#滑坡位于沖溝東側，平面形態呈簸箕形，滑坡后緣可見陡坎，東側緣為一小型沖溝，西側緣為沖溝，滑坡周界較清晰。

局部滑塌區位于1#滑坡西側，為順層巖質結構，滑塌物質為全風化板巖，較破碎。滑坡基本特征見表1。

表1 柿子樹坪滑坡基本特征Table 1 The basic characteristic of Shizishuping landslide

1.3 滑坡的結構特征

根據鉆探揭露及野外現場鑒別，滑體地層上部為第四系全新統殘坡積含碎石粉質黏土；下伏寒武系水井沱組(∈1s)板巖，工程地質剖面圖如圖5所示，各巖土層的特征分述如下。

圖5 工程地質剖面圖Fig.5 Landslide section

第2層 強風化板巖(∈1s)：黑灰-灰色，細粒結構，層狀及片狀構造，巖芯呈碎塊狀、短柱狀，層理發育，節理及風化裂隙發育。厚度0.50～4.50 m，層底深度5.50～27.70 m，層底高程883.75～908.75 m。

第3層 弱風化板巖(∈1s)：黑灰-灰色，細粒結構，層狀及片狀構造，屬較硬巖，巖芯呈短柱狀、長柱狀，層理發育，節理及風化裂隙較發育，傾向300°～305°，傾角12°～20°。具本地區同類巖石試驗結果，飽和單軸抗壓強度15～35 MPa。此層未揭穿，最大揭露厚度6.5 m。

1.4 滑坡影響因素分析

1.4.1 地形地貌

本區地處中山區，微地貌為山前斜坡地段，坡面主要為林地和耕地。滑坡整體坡度約40°，其中滑坡前緣坡度為35°～40°，中部、后緣坡度為35°～50°。坡體下部因削坡建房改變了坡體的原始地形，形成了高2～10 m、坡度近70°的陡坡，大面積的開挖使坡體局部應力集中，地形因素有利于滑坡的形成。

1.4.2 地層巖性

兩處滑體地層巖性為第四系全新統殘坡積含碎石粉質黏土，碎石含量約占20%，局部地段含量較少，局部滑塌區滑體為全風化板巖，滑床均為強風化-中風化板巖。滑體坡體上部巖土體接受降水時，水體下滲，在有黏性土隔離處形成上層滯水，無黏性土處水體直接下滲到滑床附近。降雨使表層土體含水率增大，下滲的水使滑帶土飽和，重力增加而強度顯著降低，不利于坡體穩定。

1.4.3 水的作用

根據鎮坪縣氣象站1960—2016年降雨資料顯示，雨季一般在5月上旬開始，終止一般在10月上旬，平均為150 d，平均降水量為760 mm，占全年降水量的75%。在暴雨、連陰雨情況下，坡面土體很快飽和，降水快速直接入滲至基巖面形成地下水，表層土體含水率增大，重度增加；下滲的雨水沿基巖面及其上部的碎石土層聚集，造成坡體內水力坡度增大，產生靜水壓力，隨著土體強度降低，進一步降低斜坡的整體穩定性，在重力作用下，易沿產生蠕動變形，是導致滑坡災害發生重要因素。

2 滑坡穩定性分析

近年來受雨季持續性降雨以及當地基礎建設的影響，該滑坡在2014年、2015年發生了多次小規模滑塌，滑坡中部及下部由于道路建設形成局部陡坎。目前滑坡威脅范圍包括：村民33戶156人，通村公路約320 m，直接經濟損失超過1 500萬元，潛在經濟損失過3 000萬。根據野外勘察結果，滑坡滑體為殘坡積含碎石粉質黏土；滑動面位于飽和含碎石粉質黏土底部附近；滑床為強風化板巖，地層界線起伏變化較大，為基巖面松散堆積層滑坡。

通過傳遞系數法對該滑坡的穩定性進行計算，其中滑體抗剪強度參數(見表2)。通過對平面圖(圖6)中6個剖面均進行了天然工況和暴雨工況下穩定性計算，計算結果見表3。

圖6 滑坡區工程地質平面圖Fig.6 Geological plan of landslide project

表2 滑體強度參數Table 2 Mechanical parameters of landslides

表3 柿子樹坪滑坡穩定性計算結果表Table 3 The calculation results about Shizishuping landslides stability

由此可以看出，天然工況下，滑坡基本上處于穩定狀態，但在暴雨工況下，滑坡便處于不穩定或者欠穩定狀態。

3 變形特征分析

由以上的計算結果可以看出，在暴雨工況，滑坡處于不穩定或者欠穩定的狀態，進而反映出降雨是導致該滑坡發生滑動的重要因素。以下將通過PFC2D對飽和狀態下柿子樹坪滑坡的變形破壞特征進行分析研究。

3.1 細觀力學參數標定

首先需要建立細觀力學參數與粉質黏土的宏觀力學參數之間對應關系，故通過雙軸實驗直接對飽和狀態下的含碎石粉質黏土細觀力學參數進行標定。

在PFC中，生成不同粒徑、不同比例的顆粒是可以通過體積分配實現的。由于碎石的粒徑大于20 mm與粉質黏土的粒徑在0.005～0.075 mm差著104的數量級，如果按照粒徑級配曲線生成相對應的顆粒，小顆粒數目過大，計算機運算難度也隨之增大。因此本模型中在生成足夠多的顆粒同時，按照粒徑0.08～0.10 mm(20%)、0.1～0.15 mm(28%)、0.15～0.20 mm(31%)、0.20～2 mm(1%)、2～2.2 mm(20%)進行分配。為避免尺寸效應，建立雙軸實驗的模型尺寸為14 m×28 m，共生成9 359個顆粒。

進行雙軸加載試驗前，需根據實際加載圍壓對模型進行伺服[18-19]，伺服采用的圍壓分別為50、150、300 kPa，通過控制變量法獲取相應的法向接觸強度，切向接觸強度以及摩擦系數。根據不同圍壓下的應力應變關系圖(圖7)，繪制摩爾-庫倫破壞包絡線(圖8)，從而得到相應的滑體的細觀力學參數(見表4)。

表4 滑體顆粒細觀力學參數Table 4 Particle parameters of landslide

圖7 不同圍壓下雙軸壓縮應力應變曲線Fig.7 Stress-strain curve under different confining pressures

圖8 摩爾-庫倫破壞包絡線Fig.8 Mohr-Coulumn damage envelope

3.2 構建滑坡模型

導入滑坡剖面，生成墻體，按照比例對滑坡體進行顆粒填充并賦值，使得設定模型計算至顆粒平衡，從而得到滑坡的初始模型(圖9)。為了更直觀的研究該滑坡的破壞和運動特征，分別在坡體的前、中、后部布置一定數量的測量圓，坐標從左至右分別為：B1(51.750 3，105.274 8)，B2(70.773 3，89.052 9)，B3(89.796 3，73.680 0)，B4(108.819 3，59.147 3)，B5(127.842 3，45.184 9)，B6(146.865 3，30.655 9)，B7(165.888 3，14.568 4)，在滑坡運動過程中通過history對各監測點位移、應力以及速度的變化進行監測。

圖9 滑坡初始模型Fig.9 Initial model of landslide

3.3 滑坡變形過程分析

滑坡不同階段的位移云圖如圖10所示，其中利用Geometry生成的邊界線不同于wall，沒有賦予強度參數指標，只是為了方便觀察滑坡體的變形。

運行至5 000時步[圖10(a)]，從位移顯示結果來看，斜坡土體并沒有產生較大的形變，形變量集中在6～12 cm，基本上以蠕滑變形為主。坡體中部厚堆積層處位移相較于其他部分更大一些，是因為此處堆積層較厚且突出，表層的土體易出現臨空面發生形變。

運行至15 000時步[圖10(b)]，從整體的云圖顯示結果來看，相較于上一階段，滑坡體在自重應力場的作用下，繼續發生著隨時間的形變，表現出“時效變形階段”的特征[20]。此時的坡體形變以中部厚堆積層為界大致分成了兩段，表現形式較為統一，從高處到低處，形變量逐漸增加。坡頂位置開始發育張拉裂縫，斜坡上的土體在重力作用下，對坡腳的擠壓效應也愈發明顯，坡腳處的形變位移較5 000時步已有了明顯變化，并伴有零星的顆粒從邊界溢出。

運行至100 000時步[圖10(c)]，此時坡體隨著“時效變形”的發展，整體的位移云圖也出現了明顯的分段階梯式變形。坡腳處的形變位移最為明顯，已經順坡腳剪出，并不斷后上方發展，已表現出整體滑動的趨勢，開始向“破壞發展階段”[21]過渡。其中坡體中部的堆積層也從開始的表層土體形變發展到深部形變。坡頂處的裂縫寬度逐漸擴大，出現了新的臨空面，此時裂縫處的土體應力釋放，失去支撐力，開始出現顆粒塌落現象。

運行至2 500 000時步[圖10(d)]，滑面貫通，滑體已順滑面滑動，此時坡腳土體的位移量已經達到11 m。

圖10 滑坡位移云圖Fig.10 Landslide displacement map

根據以上特征可知該滑坡在滑移初始階段時，主要產生蠕滑變形，隨著滑坡體在自重作用下變形量的逐漸增加，坡腳的應力不斷累積，坡腳順滑面剪出，并不斷向后方牽引，導致整體滑動，表現出典型的牽引式破壞。

3.4 滑坡運動特征模擬

該滑坡體發展到“時效變形”階段時，便具備了潛在的破壞性。此時對其進行工程防治干預是很有必要的。反之，滑面可能會因持續的變形而逐漸喪失強度，從而加速形變，進入“破壞發展階段”，最終導致邊坡發生失穩破壞。

圖11為各監測點記錄的位移-時步曲線。初期階段(5 000時步之前)，位移變化速率較大，此時是因為切坡建模時給顆粒賦值，原始的平衡狀態被打破，顆粒為了重新達到平衡狀態而快速運動。中期階段(5 000～15 000時步)，7個測量圓記錄到的位移增速明顯變緩，整體上較為同步，印證了滑坡發生初期飽和土體的確存在微小蠕變，但B7測量圓監測到的坡腳處的位移要大于其他部位的監測數據，是因為此時坡腳處收到后方土體的推擠也較為明顯，同時與位移云圖的結果相呼應。后期(15 000時步之后)，B1號測量圓記錄的位移數據明顯在增勢上緩于其測量圓，位移差也越來越大，這是因為斜坡前緣滑動較快，使后緣產生張拉裂縫，隨著裂縫逐漸擴張，后緣失去了前部的牽引作用，僅靠重力驅動向下滑塌。因此B1號測點區域土體與前部的位移差越來越大。

圖11 位移-時步曲線Fig.11 Displacement-duration history curve

圖12為各監測點記錄的應力-時步曲線當模型開始加載，平衡被打破的時候，斜坡土體開始運動，初始階段，所有測量圓記錄到的應力隨時間以指數形式增大，且在100 000時步前后都達到了初始階段的應力峰值。由此可見，在變形初期，滑體蠕變量的增加，各監測部位所受的壓應力逐漸增加，由于自重驅動，應力不斷累積，達到峰值出現應力回落現象，此時坡腳開始出現零星滑塌現象。由B7號測量圓記錄的應力數據來看，坡腳滑面在并未貫通的情況下，斜坡土體在自重作用下不斷發生形變，持續對坡腳產生作用力，坡腳應力集中較為明顯且較其他測量圓測量的數據大，當應力積累到一定程度時，坡腳發生剪切破壞，滑體沿著滑面滑出，應力出現回落。B1號測量圓記錄的應力數據是最早出現應力回落現象的，是因滑體在持續變形后，坡體后緣出現張拉裂縫，應力釋放。B4號測量圓在運行至110 000步時，應力持續增大，此時是因坡體后緣變形加快，坡體中部原后堆積層的的鎖固作用，暫時減緩了后方土體的滑移，此時中部堆積厚度進一步增大，土體的壓應力產生的效果明顯，在應力-時步曲線上表現出其他測量圓記錄的數據都出現應力回落現象后，此處應力在一小段時間內持續增大。隨后由于中部鎖骨段被剪斷，發生滑移，使得該處測量圓記錄的應力之后也開始逐漸降低。

圖12 應力-時步關系曲線Fig.12 Stress-duration history curve

從整體運動位移及應力變化趨勢上看，該滑坡體應力和位移總體上呈現出從前到后，從下到上逐漸遞減的變化模式。應力上因前端的滑體變形較快，導致后端出現拉張裂縫，應力釋放，最先出現應力回落現象。位移上由于滑動過程中后方的土體受到前部土體阻擋，前端形變位移較后方位移大。直至滑坡停止，各處的位移將達到最大值。此過程坡體也體現出前緣沿接觸面滑移、坡頂出現張拉裂縫、中部鎖固段剪斷的階段性變形特征[21]。

4 結論

通過現場調查以及勘察結果對滑坡的穩定性進行分析，并結合顆粒流PFC2D對陜西省安康市鎮坪縣柿子樹坪滑坡的變形破壞特征進行了模擬研究，得出以下結論。

(1)天然工況下，滑坡基本上處于穩定狀態，但在暴雨工況下，滑坡的穩定系數均小于1.0，整體處于不穩定的狀態。

(2)該滑坡在滑體飽和狀態下，首先發生蠕滑變形，隨著滑坡體在自重作用下變形量的逐漸增加，坡腳的應力不斷累積，坡腳順滑面剪出，并不斷向后方牽引，導致整體滑動，表現出典型的牽引式破壞。

(3)從整體運動位移及應力變化趨勢上看，該滑坡體應力和位移總體上呈現出從前到后從下到上逐漸遞減的運動特征。也呈現出前緣沿接觸面發生滑移、坡頂出現張拉裂縫、中部鎖固段剪斷的階段式的變形特征。

(4)模擬結果與現場推斷情況較吻合，說明PFC2D適用于該滑坡變形破壞特征的分析研究，對該滑坡的預防及治理具有一定的參考價值。