地震作用下滑坡动力响应及稳定性模拟分析

肖 衍，袁宝远，程 岩

(河海大学地球科学与工程学院，南京 211100)

斜坡受各种自然和人为因素的作用会发生多种形式的变形和破坏，其中滑坡是普遍的破坏方式，滑坡会给工程建筑和生命安全造成难以避免的巨大损失。我国中西部以高山、丘陵为主，斜坡的稳定问题突出，特别是地震作用下，滑坡破坏频发，5·12汶川大地震，沿着约两百公里长的龙门山断裂带，造成了许多的崩塌滑坡，在此次灾难中，其中近3成人员伤亡跟崩塌滑坡有密切联系。

滑坡在地震作用下稳定性分析方法主要有极限平衡分析方法[1-3]和数值分析方法[4-5]。拟静力作用[6]下极限平衡分析法应用很多，拟静力作用下数值分析法也有不少应用。但考虑地震作用下的动力响应[7-8]，进一步分析滑坡的地震动力作用稳定性的实例还较少。

白龙江流域的下游段，两岸山体雄厚，高程为1 500～3 000 m，相对高出河水位600m以上绝大多数为基岩岸坡，整体稳定性好，影响库岸稳定的主要是松动、倾倒变形体和滑坡。郭家滑坡为一古滑坡体，滑坡的纵向长度870m，平均宽400m，面积0.34km2，滑坡体厚度大于57m，体积2.35×107m3。其下游3km为一在建大中型水电站，其稳定性会对下游水电站坝基稳定以及库岸稳定性带来直接的威胁，该区50a超越概率10%时的基岩水平峰值加速度0.195g。因此，在该区进行地震安全性分析评价意义重大。

1 滑坡体工程地质条件

滑坡体位于甘肃省文县口头坝乡，该区属陇南构造剥蚀及侵蚀的中高山区，地形总体上西高东低，并由东、西两侧中高山向白龙江河谷倾伏。白龙江属降雨补给型河流，库区位于中游河段，区内年降雨量400～500mm，多阵性暴雨，长历时降雨少，年降雨量多集中在7、8、9月份。

滑坡体前缘高程800m，大约高出河水位90～100m，后缘高程1 400m左右，地形上缓下陡，上部24～30°，中部34°，下部41°，从平面形态上来看，滑坡呈长舌状，长轴方向SE155°，与主滑方向基本一致。根据2017年8、9月份现场勘察与取样分析，晚更新世的风成黄土布满了郭家滑坡的后部边界地带，滑坡的中以及前半区域在坡面上能够划出两个区，上游区主要是凝灰岩和局部瓦解的板岩构成，下游区物质主要为板岩碎块，从剖面来看，具有清晰的分层性，由表及里共可以分为6层:滑坡体的风积黄土主要分布在后缘的缓坡和1 200m高程以上地带；变质凝灰岩碎石层主要见于滑坡的浅表部，分布高程低于1 200m，主要是块石，其所占比例可达80%，垂直向的厚度23～35m；保留原岩层序的变质凝灰岩块石层大部分见于中浅部，依然维持着原岩的层序结构，根据平硐揭露水平向厚度60～65m；保持层状结构的滑体,主要分布在滑坡体中深部，厚度40～45m；板岩强烈揉皱层，厚度25～30m；滑带土层，厚度7.5m左右，灰白色的岩粉、土黄色—深褐色泥质物是其主要成分，含有角砾，角砾粒径0.3～1cm，较密实，有潮湿感，手捻有砂感。

河流最高水位800m，低于滑坡前缘高程。

2 动力计算模型

选用FLAC3D数值分析软件进行地震动力作用下滑坡稳定性数值分析计算。计算的模型见图1，模型沿y向取一定的厚度进行建模。模型计算底部和左右边界是固定约束，上部表面是可以自由移动的。根据所取样品进行的室内试验结果，参考收集的白龙江下游苗家坝工程的试验资料，综合分析得出计算所取参数见表1。材料本构模型采用摩尔库伦模型。

表1 材料计算参数

图1 滑坡数值计算剖面模型Figure 1 Landslide numerical computation section model

动力作用计算选取加速度时程从模型底部输入达到加载目的，加速度时程根据人造波合成最大加速度为73cm/s2，持续40s，两个加速度之间相隔为0.02s，其时程曲线如图2所示。

时间/×10s图2 地震波时程曲线Figure 2 Seismic wave time-history curve

边界条件的设置有静态边界和自由场边界两种方式。静态边界添加了阻尼器来接受输入的波，自由场边界通过把一二维网格施加到所计算的物体侧围来添加这种边界。本文计算选择自由场边界。介质之间相互摩擦相挨的表面滑动都会有力学阻尼，FLAC3D所提供阻尼有瑞利阻尼、局部阻尼和滞后阻尼。本次计算选用瑞利阻尼来求解。

3 动力响应分析

加速度可以用来表征地震力，通过在坡面和坡体内布置一系列监测点，来观察滑坡体对地震波的响应规律，其中所布监测点高程分别为805、863、936、1 041、1 120、1 201、1 280、1 358m，滑坡体内某点所达到的最大加速度与加载的最大加速度相比所得加速度放大系数，图3和图4是记录的高程在1120 m坡表面和坡体内部特征点加速度随时间变化曲线。

图3 特征点(287,0,758)加速度时程曲线Figure 3 Characteristic point (287, 0, 758) acceleration time-history curves

图4 特征点(231,0,767)加速度时程曲线Figure 4 Characteristic point (231, 0, 767) acceleration time-history curves

由图3、图4可得坡面特征点最大值3.39m/s2，加速度与坡底相比放大了4.65倍，坡内特征点最大值6.82m/s2，与坡底相比放大了9.65倍。

定义坡内和坡面的加速度放大系数分别为α和β，加速度放大系数随高程变化如图5所示，可以看出加速度向上传播时候均有高程放大效应，不同高程处放大系数不同。

4 稳定性分析

利用FLAC3D软件计算滑坡的稳定性。 首先计算滑坡在重力和地下水作用下的稳定性。图6和图7是剖面在自重条件下最大主应力和最小主应力分布图，由图6、图7可知，在自重下坡体内应力均为负值，所以拉应力不存在，第一主应力最大20MPa，第三主应力最大5.1MPa，均出现在底部。

图5 加速度随高程动力响应曲线Figure 5 Dynamic response curve of accelerationin pace with elevation

图6 剖面最大主应力分布Figure 6 Section maximum principal stress distributions

图7 剖面最小主应力分布Figure 7 Section minimum principal stress distributions

图8是剖面在自重条件下所得剪应变增量图，剪应变增量最大值在0.012，根据FLAC3D自带求解安全系数的方法，求解的剖面在自重下安全系数为1.13。图9为求解安全系数下剪应变增量图， 从图中可以看到强度减小后剪应变增量相比之前数值突增，贯通的剪应变增量在0.02～0.01。

图8 剖面剪应变增量Figure 8 Section shear strain increment

在计算滑坡在重力和地下水作用下的稳定性的基础上，计算滑坡加上地震动力作用下的稳定性。图10是加40s地震波后的剪应变增量图，图11是地震工况下剖面塑性区分布,图12是折减系数1.05剪应变增量图，当系数折减1.05时，剪应变最大值已变为0.1，剪应变增量贯通区在0.03～0.01，可以认为滑坡稳定性折减系数在1～1.05。

为了地震动力数值分析能与极限平衡方法进行对比，本文利用拟静力极限平衡法中的4种典型计算方法计算了天然工况(考虑重力和地下水)、地震工况(地震加速度取0.125g)和地震工况(地震加速度取0.195g)，计算得到的稳定系数见表2。

图9 安全系数下剖面剪应变增量Figure 9 Section shear strain increment under safety margin

图10 地震动力作用下剪应变增量Figure 10 Shear strain increment under earthquake dynamic action

地震动力数值分析所加的地震波量值上相当于拟静力极限平衡方法中地震力0.195g，对比二种分析方法，可以看出地震动力数值分析得到的稳定性结果明显低于拟静力极限平衡方法，究其原因有两个，一是静荷载与动荷载的作用方式不同会导致结果差异，另一个是滑坡体部位由于动力响应产生加速度的放大作用使滑坡体受到的地震力比底部增大。

表2 极限平衡方法计算成果

图11 地震作用下剖面塑性区分布Figure 11 Section plastic range distribution under earthquake action

图12 折减系数1.05剪应变增量Figure 12 Reduction factor 1.05 shear strain increments

5 结论

1)地震作用下的滑坡体的稳定性比天然工况明显减低，在工程分析中应引起足够重视；

2) 拟静力极限平衡分析法与地震动力数值分析法计算的滑坡稳定性结果差异明显，地震动力数值分析方法更符合滑坡体的受力机理；

3)地震作用下滑坡动力响应分析结果表明，地震动力作用下滑坡体不同部位的动力响应差别明显，一般表现为加速度的放大作用，降低了滑坡体的稳定性；

4)郭家滑坡的地震作用下动力数值分析结果表明，在地震设防烈度地震动力作用下，滑坡体稳定性达不到规范要求，应引起有关部门的重视。