基于PFC3D的魚(yú)鰍坡滑坡運(yùn)動(dòng)過(guò)程分析

張家勇，鄒銀先，楊大山

（貴州省地質(zhì)環(huán)境監(jiān)測(cè)院，貴州 貴陽(yáng) 550001）

0 引言

貴州省地處西南腹地，是無(wú)平原支撐的省份，地形切割劇烈，地質(zhì)環(huán)境脆弱，地質(zhì)災(zāi)害頻發(fā)[1?3]。2016年5月23日，貴州省開(kāi)陽(yáng)縣龍崗鎮(zhèn)大石板村魚(yú)鰍坡組突發(fā)滑坡災(zāi)害，造成下方1 戶居民2 人死亡，5 間房屋毀壞，10 余畝耕地?fù)p毀，直接經(jīng)濟(jì)損失90 萬(wàn)元（圖1）。魚(yú)鰍坡滑坡為堆積層失穩(wěn)，是貴州省內(nèi)典型突發(fā)性災(zāi)害類型之一。貴州及我國(guó)西部山區(qū)常見(jiàn)堆積體，在自然、人為因素影響下，容易誘發(fā)滑坡地質(zhì)災(zāi)害，給人民的生命財(cái)產(chǎn)安全造成了嚴(yán)重的威脅。黃潤(rùn)秋[4]通過(guò)大量的研究表明：查清滑坡變形破壞的地質(zhì)力學(xué)模式是滑坡地質(zhì)災(zāi)害防治的基礎(chǔ)所在。

圖1 滑坡發(fā)生前后影像圖Fig.1 Image before and after landslide

當(dāng)前，針對(duì)滑坡機(jī)理分析主要是借助地質(zhì)力學(xué)方法開(kāi)展定性分析，判斷破壞模式。近年來(lái)，隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值模擬方法在滑坡領(lǐng)域得到了廣泛應(yīng)用。目前的數(shù)值模擬方法中，主要包括確定性分析和非確定性分析兩類方法，而確定性分析又可分為連續(xù)介質(zhì)分析、非連續(xù)介質(zhì)分析方法以及近年來(lái)李世海等[5]提出的基于連續(xù)介質(zhì)力學(xué)的離散元方法(CDEM)。其中連續(xù)介質(zhì)數(shù)值分析方法有有限單元法、邊界元法、有限差分法等，非連續(xù)介質(zhì)分析方法有塊體離散元法、顆粒離散元法、關(guān)鍵塊體理論、不連續(xù)變形分析(DDA 法) 等。CDEM 方法將有限元與離散元進(jìn)行耦合，在塊體內(nèi)部進(jìn)行有限元計(jì)算，在塊體邊界進(jìn)行離散元計(jì)算，通過(guò)塊體內(nèi)部及塊體邊界的斷裂，不僅可以模擬材料在連續(xù)狀態(tài)下及非連續(xù)狀態(tài)下的變形、運(yùn)動(dòng)特性，更可以實(shí)現(xiàn)材料由連續(xù)體到非連續(xù)體的漸進(jìn)破壞過(guò)程。

離散單元法，它的基本思想是把滑體介質(zhì)看作由一系列離散的、獨(dú)立運(yùn)動(dòng)的單元所組成。20 世紀(jì)90年代以后，離散單元法發(fā)展十分迅速，基于不規(guī)則形狀塊體單元，如UDEC、3DEC 等軟件，以及基于球形和圓盤形散體單元的軟件，如PFC、MatDEM 等，得到了不斷的完善，使其在許多的領(lǐng)域都得到了廣泛的應(yīng)用。顆粒流法（PFC）作為離散單元法的一種，通過(guò)圓形離散單元（Ball）來(lái)模擬顆粒介質(zhì)的運(yùn)動(dòng)及其相互作用，克服了傳統(tǒng)連續(xù)介質(zhì)力學(xué)的宏觀連續(xù)性假設(shè)，能真實(shí)直觀地模擬滑坡的失穩(wěn)破壞過(guò)程[6]。目前一些學(xué)者逐步將其應(yīng)用于滑坡破壞運(yùn)動(dòng)分析之中。吳劍[7]采用PFC 對(duì)滑帶剪切過(guò)程進(jìn)行了模擬研究，張龍等[8]、施鳳[9]、李冬冬[10]、王畯才[11]采用PFC3D對(duì)多個(gè)高速遠(yuǎn)程滑坡的運(yùn)動(dòng)過(guò)程進(jìn)行了模擬，深入討論了位移、速度、堆積形態(tài)等各方面的因素。杜永彬[12]、王宇等[13]、胡江春等[14]通過(guò)采用PFC2D對(duì)滑坡形成演化過(guò)程進(jìn)行模擬，得到了良好的效果。趙洲等[15]、陳達(dá)等[16]、閆曉娟等[17]、汪華安等[18]采用PFC2D對(duì)堆積層滑坡形成演化過(guò)程進(jìn)行模擬，并對(duì)堆積層滑坡的機(jī)理進(jìn)行了總結(jié)[19?21]。

因此，本文利用顆粒流程序PFC3D模擬研究魚(yú)鰍坡滑坡的破壞運(yùn)動(dòng)過(guò)程，選用Ball-Wall 建模方法，引入顆粒流(PFC3D)程序中平行黏結(jié)模型，進(jìn)行滑坡模型的建立，對(duì)滑坡不同關(guān)鍵部位顆粒進(jìn)行位移、速度監(jiān)測(cè)，模擬分析魚(yú)鰍坡滑坡破壞運(yùn)動(dòng)過(guò)程。研究成果可為對(duì)該類滑坡影響范圍預(yù)測(cè)，以及工程措施的制定具有一定的參考意義。

1 滑坡概況

魚(yú)鰍坡滑坡位于開(kāi)陽(yáng)縣龍崗鎮(zhèn)大石板村魚(yú)鰍坡組。研究區(qū)屬亞熱帶季風(fēng)氣候，冬季干冷少雨，夏季受季風(fēng)影響，暖濕多雨。四季降雨分布不均，且降雨多集中在汛期4—9月，占全年降雨量的75.8%，滑坡發(fā)生前10 天的累計(jì)降雨量達(dá)到了70 mm?；滤巺^(qū)域?qū)傩逼碌孛?，總體地勢(shì)為東北高，西南低，相對(duì)高差約210 m，平均地形坡度25°，地形起伏變化較大?；聟^(qū)地層巖性主要為第四系殘坡積層（Qel+dl）含碎石黏土和寒武系下統(tǒng)金頂山組(∈1j)粉砂巖、泥巖。魚(yú)鰍坡滑坡地處大翁林背斜南東翼，靠近背斜核部，滑坡前緣靠近一條北西-南東向平移斷層，左側(cè)緣靠近一條北東-南西向的正斷層。地下水分為碎屑巖裂隙水與松散層孔隙水，根據(jù)資料顯示該區(qū)地下水主要富存在第四系松散堆積層中，其透水性及含水性較好。地下水以大氣降水補(bǔ)給為主，地表水補(bǔ)給次之。

魚(yú)鰍坡滑坡在平面上的投影大致呈上小下大的塔式形狀，滑坡在剖面上呈現(xiàn)出近似直線狀，主滑方向244°，滑坡斜長(zhǎng)約140 m，平均橫寬約60 m，平均厚度約4 m，總體積約40 000 m3（圖2）?；峦馏w為寒武系牛蹄塘-金頂山組地層風(fēng)化形成的殘坡積土，碎石含量較少，碎石粒徑較小（圖3）。

圖2 魚(yú)鰍坡滑坡平面圖Fig.2 Plan of Yuqiupo Landslide

圖3 滑坡剖面圖（I?I′）Fig.3 Section of landslide (I?I′)

魚(yú)鰍坡滑坡的變形破壞是其獨(dú)特的內(nèi)外因共同作用的結(jié)果，坡體臨空條件，上部疏松下部致密的巖土體組合，不良地表排水條件等為內(nèi)因，短期強(qiáng)降雨入滲為外因。降雨為魚(yú)鰍坡滑坡的直接誘發(fā)因素，斜坡變形破壞模式為上部的土體沿基巖面先出現(xiàn)變形，對(duì)下部滑體產(chǎn)生推擠作用，直至坡腳產(chǎn)生剪切破壞，并向上牽引發(fā)展，導(dǎo)致整體滑動(dòng)。

2 滑坡破壞過(guò)程數(shù)值模擬

2.1 計(jì)算模型的建立

本文根據(jù)滑前1∶1 萬(wàn)地形圖、滑后無(wú)人機(jī)航測(cè)DEM，生成滑體和滑床的幾何邊界模型，并導(dǎo)出為PFC3D能識(shí)別的STL 格式，通過(guò)“Geometry import”命令導(dǎo)入PFC 中作為兩個(gè)“Geometry”。選擇Ball-Wall 模型建模，其中滑床的“Geometry”作為Wall 模型。然后通過(guò)有限步循環(huán)計(jì)算消除顆粒間的內(nèi)力，得到滑坡最終計(jì)算模型（圖4）。

圖4 滑坡數(shù)值計(jì)算模型Fig.4 Numerical calculation model of landslide

該滑坡模型長(zhǎng)230 m，寬170 m，高115 m，顆粒最小半徑0.25 m，最大半徑0.5 m，顆?？倲?shù)為2 447 個(gè)。

為了從整體和局部?jī)煞矫婀餐盐栈w的運(yùn)動(dòng)特征，此次模擬在滑坡前緣、滑坡中部、滑坡后緣共設(shè)置了6 個(gè)監(jiān)測(cè)點(diǎn)，分別用于監(jiān)控滑坡過(guò)程中滑體的速度、位移的變化（圖4）。

2.2 細(xì)觀力學(xué)參數(shù)的確定

PFC 滑坡模擬參數(shù)主要包括本征參數(shù)和接觸參數(shù)，其中本征參數(shù)主要包括剪切模量、泊松比和密度，可通過(guò)室內(nèi)巖土體物理力學(xué)試驗(yàn)和工程類比確定，顆粒單元間的接觸參數(shù)主要通過(guò)以下兩個(gè)途徑獲取[12]：（1）通過(guò)PFC 進(jìn)行參數(shù)標(biāo)定確定物理力學(xué)參數(shù)。（2）通過(guò)現(xiàn)場(chǎng)參數(shù)反演，即對(duì)所需微觀參數(shù)進(jìn)行賦值使顆粒流的運(yùn)動(dòng)狀態(tài)與實(shí)際狀態(tài)基本一致，可認(rèn)為參數(shù)為合理值。本文主要選取第二種現(xiàn)場(chǎng)反演的方法確定顆粒單元的接觸參數(shù)，具體取值見(jiàn)表1。

表1 滑坡巖土體顆粒細(xì)觀參數(shù)Table 1 Mesoscopic parameters of rock and soil particles in landslide

2.3 破壞過(guò)程模擬及分析

模擬計(jì)算到1 s 時(shí)，坡體出現(xiàn)變形，邊坡整體上位移較小，滑坡處于蠕滑變形階段，見(jiàn)圖5（a）。3 s 時(shí)，坡體變形進(jìn)一步增大，且坡體中后部顆粒的位移較中前部顆粒大，見(jiàn)圖5（b）。到達(dá)5 s 時(shí)，滑坡前緣產(chǎn)生剪切破壞，坡腳顆粒開(kāi)始剪出，見(jiàn)圖5（c）。7 s 時(shí)，滑坡滑動(dòng)帶已基本貫通，在坡體壓力作用下滑坡開(kāi)始整體滑移，見(jiàn)圖5（d）。到達(dá)8 s，滑坡整體的速度差進(jìn)一步減小，滑坡處于整體加速階段，見(jiàn)圖5（e）。到達(dá)10 s 時(shí)，速度峰值12.4 m/s，此時(shí)，近一半的滑體從滑源區(qū)滑出，見(jiàn)圖5（f）。到達(dá)15 s時(shí)，由于地形原因，前緣顆粒的速度逐漸減小，見(jiàn)圖5（g）。30 s 時(shí)，滑坡速度進(jìn)一步減小，前緣的滑坡已基本停止運(yùn)動(dòng)，見(jiàn)圖5（h）。到50 s 時(shí),滑坡停止運(yùn)動(dòng)并堆積于坡腳平緩處，見(jiàn)圖5（i）。

圖5 不同時(shí)步速度云圖（單位：m）Fig.5 Velocity cloud map of different steps (unit: m)

綜上可知，該滑坡在破壞初始階段以蠕滑變形為主，隨著變形量的增加，滑坡體不斷擠壓坡腳，滑坡巖土體到達(dá)應(yīng)力平衡極限，坡腳產(chǎn)生剪切破壞，并向上牽引發(fā)展，滑坡發(fā)生整體滑動(dòng)，整個(gè)滑坡歷時(shí)50 s，整體運(yùn)動(dòng)距離約為80 m，滑坡運(yùn)動(dòng)距離較短。整個(gè)滑坡破壞模式表現(xiàn)為牽引式破壞。模擬結(jié)果與實(shí)際變形情況基本一致。

2.4 運(yùn)動(dòng)特征分析

通過(guò)觀察速度?時(shí)間曲線（圖6），整體上坡體各部分到達(dá)速度峰值的時(shí)刻基本一致，均在10 s 到達(dá)峰值，各點(diǎn)速度峰值介于5～11.6 m/s，速度峰值最高顆粒是位于滑體前緣的3 號(hào)顆粒，峰值11.6 m/s。速度峰值最低顆粒是位于滑體后緣的6 號(hào)顆粒，峰值約7 m/s。6 個(gè)顆粒速度時(shí)程曲線整體上變化趨勢(shì)基本一致。

圖6 速度時(shí)程曲線Fig.6 Time-history curve of velocity

通過(guò)觀察位移–時(shí)間曲線（圖7），監(jiān)測(cè)顆粒均在45～50 s 達(dá)到位移峰值。顆粒最大位移94 m，為滑體前緣2 號(hào)顆粒；最小位移60 m，為坡后緣6 號(hào)顆粒。總的趨勢(shì)是滑坡前部顆粒較后部顆粒位移大，呈現(xiàn)出牽引式破壞特征。

圖7 位移時(shí)程曲線Fig.7 Displacement time-history curve

通過(guò)觀察平均速度時(shí)間曲線與平均位移時(shí)間曲線(圖8)，結(jié)果表明：0～7 s 時(shí)平均速度呈近似線性增大，約在10 s 時(shí)達(dá)到速度峰值7 m/s，到達(dá)峰值后，平均速度同樣呈近似線性減小；30 s 時(shí)以后降速更緩，漸趨于穩(wěn)定，但仍有較小速度，約0.3 m/s；到50 s 后，平均速度趨于0。

圖8 平均速度、位移時(shí)程曲線Fig.8 Average velocity and displacement time-history curves

平均位移時(shí)間曲線與平均速度時(shí)間曲線分析各階段速度變化情況一致。0～7 s 平均位移呈逐漸增長(zhǎng)，斜率逐漸增大，為加速階段；7～10 s 近似呈線性增大，10～30 s 時(shí)步斜率逐漸減小，為減速階段；50 s 后平均位移到達(dá)峰值80 m。

3 結(jié)論

（1）魚(yú)鰍坡滑坡是在內(nèi)外因共同作用下的結(jié)果，滑體物質(zhì)為碎石土，物理力學(xué)性質(zhì)差，降雨為魚(yú)鰍坡滑坡的直接誘發(fā)因素。

（2）PFC 數(shù)值模擬結(jié)果表明：該滑坡在破壞初始階段以蠕滑變形為主，滑坡體不斷擠壓坡腳，坡腳產(chǎn)生剪切破壞，并向上牽引發(fā)展，滑坡發(fā)生整體滑動(dòng)。整個(gè)滑坡破壞模式表現(xiàn)為牽引式破壞。模擬結(jié)果與實(shí)際變形情況基本一致。

（3）該滑坡歷時(shí)約50 s，0～10 s 為整體加速階段，滑坡最高時(shí)速12.4 m/s，11～50 s 為整體減速階段。整體位移約80 m，滑坡前部位移要比后部位移大，整體速度峰值介7～12 m/s 之間，速度和位移變化總趨勢(shì)為初始高程越低，達(dá)到的最大速度和位移越大。

（4）通過(guò)數(shù)值模擬發(fā)現(xiàn)，PFC3D軟件用于模擬堆積層滑坡運(yùn)動(dòng)過(guò)程以及形態(tài)具有較好的適用性，對(duì)此類滑坡的預(yù)防和后期治理具有一定的參考價(jià)值。