地震诱发滑坡形成机理与动态模型研究

2015-01-03苏燕兰斯梅丁建辉

福州大学学报（自然科学版） 2015年3期

苏燕，兰斯梅，丁建辉

(福州大学土木工程学院，福建福州 350116)

据不完全统计，近几十年来地震滑坡造成了数万人丧生和几十亿美元的损失．强烈地震诱发的滑坡灾害，特别是在山岳地区，其危害比地震直接造成的损失还大［1］．国内外对地震型滑坡的研究多为统计意义上的工作．Keefer［2］着力研究历史地震事件，得出震级、地震烈度、震中距等地震参数与滑坡的关系;毛彦龙等［3］认为，地震时坡体波动振荡在斜坡岩土体变形破坏过程中产生三种效应:累进破坏效应、启动效应和启程加速效应;郑颖人等［4］基于汶川地震边坡调研，运用FLAC动力强度折减法对地震边坡破坏机制进行数值模拟，得出失稳是边坡潜在破裂区上部受拉破坏与下部剪切破坏共同造成的．这些研究成果有助于制定紧急救灾方案及减小震后生命和财产的损失，但不能直接用于提高大型工程结构的抗震能力和改进设计水平．从总体上看，对于地震诱发滑坡的研究，国内外基本上集中于面上调查层面．

福建省靠近欧亚板块和菲律宾海板块，地壳活动性较大，是中国东南部地震较多的省份之一，也是华南唯一发生过8级特大地震的地区．因此必须加强地震诱发滑坡的预测和预防．基于上述背景，根据福建某山区公路滑坡情况，利用基于散体介质理论的离散元程序PFC2D，对地震荷载作用下的边坡进行动态模拟，探求地震型滑坡的形成机理．

1 数值模拟方法

离散元数值模拟分析中为了建立波在连续介质与离散介质中的传播等效关系，将离散颗粒连成串看成一根截面半径为2R(R为颗粒的半径)的弹性棒，如图1所示，通过分析波在弹性棒上的传播情况来模拟分析颗粒的受力情况．当应力波作用导致岩体发生破坏时，颗粒离散元方法既可体现出大变形亦可模拟其破坏过程［5］．

图1 波在离散介质与连续弹性介质中的等效传播Fig．1 Equivalent transmission in the discrete medium and the continuous medium

1．1 输入地震荷载

根据福建省地震设防要求，选用适用于Ⅲ类场地，按7度进行设防的典型波Elcentro波(历时22．12 s)东西方向进行分析，将地震波的峰值按同比例调整到0．981 m·s－2，相当于0．1g(如图2)．由于离散元程序PFC2D中墙体只能加速度而不能加力，将加速度时程曲线进行积分得到速度时程曲线，如图3所示，并采用颗粒边界条件，将地震波转化成力加载到颗粒边界上．输入的地震波传播方向自下而上，垂直向荷载振幅取为水平方向的2/3．

图2 加速度时程曲线(峰值0．1g)Fig．2 Acceleration versus time

图3 速度时程曲线Fig．3 Velocity versus time

1．2 确定参数

颗粒流是以圆球或圆盘来描述实际现象，颗粒间的相互作用服从经典力学理论，程序中采用接触本构模型，通过颗粒间的接触变形与应力的关系来模拟散体材料特有的作用形式．对散体介质而言，颗粒间的细观力学特性决定了其宏观力学性质．室内试验得到的大多是宏观参数，在PFC2D数值模拟中，通常借助双轴压缩试验，不断调整细观参数使其宏观力学特性与室内试验或现场试验现象一致［6］．本文基于对粘土平面应变试样的颗粒流模拟研究成果［7］，为了更好地逼近土体微观上的各向异性和不均匀性，在生成PFC2D试样时由不同半径的颗粒单元组成，考虑相似比，经大量试算，控制参数如表1所示．

表1 数值模型计算参数Tab．1 Parameters of numerical model

1．3 建立模型

根据福建省山区滑坡情况，选择几何相似比为20，边坡简化成如图4所示模型(模型高100 cm，长140 cm，顶部宽50 cm，坡度为53°)，对地震条件下边坡发生滑动破坏进行数值模拟．考虑到离散元软件的计算量较大，现有的计算机硬件水平难以完成颗粒的等比例模拟，在满足精度要求的前提下放大颗粒尺寸，并利用PFC2D内置函数按均匀分布粒径生成试样．在坡体内部不同位置设置测量圈(如图4)，追踪测量圈内的土体平均应力、孔隙率、接触数等参数的变化规律．

图4 滑坡数值模型与测量圈布置示意图Fig．4 Landslides numerical model and measurement circles

2 结果分析

2．1 位移场分析

使用Fish语言实现地震作用下颗粒位移云图的实时绘制，采用10种颜色来代表颗粒位移的大小(见图5)，通过地震作用下边坡不同时刻的颗粒位移云图(见图6)可判断边坡的破坏情况．边坡在地震作用下，临空面位移较大，首先出现松动，见图6(a)灰色部分;接着沿临空面向里发展，见图6(b)，灰色部分向边坡内部扩大;有顶角的坡先是顶角坡面出现裂纹(见图6(c)、(d))，导致边坡从顶角开始滑动;边坡开始是局部小范围松动滑坡，但最后呈现出大体积下滑的现象，滑动面为圆弧形，如图6(e)～(g)所示．

图5 不同位移的颗粒颜色示意图Fig．5 Particle displacements in different colors

图6 地震作用下边坡不同时刻的颗粒位移云图Fig．6 Particle displacement of slope during seism

2．2 应力场分析

2．2．1 粘结力分析

地震作用下，边坡土体内局部出现粘结力丧失现象，且有规律地出现裂缝，形成圆弧面，该裂缝随地震时程的增加而增大并贯通，最后坡体沿着圆弧面整体下滑，产生地震滑坡现象(如图7所示)．

图7 地震作用下边坡粘结力示意图Fig．7 Variation of slope cohesive force during seism

2．2．2 切应力及竖向应力分析

从数值模型位移和粘结力随时间变化的过程可以看出，圆弧滑裂面上参数变化较大，故选取A4组测量圈(见图4)进一步分析．测量圈的切应力见图8，测量圈内切应力突变后减小，减小越快说明测量圈对应处的土体颗粒滑动速度越快，测量圈15甚至出现负值，说明该处土体颗粒在地震作用下发生滑动．测量圈的竖向应力见图9，由竖向应力的变化规律知，竖向应力发生稳定上升或稳定下降是滑坡发生的地方，竖向应力值变大说明该测量圈对应处或附近坡体发生下滑，竖向应力值变小说明该测量圈处或附近是下滑坡体堆积的地方．

图8 A4组测量圈的切应力Fig．8 Measure shear stress of group A4

图9 A4组测量圈的竖向应力Fig．9 Measure vertical stress of group A4

2．3 细观结构分析

孔隙率变化可以反映土体的运动趋势和破坏过程．从地震作用下各组测量圈内的孔隙率时程曲线(图10)可以看出，在地震作用很短时间内，孔隙率存在一个短时间突变(上升后下降)过程，这是由于在地震荷载作用前颗粒在自重作用下已经达到初始平衡，地震荷载作用后土体发生扰动，导致孔隙率增大．对于孔隙率发生显著上升的，说明该测量圈对应位置在地震作用下发生滑坡，且形成一个临空面;若发生明显下降，说明该处是滑坡后坡体堆积的地方．

图10 A4组测量圈的孔隙率Fig．10 Measure porosity of group A4

图11 A4组测量圈的滑动比Fig．11 Measure sliding fraction of group A4

滑动比是指在测量圈内发生滑动的接触点数与总点数的比值．研究滑动比的变化规律可以进一步判断滑坡过程中颗粒与颗粒发生位移和旋转的剧烈程度(图11)．从图11可知，滑动比随地震时程的增加而发生激烈的跌宕起伏，5号测量圈从数量级上看是百分之零点几，可以认为几乎没有变化．13号测量圈的滑动比由百分之零点几增长到百分之四点多，说明地震过程中13号测量圈内颗粒接触发生较大滑移．通过上面分析可知:一般情况下，测量圈内滑动比增大速率越大，该处颗粒发生滑动量越大，若处于一直增大的过程，说明该测量圈所在处可能已发生滑坡．

3 影响因素分析

3．1 粘结刚度

图12为同一边坡在典型Elcentro波作用下，采用不同刚度，边坡滑坡情况的示意图．由图可知，随着粘结刚度的增加，滑裂面越来越陡．由图12(b)～(d)看出，刚度大到一定范围，其在滑坡体出现裂缝增多且明显，呈现出脆性状态，且土体内部颗粒的位移也在减小．当刚度增到一定值后坡面不再破坏，临空面的颗粒还是有一定的位移量，但只呈现出要破坏的趋势．

图12 地震作用下不同粘结刚度坡体位移云图Fig．12 Landslides displacement of different bond stiffness during seism

3．2 坡度

使边坡各参数相同，改变坡面坡度，在典型的Elcentro波作用下滑坡情况如图13所示．从图13(a)可以看出，坡面坡度为37°时，地震作用下坡脚颗粒位移较大，坡脚发生局部滑动概率较大;图13(b)表明，对于45°边坡，其颗粒发生位移较大的地方从坡脚开始沿坡面呈现出圆弧状，但没有贯通到坡顶;而图13(c)有完整的滑坡圆弧面且出现较明显的滑坡现象．由此可以得出结论:对于材料和其他性质相同的边坡，在相同的地震和地震时程作用下，坡度越小，地震作用下颗粒发生大位移(模型中灰褐色的地方)的颗粒范围越小，颗粒位移较大处发生滑坡现象的概率相对其他颗粒较大，即坡度越小的坡越稳定．