滑坡渐进破坏运动过程的颗粒流仿真模拟

王 宇，李 晓，王声星，侯文诗

滑坡渐进破坏运动过程的颗粒流仿真模拟

王 宇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。