基于PFC2D的柿子树坪滑坡变形特征模拟研究

2023-03-15惠航刘锋冯兵赵法锁

科学技术与工程 2023年3期

惠航，刘锋，冯兵，赵法锁

(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]。