基于颗粒流方法的滑坡灾变过程及冲击致灾预测研究

2020-07-16陈鑫，赵洲，杨鹏

水力发电 2020年4期

陈鑫，赵洲，杨鹏

(1.西安科技大学地质与环境学院，陕西西安 710054；2.西安地质矿产勘查开发院有限公司，陕西西安 710100；3.中国有色金属工业西安勘察设计研究院有限公司，陕西西安 710054)

0 引言

滑坡是指由于受到地层岩性、水、地震及人类工程活动等因素影响，坡体沿贯通的剪切破坏面或带，以一定加速度下滑移动的不良地质现象[1]。随着我国经济的发展及各项基础设施建设向城市临近的山区不断拓展，滑坡问题越显突出，对人类生命和财产安全造成极大危害。据统计，全国平均每年地质灾害直接经济损失达200多亿元，其中滑坡造成的经济损失达2/3。因此，滑坡研究和防治工作一直以来是国家关注的重点。近年来，关于滑坡破坏运动过程及冲击致灾研究受到许多国内外专家学者的普遍关注[2- 6]。黄润秋等[7]采用有限元数值模拟方法，对膨胀土滑坡的特点和发生机理分别从强度分层、裂隙影响和蠕变效应进行的研究得出，膨胀土路堑边坡滑坡受上述3种因素共同作用，具体表现为浅层、小规模、渐进性和滑坡后缘通常出现在坡腰部位。亓星等[8]通过对甘肃黑方台滑坡野外调查和位移监测，探讨了滑坡变形破坏过程，并对滑坡的发生机制进行了详细解释。吴越等[9]利用光传感器和计时器组合而成的摩擦系数测量系统，根据相似原理建立滑坡滑体下滑冲击模型，探讨了滑体下滑和冲击致灾全过程中的能量耗散规律，为准确计算滑体冲击能奠定了基础。薛强等[10]采用三维数字高程，提出了滑坡强度定量评估技术方法，并对黑方台南缘区32处滑坡强度进行了评估，为滑坡风险评估提供了新思路。党杰等[11]根据物理学中内能转化原理，以地质灾害发生时释放的能量作为强度评估的依据，建立了滑坡强度与能量关系，提出了单体滑坡强度评估的指标和方法。以上滑坡研究方法大多处于定性或半定量研究，而定量研究方法仍然处于缓慢发展阶段。

“5·12”地震后，受地震影响，楼子沟滑坡土体产生松动变形，地面隆起，房屋裂缝变形加剧，少数房屋即将倒塌，滑坡后部出现弧形裂缝，醉汉林密布，严重影响21户、90人、100余间房及百余亩良田的安全，可造成经济损失达400万元以上。为此，本文采用PFC2D数值模拟方法，以楼子沟滑坡为研究对象，仿真动态模拟滑坡发生全过程，对滑坡冲击速度、冲击距离、冲击力大小等指标进行定量分析，同时对滑坡冲击致灾程度进行定量化评价。

1 滑坡概况

楼子沟滑坡位于陕西省汉中市勉县老道寺镇楼子沟村一组，属北部山区南缘丘陵区。滑坡平面呈一定弧长的长条形，主滑方向90°，滑坡体长约120 m，前缘宽度约360 m，后缘宽度约300 m，平均厚度约5 m，体积约2.16×105m3，属中型膨胀土滑坡。滑坡后缘以陡坡为界，陡坎高度约5～8 m；滑坡左、右两侧以小冲沟为界；滑坡前缘为村民水田，在边界外侧未发现变形迹象，为原始地貌；滑坡中后部地表可见明显下挫，边界清晰；滑坡后缘有泉水出露，且有一灌溉渠，渠中常年有水，水量较小。滑坡平面见图1。

根据收集已有资料和勘查结果，地层自上而下为第四纪堆积粉质粘土、下更新统冲洪积粉质粘土、第三系泥岩。滑体物质主要为滑坡堆积粉质粘土，黄褐色，土体结构松散，节理裂隙十分发育，土质呈斑块状，掰开后形如石榴粒，可见明显的灰褐色膨胀土矿物，网格状裂缝十分发育，最大宽度3 mm，灰色粘土薄膜附于裂隙表面，表面擦痕明显，裂隙面光滑，手捻有滑感。滑带土由灰褐色粉质粘土组成，呈可塑状，滑腻感很强，滑面光滑，可见擦痕。

图1 滑坡平面

2 滑坡模型构建

2.1 滑坡岩土体宏细观力学参数标定

颗粒流(PFC)方法采用中心差分法，把虎克定律和牛顿第二定律遍历整个颗粒集合体，模拟颗粒间的运动及相互作用，运用软件特定的粘结接触模型将其集合起来模拟岩土体的完整性，不受变形量的限制，可以方便地处理非连续介质力学问题，有效地模拟介质开裂、分离等现象[12]。在实际工程中，通过室内物理试验得到的宏观力学参数(如变形模量、内摩擦角、粘聚力等)可以直接在工程实际中运用计算，而在颗粒流数值模拟方法中，岩土体材料的宏观力学性质是通过颗粒总体的运动轨迹来反映的，而其运动受颗粒细观力学参数(如颗粒接触刚度、颗粒粘结强度、颗粒粘结刚度等)的影响。因此，细观力学参数的确定是颗粒流数值模拟分析的关键。鉴于关于颗粒细观力学参数和材料宏观力学参数之间的定量关系还没有确定的公式[13]，故PDC2D中需要采用双轴压缩试验模拟室内物理试验，以此得到细观力学参数与宏观力学参数之间的对应关系[14]。

本文采用的双轴压缩试验模型见图2。模型尺寸为10 m×5 m，颗粒粒径为0.1～0.15 m，生成颗粒总数为885个。滑坡岩土体颗粒细观力学参数见表1。在表1的细观力学参数下，施加围压分别为50、100 kPa和150 kPa，得到应力-应变关系曲线(见图3)，再根据M-C理论绘制M-C包络线(见图4)，最后通过计算得到双轴压缩试验中材料的宏观力学参数。综上，通过数值模拟得到的岩土体宏观力学参数粘聚力c=13.3 kPa，内摩擦角φ=12°，与室内物理试验中c=11.5～14.7 kPa，φ=11.0°～15.3°对比可知，数值模拟所得细观力学参数符合要求。

表1 滑坡岩土体颗粒细观力学参数

图2 双轴试验模型(单位：m)

图3 应力-应变关系

图4 M-C包络线

2.2 滑坡PFC2D模型构建

PFC中提供了10种内嵌的接触模型，本次数值模拟采用平行粘结模型，因为其既可以传递颗粒间的作用力，也可以传递颗粒间的作用力矩，从而更好的模拟颗粒间粘结破坏。本文采用PFC中ball-ball方法建立滑坡模型[15]，具体步骤如下：①根据滑坡主剖面图，按1∶1在PFC中建立矩形区域并生成颗粒；②通过命令导入AutoCAD中建立的滑坡剖面图并生成滑坡边界墙体使颗粒分组(见图5)；③利用PFC中ball delete range geometry命令删除削坡部位颗粒；④根据表1中颗粒细观力学参数对颗粒进行赋值，同时对颗粒施加重力作用使滑坡在重力作用下达到初始应力平衡状态。

图5 未削坡前滑坡模型

根据上述步骤得到楼子沟滑坡PFC2D模型，见图6。模型长165 m，高18.5 m，颗粒总数28 270个。其中，滑体颗粒总数7 058个，滑床颗粒总数21 212个。

图6 滑坡PFC2D模型

3 滑坡灾变过程及冲击致灾预测

3.1 滑坡灾变过程分析

为了更好地分析楼子沟滑坡变形特征及运动过程，根据滑坡模型尺寸及测量圆设置原则[16]，共建立4个测量圆，以此来监测滑坡滑体部位应力变化情况。滑坡位移见图7。从图7可知，当运行至4 500步时，滑坡坡脚出现剪切破坏，滑坡后缘出现拉张破坏，这主要是由于降雨导致滑体自重不断增加，当其超过锁固段强度时，滑坡发生剪出破坏；当运行至180 000步时，由于颗粒间相互摩擦碰撞消耗能量，滑坡动能逐渐减少，最终停止运动堆积在坡脚处。通过测量，模拟中颗粒最大位移为58.8 m，滑坡实际位移为36.9 m。

图7 滑坡位移

滑坡开始发生时，坡脚处剪切破坏。应力变化见图8。从图8可知， 0～13 s内1号测量圆中应力持续增大，由于滑坡体坡度较缓，故x方向的应力增量大于y方向，且近似为2倍；20 s后滑坡运动停止，应力达到平衡。12 s时2号测量圆内x方向应力达到峰值110 kPa，约为y方向应力的2.5倍，此时y方向应力来源主要为上部滑体下滑堆积挤压产生；3号和4号测量圆内应力逐渐减小，最终在11 s和4 s时分别归0，这是由于滑坡整体发生破坏，3号和4号测量圆位置处滑体部位由滑床滑出，故应力为0。

图8 应力变化

滑坡速度及位移分析是研究滑坡破坏运动过程的重要途径。楼子沟滑坡全长约120 m，本文以30 m为1段监测部位，每个部位分上、中、下3排，每排设置监测颗粒5个，对滑坡速度及位移进行定量描述，以加深对楼子沟滑坡破坏运动过程的认识。滑坡监测结果见图9。从图9可知：

图9 滑坡监测结果

(1)各个监测点运动主要分为加速阶段、过渡阶段、减速阶段。滑坡开始破坏时发生整体同等加速运动，滑坡体0～30 m内颗粒最大速度为3.6 m/s，随着滑坡前端剪出破坏，30～120 m范围内滑体由于前部失去支撑而向下运动，颗粒速度逐渐增大，最大速度达6.6 m/s。滑坡体0～30 m监测点颗粒减速阶段的历时远大于其他范围内监测点，这是由于前部颗粒不仅受30～90 m范围内颗粒的推动获得能量，同时还传递90～120 m范围内颗粒的能量，而90～120 m范围内颗粒由于没有后续能量补给，无法维持远距离运动，故在减速阶段迅速停止运动。

(2)滑坡体0～30 m范围内监测点颗粒最大位移为34.8 m；30～120 m范围内颗粒位移呈现依次增大特点，最大颗粒位移为58.8 m。该运动特征符合滑坡下部剪出破坏，逐级牵引滑坡体上部向下滑动的复合式滑坡规律。

综上所述，楼子沟滑坡运动过程可分为0～4 s等速变形阶段、4～15 s加速变形阶段、15 s之后堆积静止阶段3个阶段。滑坡破坏初期主要以蠕滑变形为主，随着形变量的增加，滑坡锁固段破坏，滑坡下部剪出，逐级牵引上部滑体顺破裂面下滑，最终堆积于坡脚处，滑坡全过程符合典型牵引式滑坡的运动特征。

3.2 滑坡冲击致灾预测

滑坡冲击致灾研究是滑坡风险评价的关键环节[17]，冲击力的大小直接反映了建筑物等承载体损毁的严重程度。因此，本文通过在滑坡不同坡脚距位置设置刚性挡墙代替建筑物，建立滑坡冲击模型(见图10)，以此统计滑坡下滑时对挡墙的冲击力，为滑坡冲击致灾分析提供依据。

图10 滑坡冲击模型

滑坡冲击力见图11。从图11可知，滑坡冲击挡墙过程可分为3个阶段：①自由下滑阶段。滑坡启动但未接触挡墙。②冲击阶段。滑坡不断冲击挡墙，滑体颗粒间及颗粒于挡墙相互碰撞耗能使滑坡能量逐渐减小。③冲击静止阶段。此时滑坡运动基本停止，滑体堆积于坡脚和挡墙之间，冲击力来自滑体堆积后对挡墙的静止土压力，为恒定值。

图11 滑坡典型冲击力

滑坡在冲击阶段内的冲量E可以根据冲击曲线和时间轴的面积得出，再根据等效力在相同时间内冲量相等的原理求出等效冲击力F，具体关系如下

(1)

对滑坡破坏全过程模拟分析可知，滑坡最大位移为36.9 m，37 m处冲击力为0。故从坡脚开始，每隔2 m设置刚性挡墙，共建立滑坡冲击模型19个，以此监测统计滑坡下滑对不同坡脚距处挡墙的冲击力。利用Origin数据分析软件计算得到不同坡脚距处等效冲击力数据，最后线性拟合绘制出滑坡坡脚距-冲击力关系(见图12)。由图12可知，滑坡等效冲击力随坡脚距的增大而逐渐递减，具体关系为F=-15 188.608 96x+543 387.274 59(R2=0.969 99)。滑坡对建筑物等承灾体的损毁程度随坡脚距增大而逐渐减弱。

图12 滑坡坡脚距-冲击力关系

根据现场实际调查资料，滑坡区内房屋大多以砖混结构为主，由于砌体结构主要发生脆性变形，故建筑物等承灾体不满足设计要求时认为是失效状态，等效冲击力作用位置由滑坡冲击结束时冲击高度及墙体所受荷载分布形式共同决定，由此建立建筑物失效模型，参考吴越等[18]承灾体易损性研究可定义滑坡冲击致灾指标R，具体关系如下

(2)

式中，F′为建筑物所能承受的极限荷载值。

由GB 50003—2011《砌体结构设计规范》可知，受弯构件承载力及受剪承载力，应按下列公式计算

M≤ftmW

(3)

V≤fvbz

(4)

Z=I/S

(5)

式中，M为弯矩设计值；ftm为砌体弯曲抗拉强度设计值；W为截面抵抗矩；V为剪力设计值；fv为砌体结构抗剪强度设计值；b为截面宽度；Z为内力臂；I为截面惯性矩；S为截面面积矩。

本文取单面墙宽，分别计算墙体受弯矩控制的荷载值F弯和剪力控制的荷载值F剪并取其中最小值作为建筑物所能承受的极限荷载值F′，即F′=min(F弯，F剪)。设置墙体距离坡脚20 m，根据图11方法计算得到此时滑坡等效冲击力F=2.45×105N，根据勘察资料，墙体材料为MU10普通烧结砖和M7.5水泥混合砂浆，墙体采用“三七墙”方式砌筑，墙体厚度365 mm，由式(3)、(4)、(5)计算得到墙体所能承受极限荷载值F′=min(F弯，F剪)=F弯=2.04×103N，将F、F′代入式(2)计算得到滑坡冲击致灾指标R=120.1。根据致灾程度分析可知，滑坡发生时墙体无法承受冲击损毁，会造成严重后果，对该滑坡应进行加固治理。