滑坡漸進破壞運動過程的顆粒流仿真模擬

王 宇，李 曉，王聲星，侯文詩

滑坡漸進破壞運動過程的顆粒流仿真模擬

王 宇1，2，李 曉1，王聲星1，2，侯文詩1，2

（1.中國科學院地質與地球物理研究所工程地質力學重點實驗室，北京 100029；2.中國科學院大學，北京 100049）

采用顆粒流離散元模擬滑坡的變形破壞全過程，不需要假定滑動面的位置和形狀，顆粒根據所受到的接觸力調整其位置，最終從抗剪強度最弱面發生剪切破壞。白河滑坡是受強降雨影響誘發的上部順層、下部微切層的大型古基巖滑坡，滑坡的發生是其獨特內外因共同作用的結果。據現場調查資料，總結了滑坡的成因與形成機制。運用基于離散元法的顆粒流（PFC）程序，引入平行粘結模型，通過數值模擬，確定滑帶土細觀參數與宏觀力學性質的關系；據此建立滑坡模型；然后，運用顆粒流離散元法對滑坡運動漸進全過程進行數值模擬，對滑坡不同關鍵部位進行位移、孔隙率變化及應變監測，闡明其漸進發展過程，明確其時空演化規律。模擬結果表明：上部的土體沿基巖面先出現開裂，導致中部滑體產生剪脹，對下部古滑體產生推擠作用，導致新老崩積堆積區及白河中學食堂地區的失穩，形成整體滑動。通過對白河滑坡發生→發展→破壞的動態演化模擬，可為工程決策提供理論依據。

滑坡；顆粒流（PFC）；數值模擬；漸進破壞；平行粘結模型

1 研究背景

滑坡破壞通常并不是一蹴而就的，而是一個巖土體隨時間而緩慢變形和突發性崩塌的過程，該過程具有時空漸進性。邊坡工程界對這種現象引起重視是在20世紀60年代末，由Skempton［1］首次提出了邊坡漸進破壞的概念，并指出邊坡失穩現象往往并不是在同一時間滑面上各點一起達到極限狀態的，最可能的模式是從一個局部開始破裂，然后向其他各部分發展，最后形成整體滑動。由于滑坡漸進破壞過程的復雜性，較真實地反演邊坡漸進從局部發生到貫穿整個邊坡的過程相對比較困難，國內外學者多是采用數值模擬手段［2－6］，借助先進軟件如FLAC，DDA及UDEC進行模擬分析。然而，滑坡體作為一種復雜的地質體（屬于非連續介質），坡體的不同部位其力學性質、應力狀態、位移的規律等是不同的。應用連續介質分析程序（如FLAC）進行模擬分析并不是最優的選擇；尤其在對主要成分為土石混合體的堆積層滑坡進行模擬分析時，土石混合體由于具有高度非均質、非連續、非線性等特點，在力學特性上表現出顯著的獨特性，導致堆積層滑坡與土質或巖質滑坡有不同的演化模式與變形破壞特征，應用塊體離散元程序（UDEC，DDA）進行分析準確性差，甚至不能準確地反映其變形演化過程。

鑒于上述情況，本文采用二維顆粒流程序（Particle Follow Code PFC2D）數值模擬新技術，研究土石混合體滑坡的變形過程應更為恰當。PFC2D程序可用于顆粒團粒體的穩定、變形及本構關系等力學性態的分析，專門用于模擬固體力學大變形問題。它通過圓形（或異型）離散單元來模擬顆粒介質的運動及其相互作用。本文首先采用顆粒PFC2D程序，通過雙軸壓縮數值試驗，確定滑帶土細觀參數與宏觀力學性質的關系，得到給定細觀力學參數試樣的宏觀力學反映；然后，引入平行粘結模型，通過數值試驗建立滑坡地質模型，模擬分析白河滑坡發生、發展、漸進破壞的全過程。

2 顆粒流程序PFC2D簡介

PFC2D程序是基于離散單元法來模擬圓形顆粒介質的運動及其相互作用的。土的結構可視為一個由單粒、集粒或凝塊等骨架單元共同形成的結構體系，這些體系之間的黏結關系，即嵌固咬合狀態、組合特征及散粒體的強度性質、大小等影響著材料的應力和變形的非線性關系。這些離散的顆粒僅在接觸部分受力，當接觸點受到的作用力大于接觸強度時，這些顆粒相互分離，使模型物體發生變形或位移［7－9］。PFC顆粒流離散元是位移分析法，顆粒流理論在整個計算循環過程中，交替應用力－位移定律和牛頓運動定律，通過力－位移定律更新接觸部分的接觸力，通過運動定律，更新顆粒與墻邊界的位置，構成顆粒之間的新接觸。

3 平行粘結接觸本構模型

在PFC2D中，材料的本構特性是通過接觸本構模型來模擬的。顆粒的接觸本構模型有：①接觸剛度模型；②滑動模型；③粘結模型。接觸剛度模型提供了接觸力和相對位移的彈性關系，滑動模型則強調切向和法向接觸力使得接觸顆粒可以發生相對移動，而粘結模型是限制總的切向和法向力使得在粘結強度范圍內發生接觸。PFC2D提供了2種粘結模型，即：接觸粘結模型（Contact-bond model）和平行粘結模型（Parallel-bond model）。接觸粘結認為粘結只發生在接觸點很小范圍內，而平行粘結發生在接觸顆粒間圓形或方形有限范圍內。接觸粘結只能傳遞力，而平行粘結同時能傳遞力矩［10－11］。

平行粘結模型是由法向剛度ˉkn、切向剛度ˉks、法向強度ˉσc、切向強度ˉτc和粘結半徑ˉR這5個參數定義的。與平行粘結相對應的總接觸力和力矩用ˉFi和ˉM3表示，根據習慣，總的接觸力和力矩表示平行粘結在顆粒B上的作用，如圖1所示。可以將總的接觸力沿接觸面分解為切向分量和法向分量，即

圖1 顆粒間平行粘結模型Fig.1 Parallel-bond model for particles

當粘結形成時，ˉFi和ˉM3均初始化為零，以后在接觸處由位移增量和旋轉增量引起的彈性力和力矩的增量將迭加在當前值中。在一個時步Δt內，轉化為彈性力和彈性力矩增量，如式（3）和式（4）所示。

式中：ˉkn，ˉks分別為平行粘結的法向剛度和切向剛度；ni為接觸點法向矢量；A為平行粘結截面面積；I為接觸面相對于過接觸點沿Δθ3方向軸的慣性矩；V為粘結速度。有了力和力矩的初值及增量，得到新的力和力矩力為

根據梁理論分析，在平行粘結周圍分布的最大拉應力和最大剪應力分別為

當最大張應力和最大剪應力分別超過法向與切向粘結強度時，平行粘結發生破壞。相應地，在滑坡破壞過程中表現為張拉破壞和剪切破壞。

4 工程實例分析

4.1 滑坡工程地質概況

白河滑坡位于河南省嵩縣白河鄉白河街南山北坡上，為一大型古基巖滑坡。在上世紀20年代曾發生過滑動，隨后在20世紀70—80年代又相繼出現裂隙與變形。受2010年7月3日至24日強降雨的影響，斜坡中下部古滑坡部分復活［12］。滑坡體平面總體上呈圈椅狀，后緣呈弧形，側緣呈扇形展開，前緣呈弧形。滑坡后緣高程約655 m，前緣剪出口高程約551 m，滑坡前緣寬200 m（不包括影響區），后緣寬600 m。滑坡平均坡度約30°，滑動主滑方向45°，縱向長度約710 m，面積40 426 m2，厚度7～25 m，按平均厚度15m計算，滑坡體積約580 000× 104m2。滑坡區出露基巖為上元古界寬坪群角山片巖，上部覆蓋層地層巖性為新生界第四系上更新統粉質黏土夾碎塊石。邊坡工程地質剖面圖見圖2，室內土工試驗結果見表1。

4.2 滑坡影響因素及成因機制分析

圖2 白河滑坡地質縱剖面示意圖Fig.2 Geological profile of the Baihe landslide

白河滑坡的發生是其獨特內外因共同作用的結果，坡體臨空、巖層產狀上陡下緩、順向坡地質結構、巖性軟硬相間、棋盤格式的裂隙網絡、不良的地表排水條件等為內因，長期持續的強降雨入滲則是其外部誘因。斜坡變形破壞模式屬于孔隙水壓力誘發的平推式地質模式，強降雨是滑坡的觸發因素。白河滑坡具有“推落式”力學性質與特點，表現為：

（1）滑坡變形最先發生于坡體的上部，而后逐漸由上向下傳遞或擴展，最終導致古滑坡的復活。根據監測資料，后緣裂隙不斷擴展。

（2）坡體發育于“靠背椅狀”順向坡地質結構。滑坡的“動力”主要源自坡體下部“阻滑段”自重“喪失”、抗滑力“損失”及中上部“下滑段”的推動。

（3）從滑坡整體來看，滑坡區存在層間錯動帶和前緣緩傾角裂隙性斷層，這些因素會使滑坡區形成貫穿滑動面，持續的降雨不斷入滲滑體內，滲透水壓力不斷增大，對坡體起平推作用。

4.3 數值試樣試驗

由于沒有具體的數學表達式來建立巖土宏觀力學參數和細觀力學參數間的關系，而宏觀參數又不能直接應用于PFC2D中，而雙軸壓縮數值試驗可以得到給定細觀力學參數試樣的宏觀力學反映。降雨是白河滑坡的主要誘因，因滑坡由土石混合體構成，混合體空隙大，降雨使滑坡自重增加；另外，滑帶土在飽水作用下發生軟化，抗剪強度降低，使坡體不斷產生蠕滑－拉裂變形至失穩。因此，獲得滑帶土飽和殘剪抗剪強度指標（表1）相當的細觀力學參數是模擬分析工作的關鍵。建立雙軸壓縮試驗模型，模型尺寸為40 mm×20 mm，見圖3。在雙軸數值模擬試驗中，設定顆粒半徑由不同半徑的顆粒單元組成。半徑R的分布采用從Rmin到Rmax的高斯分布，輸入參數接觸強度分別取通過各種比較模擬，反復調試顆粒微觀參數，使其能夠反映巖土體的宏觀力學參數，與室內滑帶土飽水殘剪試驗相吻合。以粉質黏土數值試驗為例，PFC2D模型輸入參數見表2。

圖3 雙軸壓縮數值試驗模型Fig.3 The biaxial compression testmodel

表2 PFC2D模型輸入參數Table 2 PFC2Dmodel parameters of landslide

根據表2參數，由程序模擬可以得到典型粉質黏土在輸入圍壓下的全應力－應變曲線，典型應力－應變曲線如圖4所示，與室內試驗相吻合。

5 滑坡漸進破壞數值模擬分析

5.1 計算模型的建立

采用由數值試樣試驗得到的細觀參數（見表3）建立計算模型（圖5），顆粒間采用顆粒元程序中的平行粘結模型，該模型不但能抗拉、抗剪，還能承受彎矩，能很好地模擬具有黏聚力的巖土材料。強、弱風化角閃片巖層的細觀參數通過工程類比及試算獲得。不同顆粒間的接觸面可以看作是巖體中的節理。邊界墻體法向剛度取kn＝1.0×108N／m，切向剛度ks＝1.0×108N／m，摩擦系數為0.5，粉質黏土天然重度γ＝19.8 kN／m3，讓顆粒在自重作用下計算達到平衡狀態，從而模擬得到滑坡的初始應力場（見圖6）。黑色線代表顆粒間壓應力，線愈粗，代表壓應力愈大，最大應力位于滑坡底部，經試算和實際相況基本一致。

表1 土工試驗參數Table 1 Physical-mechanical parameters of soil specimens

圖4 試樣全應力－應變曲線Fig.4 Complete stress-strain curve of soil specimens

表3 巖土體微觀參數取值Tab le 3 M icro-parameters of rock soilmass

圖5 數值計算模型圖Fig.5 Numerical calculation model

平均不平衡接觸力變化曲線見圖7，從圖中可以看出，經過迭代運算后，體系最大不平衡接觸力隨迭代計算的運行，逐漸逼近于0，表明體系最終達到了力平衡狀態，滑坡的初始應力場已形成。接下來就可以用表3中的與滑帶土飽和殘剪抗剪強度指標相當的細觀力學參數進行漸進破壞的模擬分析。

圖6 滑坡初始應力場Fig.6 Initial stress field of the landslide

圖7 平均不平衡接觸力Fig.7 Curve of average unbalance force

5.2 模擬結果及分析

采用PFC2D程序對滑坡的變形和位移進行計算，對滑坡發生、發展的全過程進行仿真模擬。在模擬過程中，分別對坡頂、頂中和頂底進行位移變形監測、孔隙度及最大剪應力－剪應變監測，以追蹤滑坡漸進破壞的全過程。

模擬計算到5 000時步時，如圖8（a），坡頂面出現拉張拉裂隙，并且裂隙有不斷擴張的趨勢，說明坡體開始發生破壞。到達10 000時步時，如圖8（b），因坡腳位移增大，導致上部裂隙擴大，滑移加劇，顆粒間的接觸變得不規則，滑坡中部出現鼓脹變形，坡體中出現剪切滑動，滑坡上部顆粒翻滾滑落。到達20 000時步時，如圖8（c）堆置體中部出現剪切滑動，底部顆粒受上部重壓開始錯位，整個滑體間顆粒粘結破壞效應進一步加劇，但滑坡后緣滑移長度并未發生明顯變化，說明滑坡處于蠕滑變形階段，由于坡體下部的鎖固作用，中部鼓脹變形加劇。到達31 000時步時，裂隙遍布滑坡體，分布規律是滑坡前緣和滑面裂隙密集，滑體內中部裂隙稀疏，滑體產生了大幅的滑動變形，后緣與原坡面線相比，產生明顯錯落，見圖8（d）。

圖8 滑坡隨時步變形情況Fig.8 The deform ation at different time-steps

隨計算時步的增加，滑坡變形不斷增大，滑坡破壞始于后緣的拉張裂隙，裂隙擴大在重力作用下推動下部滑體滑移，滑坡變形經歷了后緣拉裂、蠕滑、變形加劇的過程。圖9為滑坡剪切帶貫通時的位移圖，滑坡沿基巖面發生向下的滑移破壞。

圖9 滑坡位移矢量圖Fig.9 Displacement vectors of the landslide

模擬分析時顆粒間接觸模型采用顆粒元程序中的平行粘結模型，該模型不但能抗拉、抗剪，還能承受彎矩，可以很好地反映顆粒的變形情況。滑坡剪切帶漸進展過程中平行粘結力的變化情況，模擬時步較小時，坡體內平行粘結力以拉為主，剪切力較小，時步增加時，平行粘結力發展為以剪力為主，隨著計算時步的為斷加大，剪切滑動面向貫通趨勢發展。圖10為剪切帶貫通時，由3個監測圓得到的應變隨時步變化曲線圖。在滑坡漸進破壞過程中，應變由受拉變為受剪，直到剪切帶貫通。坡體中應變變化速率依次為：坡中、坡底和坡頂。由于滑坡上部土體的擠壓，坡中應變變化劇烈，顆粒旋轉、平動不斷調整位置，以致剪切帶貫通，而后緣坡頂處因顆粒不斷滑落，顆間接觸松散，應變發展不明顯。圖11為模型不同位置的監測圓中孔隙率變化過程。在初期，各個監測圓中的孔隙率都接近初始給定值。模擬初期，受重力作用，顆粒擠密，3個監測圓內的孔隙率變小，隨著滑體上部土體顆粒的滑動，材料發生剪脹效應，3個監測圓內的孔隙率值整體呈緩慢增長的趨勢，由于坡底土體鎖固作用，土體被繼續壓密。剪切帶貫通后，剪脹作用消失，孔隙率變小。其中，坡體中部位置孔隙率增長較快，伴有顆粒彈性變形、磨粒、擠密的產生與遷移，顆粒間變化更為復雜，說明滑坡土石混合體在該處變形較為劇烈。相對于坡中，坡頂和坡底處顆粒孔隙率增長幅度相對較小，變形活動相對緩和。圖12為坡頂和坡中顆粒滑移在同一個坐標系中的記錄，坡頂顆粒滑落速度明顯大于頂中，由于坡面臨空面較陡，裂隙的擴張導致坡頂顆粒迅速下滑，由于坡中顆粒的剪脹作用使此處顆粒滑移速度較慢，由于顆粒間的旋轉、剪切作用使位移變小，趨于穩定。

圖10 監測圓剪應變速率曲線Fig.10 Curves ofmonitored round shear strain rate

圖11 坡體不同位置孔隙率變化曲線Fig.11 Curves of porosity change in different positions

圖12 監測點豎向位移曲線Fig.12 M onitored curves of vertical displacement

6 結論與討論

利用基于顆粒流理論的PFC2D程序對其白河滑坡進行數值模擬，再現了滑坡土石混合體滑移、變形、破壞的整個過程，得出了以下結論：

（1）通過顆粒元平行粘結模型，采用數值試驗的手段，很好地模擬了土石混合體堆積體滑坡的漸進性破壞過程，對滑坡的變形特點及成因機制有了更深的認識，這一結果為滑坡后續加固處治方案的確定提供了理論依據。

（2）滑坡的直接誘發因素是降雨，滑坡土體強度不足造成的滑坡破壞首先從坡頂處產生裂隙、繼而發生小型的滑塌開始，而后演變成大的滑坡，破壞形式表現為上部拉裂，中部剪切，底部擠壓破壞。

（3）對于含有土石混合體的堆積體滑坡來說，滑坡具有明顯的土－石混合介質的特征，具有高度非均質、非連續、非線性等特點，在力學特性上表現出顯著的獨特性，導致堆積層滑坡與土質或巖質滑坡有不同的發育模式與變形破壞特征。下一步的工作，因在野外及室內試驗的基礎上，確定土石混合比，建立更符合實際情況的滑坡模型進行模擬分析。

顆粒流理論及其數值方法，克服了傳統連續介質力學模型的宏觀連續性假設，從細觀層面上分析巖土工程特性，并通過細觀力學參數研究分析宏觀力學行為，尤其適用于土石混合體的堆積體滑坡穩定性分析問題。

［1］ SKEMPTON AW.Long-term Stability of Clay Slopes［J］.Geotechnique，1964，14（2）：77－102.

［2］ 芮勇勤，賀春寧，王惠勇，等.層狀邊坡漸進破裂與失穩

過程數值模擬探討［J］.長沙交通學院學報，2002，18（3）：9－12.（RUIYong-qin，HE Chun-ning，WANG Huiyong，et al.Numerical Simulation Approach to Progressive Fracture and Instability of Stratiform Slope［J］.Journal of Changsha Communications University，2002，18（3）：9－12.（in Chinese））

［3］ URCIUOLIG，PICARELLI L，LEROUEIL S.Local Soil Failure Before General Slope Failure［J］.Geotechnicaland Geological Engineering，2007，25（1）：103－122.

［4］ 王志偉，王庚蓀.裂隙性黏土邊坡漸進性破壞的FLAC模擬［J］.巖土力學，2005，26（10）：1637－1640.（WANG Zhi-wei，WANG Geng-sun.FLAC Simulation for Progressive Failure of Fissured Clay Slope［J］.Rock and SoilMechanics，2005，26（10）：1637－1640.（in Chinese））

［5］ 陳亞軍，王家臣，常來山，等.節理巖體邊坡漸進破壞的試驗研究明［J］.金屬礦山，2005，35（8）：11－13.（CHEN Ya-jun，WANG Jia-chen，CHANG Lai-shan，et al.Experimental Research on Progressive Failure of Jointed Rock Slope［J］.Metal Mine，2005，35（8）：11－13.（in Chinese））

［6］ 殷坤龍，姜清輝，汪 洋.新灘滑坡運動全過程的非連續變形分析與仿真模擬［J］.巖石力學與工程學報，2002，21（7）：960－962.（YIN Kun-long，JIANG Qinghui，WANG Yang.Numerical Simulation on the Movement Process of Xintan Landslide by DDA Method［J］.Chinese Journal of Rock Mechanics and Engineering，2002，21（7）：960－962.（in Chinese））

［7］ CUNDALL P A，STRACK O D L.A Discrete Numerical Model for Granular Assemblies［J］.Geotechnique，1979（29）：47－65.

［8］ Itasca Consulting Group.PFC2DUser’s Manual（Version 3.1）［M］.Minneapolis，Minnesota：Itasca Consulting Group，2004.

［9］ Itasca Consulting Group.PFC2DTheory and Background［M］.Minneapolis，Minnesota：Itasca Consulting Group，2004.

［10］張曉平，吳順川，張志增，等.含軟弱夾層土樣變形破壞過程細觀數值模擬及分析［J］.巖土力學，2008，29（5）：1200－1204.（ZHANG Xiao-ping，WU Shun-chuan，ZHANG Zhi-zeng，et al.Numerical Simulation and Analysis of Failure Process of Soilwith Weak Intercalated Layer［J］.Rock and Soil Mechanics，2008，29（5）：1200－1204.（in Chinese））

［11］吳順川，張曉平，劉 洋.基于顆粒元模擬的含軟弱夾層類土質邊坡變形破壞過程分析［J］.巖土力學，2008，29（11）：2900－2904.（WU Shun-chuan，ZHANGXiao-ping，LIU Yang.Analysis of Failure Process of Similar Soil Slope withWeak Intercalated Layer Based on Particle Flow Simulation［J］.Chinese Journal of Rock Mechanics and Engineering，2008，29（11）：2900－2904.（in Chinese））

［12］河南省地質勘察院.河南省嵩縣白河鄉白河街村南山滑坡防治工程勘察報告［R］.駐馬店市：河南省地質勘察院，2011.（Henan Provincial Geological Survey Institute.Report of Survey on Baihe Landslide in Songxian County of Henan Province［R］.Zhumadian：Henan Provincial Geological Survey Institute，2011.（in Chinese） ）

（編輯：王 慰）

PFC Simulation of Progressive Failure Process of Landslide

WANG Yu1，2，LIXiao1，WANG Sheng-xing1，2，HOUWen-shi1，2
（1.Key Laboratory of Mechanics in Engineering Geology，Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing 100029，China；2.University of Chinese Academy of Sciences，Beijing 100049，China）

The Baihe landslide is a large ancientbedrock landslide induced by heavy rainfall，with the upper partof parallel layers and the lower part of slightly incised layers.According to field investigation，the causes and formation mechanism of the landslide are summarized.The particle flow code（PFC2D）program is employed and the parallel-bond model is introduced to determine the relation betweenmicro-parameters andmacro-mechanical properties of soils in the sliding zone through biaxial compression numerical tests.In subsequence，the landslide geological model is established and the progressive process of landslide is simulated numerically by the particle-flow discrete elementmethod.The displacement，porosity and strain at some key parts of the landslide aremonitored.The simulation results show that crack occurs first in the soils in the upper part along the bedrock plane，which leads to shear dilatancy of slip mass in the middle part，thus pushing the ancient slip mass in the lower part，and finally causes instability of the whole accumulation.The position and shape of sliding plane don’t need to be assumed，and the particles adapt their positions according to contact forces，and finally shear failure happens at the plane which has the smallest shear strength.The simulation could provide theoretical basis for decision-making of the project.

landslide；particle flow code（PFC）；numerical simulation；progressive failure；parallel-bond model

P642

A

1001－5485（2012）12－0046－07

10.3969／j.issn.1001－5485.2012.12.010 2012，29（12）：46－52

2011－10－17；

2011－11－18

國家自然科學基金資助項目（41027001）；國家重點基礎研究發展計劃（973）項目（2009CB724605）

王 宇（1985－），男，河北滄州人，博士研究生，主要從事地質災害與巖石力學等方面的學習研究，（電話）18046522933（電子信箱）good541571889＠126.com。