基于颗粒流方法的滑坡破坏机理与强度分析

2018-08-10魏江波

西安科技大学学报 2018年4期

赵洲，魏江波

(西安科技大学地质与环境学院，陕西西安 710054)

0 引言

基于数值分析的滑坡破坏过程模拟、机理分析与强度评价是滑坡灾害研究的热点之一。目前，国内外对滑坡破坏机理等方面的数值模拟研究多采用有限单元法(FEM)和有限差分法(FDM)，这2种方法均假设研究体为连续体，忽略了岩土体颗粒间的复杂相互作用及实际岩土体的高度非线性行为，对大变形大位移求解尤其土质滑坡分析也不够理想，对整个滑坡的稳定情况评价力度不足。颗粒流方法(particle flow code，PFC)是一种特殊的离散单元法[1-2]，通过刚性的圆形或者球体来模拟岩土体颗粒，以力-位移定律和牛顿第二定律为基本理论模拟颗粒的运动和相互作用，运用特定的粘结模型将其集合于一体来模拟岩土体的整体性，用滑动模型模拟颗粒之间的剪切运动及岩土体宏观大变形，通过细观颗粒的运动特征来反映宏观岩土体的变形破坏特征。同时，颗粒流方法最大优点在于不用假设滑动面位置以及可直观的观察岩土体的动态真实破坏过程[3]。

近年来运用颗粒流方法对滑坡以及边坡的研究也越来越多。周健等、廖彪应用颗粒流方法分析了土体细观参数对土质边坡破坏形式的影响，认为随着土体颗粒黏性的增大，其破坏类型从塑性破坏向脆性破坏发生过渡，并得出了颗粒流方法与传统方法对均质土坡稳定性计算结果比较接近的结论，为边坡稳定性分析开辟了一条新的途径[3-5]。张龙以鸡尾山岩质滑坡为例，对滑坡的运动过程、破坏方式及堆积规律进行了颗粒流分析，认为摩擦系数及土体强度对滑体物质在堆积区的分布具有重要影响，但对滑体的最大位移影响不大[6]。马建全对黑方台焦家崖头灌溉影响下的黄土滑坡进行了基于PFC2D的颗粒流模拟，展现了滑坡的整个运动过程，并分析了滑坡各阶段的破坏特征[7]。廖静薇基于试验分析与颗粒流理论，建立了粉质粘土宏观力学参数与细观力学参数的对应关系，提出了细观力学下粉质粘土边坡的失稳判据和边坡安全系数计算方法，并结合实例应用证明了颗粒流方法在土坡稳定性分析中的有效性[8]。N.Li等利用PFC进行了荷载作用下土质边坡的变形破坏机制研究，并与室内物理模拟试验进行了对比，取得的模拟结果十分相近[9]。Valentino等通过颗粒水槽试验模拟了颗粒流的运动过程，并采用离散元软件PFC2D进行数值模拟水槽试验，同时获取颗粒流冲出距离和冲击力频谱特性[10]。唐红梅等借助PFC2D软件对三峡库区龚家方2号斜坡破坏过程进行了数值模拟，得出模拟结果与实际破坏情况相一致的结论[11]。吴越等基于功能关系，通过物理实验与数值模拟结合对比，从滑坡力学特征的能量角度研究了滑坡运动过程中的能量耗散以及滑体对承载体的冲击作用，建立了承灾体损毁程度与滑坡冲击作用的定量函数关系[12-15]。

由此可见，目前采用颗粒流方法对滑坡的研究多集中在滑坡颗粒流模型的构建、岩土体细观参数对滑坡破坏方式的影响、滑坡运动过程的模拟、滑坡稳定性计算的对比分析等方面，而对滑坡破坏过程中的力学机理分析、滑坡速度、动能及冲击力等反应滑坡强度变化特征指标的分析和研究相对较少。

利用颗粒流软件PFC2D，以府谷县新府山滑坡为例，建立了滑坡颗粒流模型，模拟了滑坡从变形失稳到运动停止的动态破坏全过程，分析了滑坡变形破坏特征及其力学机理，研究了滑坡强度的变化特征，给出了以冲击力为表征的滑坡冲击破坏强度定量分析方法。

1 滑坡概况及岩土体参数确定

1.1 滑坡概况

新府山滑坡位于陕西省府谷县府谷镇新府山安居小区北侧斜坡(图1)。该滑坡体前后缘相对高差25 m，宽约148 m，长约40 m，厚约3～3.5 m，体积约1.1×104m3，坡体中部较缓，坡顶、坡脚部较陡，滑坡坡面总体近直线型，倾向东南。滑坡周界分明，左右两侧以冲沟为界。滑坡体为第四系全新统褐色风积粉土，下伏上古生界二叠系棕红色泥岩，产状200°∠35°。由于坡体前缘开挖使得坡体下部变陡，加上降雨下渗作用使得坡体强度降低，降雨汇集于泥岩顶面，软化泥岩结构，降低抗剪强度，从而诱发坡体沿粉土-泥岩接触面下滑，属典型的土—岩接触性滑坡(图2)。

图1 新府山滑坡平面位置及发育特征示意图Fig.1 Location and characteristics of Xinfushan landslide

图2 新府山滑坡地质剖面Fig.2 Geological section of Xinfushan landslide

1.2 滑坡岩土体宏细观参数确定

岩土体力学参数是滑坡稳定性评价的关键和难点[16]。文中根据土样试验及工程经验分析，确定的新府山滑坡相关岩土体宏观力学参数见表1.

PFC赋予岩土体的是其颗粒的细观参数，包括颗粒密度、粒径、法向刚度、切向刚度、摩擦系数等，并通过颗粒的细观参数来反映岩土体的宏观特征。

表1 新府山滑坡岩土体力学参数Table 1 Rock and soil mechanical parameters of Xinfushan landslide

为了获得宏观参数相对应的细观参数，利用PFC2D程序中的双轴试验模型方法构建了6 m×12 m的矩形模型(图3)，上下墙(wall)施加压力，左右墙模拟约束。在构建好的矩形墙内，用半径扩张法结合generate命令生成一定数量一定空隙率大小的颗粒[17-20]。为了更好的模拟岩土体的真实特征，颗粒(ball)半径越小越好，但受计算机的容量与运行速度的限制以及所需精度要求，可适当扩大颗粒半径。最后根据土体颗粒特性，选用颗粒的接触粘结模型，并用property命令生成接触粘结模型。

在双轴试验模型建立好之后，赋予颗粒一组细观参数及不同围压，进行模拟并获得对应偏应力-应变曲线，进而获得对应的摩尔应力圆。进行大量的试算模拟，通过每一组细观参数与不同的围压所获得的摩尔应力圆绘制出对应的强度包络线并反算出宏观参数与表1参数进行对比，获取与表1相匹配的那组细观参数作为实际岩土体的模拟参数。借鉴文献12中的宏细观参数的定量关系式[21]，并通过多次试算获得了新府山滑坡岩土体颗粒的细观参数(表2)。

图3 双轴试验模型Fig.3 Biaxial test model

表2 新府山滑坡岩土体颗粒细观参数Table 2 Rock and soil microscopic parameters of Xinfushan landslide

以粉土为例，赋予表2的细观参数，对粉土进行围压分别为50，100和200 kPa的双轴试验模拟，并得知50 kPa围压下峰值强度为190，100 kPa围压下峰值强度为310，200 kPa围压下峰值强度为570 kPa，则可绘制Mohr应力圆与Mohr-Coulumn破坏包络线，计算得粘聚力约为21.3 kPa，内摩擦角约为26.8°(图4)。对比表1可知，表2中的细观参数可用于后文滑坡模型模拟。

图4 Mohr应力圆与Mohr-Coulumn包络线(粉土)Fig.4 Mohr stress circle and Mohr-Coulumn line(silt)

2 滑坡颗粒流模型构建

为了更真实的模拟并反映滑坡的特性，解决传统半径扩张法存在的岩土体层间无接触问题与放大系数难确定问题，文中结合“分层下落法”的思路构建滑坡模型，该方法不需要建立坡面约束墙和层间分界墙，不用消除浮球，且存在层间接触[22]，与实际情况更加吻合。根据图1尺寸与表2参数构建滑坡模型，过程如图5所示。

在建模过程中，以单元平均不平衡力达单元平均接触力的1%，且最大不平衡力达最大接触力的1.0×10-5时作为平衡条件，即循环停止，且平衡后进行速度与位移清零。为便于区分上下层岩土体，将其进行分组(group)，并设置不同颜色。保留斜坡边界无刚度墙便于观察颗粒运动过程中的情况，其只作为参考作用，并无阻碍作用。

初始模型构建完成后，删除右侧竖直墙(wall)，并在坡脚处建一水平刚性墙，设置其摩擦系数(friction)为1.0，用以模拟坡脚处水平地面。同时为了研究分析滑坡体应力变化，在滑动面上中下不同位置设置测量圆，以监测(history)记录对应位置在坡体滑动过程中的应力变化状况。最终模型如图6所示。

图5 分层下落法建模过程Fig.5 Modeling process of stratified rain fall

3 滑坡破坏过程模拟与机理分析

3.1 基于颗粒流强度折减法的滑坡稳定性分析

颗粒流强度折减法就是将颗粒的接触粘结强度和摩擦系数同时折减一定的比例，使得滑坡刚好处于临界位移状态，以临界位移作为滑坡破坏判断依据，此时，将折减前强度参数(n_b，s_b和fric)与折减后的临界强度参数(n_bcr，s_bcr和friccr)的比值定义为滑坡稳定性系数。其公式如下

式中n_b为颗粒的法向接触粘结强度；s_b为颗粒的切向接触粘结强度；fric为颗粒间的摩擦系数；n_bcr为临界状态时的颗粒法向接触粘结强度；s_bcr为临界状态时的颗粒切向接触粘结强度；friccr为临界状态时颗粒间摩擦系数。

图6 新府山滑坡颗粒流模型Fig.6 PFC2D model of Xinfushan landslide

根据公式(1)进行滑坡强度折减模拟。起初设定折减系数为1.0，逐渐增大折减系数，并进行模拟观察颗粒位移场的变化，当折减系数增大至1.115时滑坡刚好处于临界位移状态，此时的折减系数即为滑坡的稳定性系数，与采用SLOP/W模块求得的滑坡稳定性系数为1.116的结果十分接近[23]。可知，该滑坡在自然工况下处于基本稳定状态，但在人类工程活动、强降雨等诱发因素的影响下，不排除其复活的可能性，一旦滑坡强度达到临界状态时，该滑坡将继续失稳破坏。因此，有必要对滑坡的破坏过程及其机理做进一步分析。

3.2 滑坡破坏过程模拟与破坏机理分析

3.2.1 从变形特征分析滑坡破坏机理

选取上述强度折减法所求的临界位移状态时的细观参数，即将岩土体细观强度参数折减1.115倍，采用PFC2D程序对滑坡破坏过程进行仿真模拟。

根据颗粒流模拟结果分析，滑坡在初始时段处于变形状态，应力重新分布，坡脚颗粒受上覆压力作用下变形不断积累，3.0×104步时坡脚挤压破坏(图7(a))。坡脚挤压破坏后，坡体失去支撑，变形加剧，颗粒根据自身的接触力自行调整位置，并且使得坡体中部抗剪强度最弱处产生剪切破坏。在上覆压力作用下，变形继续加剧，坡体中的剪切裂隙处产生应力集中，当裂隙尖端应力大于其抗剪强度时，裂隙尖端继续破坏，裂隙延伸，表现为坡体中剪切裂隙向上下附近的粘结最弱处延伸。变形过程中，破体向下滑移，坡顶失去支撑，5.0×104步时坡顶产生明显拉张裂缝(图7(b))。当上部拉张裂缝与中部剪切裂隙相贯通，并与坡脚挤压破坏处相连接后，坡体将失稳破坏，产生整体破坏并滑移。当滑坡整体破坏后，上部台阶处坡顶颗粒向下滑落，达到一定时步后，其台阶之上的原滑坡后壁所形成的小边坡失去下部支撑，1.8×105步时发生次生滑坡(图7(c))。

图7 不同时步滑坡破坏形态Fig.7 Failure form of landslide at different time

在滑坡变形—失稳—运动过程中观察颗粒间的粘结(即黑色短线)变化情况(图8)，可看出，在滑坡再次失稳破坏过程中，滑体颗粒间的粘结基本破坏，坡体愈加破碎，并顺势下滑堆积于坡脚。对于以位移分析法为主的颗粒流法，在滑坡运动过程中，通过位移矢量图可以直观看出颗粒自动搜索的滑动面位置(图9)。

图8 t=150e4时步接触粘结图Fig.8 t=150e4 time step contact bond

3.2.2 从应力特征分析滑坡破坏机理

根据滑坡破坏过程中应力变化曲线的分析可知(图10)，坡顶由于产生拉张裂缝和垮塌，所以1号测量圆内y方向的应力明显减小，x方向略有减

图9 t=150e4时步位移矢量Fig.9 t=150e4 time step displacement vector

小，最后当滑坡破坏后，1号测量圆所在滑坡后壁处于稳定状态，故该位置应力最后均趋于稳定(图10(a))。对于坡中由于坡体持续性向下滑移，且破碎程度不断增大，使得2号测量圆所处位置上覆岩土体不断减少，因而其所受竖向压力减小，所以2号测量圆内y方向应力明显降低，由于上部坡体破坏对下部岩土体的挤推作用与持续下滑使得x方向应力起初略有增大后又稍微减小(图10(b))。对于坡低，由于坡体下滑挤推作用及堆积作用使得3号测量圆内x和y应力都显著增加(图10(c))。

图10 应力变化模拟曲线Fig.10 Stress change simulation curve

4 滑坡强度分析

滑坡强度是表达滑坡所储备的破坏能力及对承灾体破坏程度的一项重要指标，也有学者认为滑坡强度指滑坡具有的能量大小[24]，其研究对滑坡风险评估具有重要意义。滑坡强度大小与滑坡规模、斜坡坡度、坡高、滑坡岩土体类型等滑坡发育特征均具有一定联系，然而目前仍没有一种确定的滑坡强度定量分析与评价的方法。一般认为，滑坡体的位置越高、滑体体积越大、移动速度越快，则滑坡强度也就越高，对承灾体的破坏能力也就越大。

为了分析新府山滑坡在变形与破坏过程中强度变化特征，文中从速度、动能、冲击力等指标入手，间接分析新府山滑坡强度大小的变化特征。

4.1 滑坡速度与动能变化特征分析

由于颗粒流模型中难以监测出滑坡整体的运动速度，因此文中选取滑坡破坏后颗粒最大速度vmax进行分析(图11)，以及进行实时监测滑坡动能的变化(图12)。虽然颗粒最大速度不能代表滑坡的整体运动速度，但就特定滑坡而言，仍然可以反映出滑坡能量大小的变化特征。

从图11，图12对比可以看出，滑坡速度特征与动能特征基本一致。滑坡在1.0×105时步破坏时，虽然速度最大，但滑坡破坏体积较小，所以此时的动能并非最大。滑坡破坏后，由于滑铲作用使得滑坡破坏体积增大，直至1.5×105时步时，滑坡的破坏体积和速度都处于较大状态，此时动能达到最大。随后，由于滑坡在水平地面的缓慢运动使得动能逐渐降低。这一变化过程也符合滑坡破坏的突发性和堆积的缓慢性特征。

因此，从速度和动能指标来看，新府山滑坡的强度在其破坏的初始阶段为最高，然后逐渐变小。

图11 颗粒最大速度变化曲线Fig.11 Maximum velocity variation curve of particles

图12 滑坡运动过程动能变化曲线Fig.12 Landslide kinetic energy variation curve

4.2 滑坡冲击力变化特征分析

滑坡强度研究的主要目的之一是评价滑坡对建筑物等承灾体的冲击破坏能力大小，强度越高的滑坡，对承灾体的破坏能力也越大[25-26]。从致灾后果来看，无论滑坡动能、速度等指标如何变化，若不考虑承灾体自身抗灾能力特征，就滑坡自身而言，最终反映建筑物等承灾体损毁程度的变量是滑坡冲击力大小。因此，为了分析新府山滑坡强度即滑坡对承灾体的冲击破坏能力大小变化特征，文中选择滑坡冲击力这一指标来进行滑坡强度大小分析。

从图11，图12的分析可以看出，随着新府山滑坡运动的停歇，其速度、动能等呈衰减趋势，也就是说，在水平地面上距离滑坡坡脚不同位置(文中将其称为坡脚距)处的承灾体，其遭受滑坡的冲击力大小不同，破坏程度也不同。

通过在沿滑面倾向方向上不同坡脚距(用L表示)处构建竖直刚性wall，并分别监测滑坡破坏过程中坡脚刚性wall上的冲击力随时间的变化，以此来分析滑坡对不同坡脚距处建筑物的冲击作用，从而间接反映出滑坡强度大小及其变化特征。通过上文模拟分析和计算，新府山滑坡的最大滑移距离约为16 m，因此，在坡脚距从0 m开始每间隔1.0 m分别设置滑坡冲击力测试的刚性wall，共建立17个冲击力测试模型(图13)，分别实时监测刚性wall上的水平冲击力大小(图14)。

图13 刚性wall设置示意图Fig.13 Rigid wall setup illustration

图14 滑坡冲击力随时间的变化曲线Fig.14 Variation curves of impact force with time

由图14可知，当坡脚距一定时，起初滑动阶段，由于滑坡速度与破坏的滑体体积不断增大使得冲击力随滑坡的运动而增加，当滑体减速堆积并停止运动时，冲击力逐渐趋于稳定(此时的冲击力表现为滑体对刚性墙的静止土压力)。同时，当坡脚距不同时，冲击力随坡脚距的增大明显减小。然而由于滑坡体特征及运动特征等因素影响，导致滑坡冲击力峰值基本上接近于滑坡体失稳运动后的稳定状态下堆积体对挡墙的静止土压力。

对不同坡脚距上的最大冲击力进行统计分析可知(图15)，滑坡最大冲击力Fmax与坡脚距Li之间的关系可以表达为

Fmax=761.46e-0.17/Li，(0≤Li≤16)

式中Fmax为滑坡最大水平冲击力；Li为坡脚距。

图15 滑坡冲击力随坡脚距的变化曲线Fig.15 Variation curves of impact force with Li

因此，以冲击力表达的新府山滑坡强度大小随着滑坡在水平地面运动距离的增大而减小，并在运动停止前呈现出指数衰减的特征。从滑坡破坏及运动过程分析，当滑坡体在整体发生破坏并顺着滑面运移到斜坡坡脚处时，由于具有较大的体积和速度，因而在坡脚处具有较大的冲击力即破坏强度，但随着滑坡体与水平地面的碰撞、摩擦运动等，滑体内部能耗、与地面之间的摩擦能耗增加，导致滑体颗粒运动速度变缓，冲击力减小，滑坡强度也随之降低，直到滑坡停止运动强度为零。