地震荷载作用下堆载影响机场滑坡稳定性的数值计算分析

刘 波

(四川省地质矿产勘查开发局404地质队,四川西昌 615000)

在公路、铁路、水利水电等众多领域建设过程中，经常会遇到各种地形地貌的边坡，边坡在人为因素和自然因素的作用下会发生滑坡，其中堆载是引发滑坡的重要原因之一。目前，国内外已有许多学者对引发滑坡的原因开展了一系列的研究。董夫钱[1-2]等发现堆载不仅直接改变边坡的稳定性，还会影响地下水的渗流条件。李树森[3-4]等对顺层坡发生破坏进行力学分析，建立了破坏的力学模型。陈春利[5]等发现堆载会改变坡体内部应力分布，导致坡体中下部土体剪应力逐渐增大，边坡稳定性随之降低，最终失稳。王洪兵[6]等发现中部及坡脚堆载对边坡稳定有利，坡顶堆载对其稳定极不利。孙红月[7]等发现堆载使坡体产生垂直方向的压缩变形，也使坡体向坡脚方向产生一定的挤出变形。但是，目前同时把堆积荷载和地震荷载一起作为影响边坡变形破坏的研究较少。因此，本文以贵州威宁机场13#斜坡为典型实例，采用数值模拟的方法，对在仅受地震荷载作用和堆载与地震荷载综合作用两种工况条件下斜坡的位移和塑性区进行计算，据此分析其变形破坏规律。

1 滑坡概况

威宁机场13#滑坡所在斜坡位于机场北端近中轴线部位，艳阳沟北东侧上方，斜坡右侧发育一条沟谷，整体地形受沟谷切割呈凸状地形，整体地势西高东低，最高点位于坡顶分水岭处，高程2 454 m，最低处位于冲沟下游斜坡坡脚，高程2 362 m，相对高差约92 m。整个滑坡坡面呈折线型，上下部坡度较陡，中部坡度较缓。滑坡上部大体为一直线坡，坡度约16 °，分布高程2 408.6～2 440.0 m，长约110 m，宽50～105 m；滑坡中部地形相对较缓，整体坡度约9 °，分布高程2 391.1～2 408.6 m，长约90 m，宽105～140 m；下部坡度变陡，约18 °，分布高程2 360.6～2 391.1 m，长约130 m，平均宽度约60 m。

滑坡整体呈SE向展布，滑动主方向为99 °，滑坡平面形态呈长舌型，后缘最宽约105 m，平均宽度约75 m，前缘最宽处约115 m，平均宽度约60 m，滑坡体长度约330 m，平均厚度约10 m，滑坡面积约25 533 m3，总体积约25.5×104m3。

滑坡内岩性主要为二叠系下统梁山组砂岩、炭质页岩和泥岩，石炭纪上统马平群白云质灰岩和第四系粉质黏土及粉质黏土夹碎石。斜坡岩性分为2层，上层以第四系粉质黏土及粉质黏土夹碎石为主，碎石母岩以砂岩为主；下层岩性为强风化泥岩、砂岩、炭质页岩互层，产状紊乱，连续性差。滑床为倾向坡外的中风化炭质页岩及砂岩互层，岩体较完整。根据滑坡体形态和形成条件分析，该滑坡类型为推移式岩质古滑坡。

2 地震荷载模拟方法和Flac3D计算方法

2.1 地震荷载的模拟方法

地震荷载作用下边坡稳定性分析的方法包括动力分析法和拟静力法等。动力分析方法虽然精度较高，但是由于需要较深的专业知识和技能，对工程师素质要求较高且操作复杂，其计算结果的评价也不容易。拟静力法计算方法简单、计算量小、参数易于确定，并已有丰富的使用经验，易被工程设计人员接受[8-9]。

拟静力法就是将地震荷载简化为水平方向和竖直方向不变的惯性力，然后根据极限平衡法和有限元强度折减法计算得到边坡拟静力法安全系数[10-11]。在实际情况下，竖直方向的地震加速度对许多结构的破坏要比水平方向的地震加速度大。因此，本文在考虑地震荷载作用时，只加载了竖直方向上的加速度。

本文运用拟静力法，将地震荷载简化为竖直方向加速度来模拟地震作用对滑坡稳定性的影响。根据G50011-2010《建筑抗震设计规范》和中国地震烈度区划图(1990)，该滑坡区实际地震烈度为6度，设计地震竖直加速度为0.0 637g。

2.2 Flac3D计算方法

2.2.1 三个基本条件

采用Flac3D软件进行数值计算时，必须确定三个部分：

(1)有限差分网格，网格用来定义分析模型的几何形状。

(2)本构关系和材料特性，用来表示模型在外力条件下的力学响应特征。

(3)边界条件和初始条件，用来确定模型的初始状态[12]。完成三个基本条件定义之后，即可进行模型的自重计算，接着施加地震荷载。

2.2.2 Flac3D的优点

(1)采用混合离散法来模拟材料的塑性破坏和塑性流动，比有限元法中的离散集成法更准确、合理。

(2)模拟物理上的不稳定过程不存在数值上的障碍。

(3)采用显示差分法求解微分方程，不形成刚度矩阵且每一时步变形小，可以大大节约内存和时间。

3 数值计算分析

3.1 计算工况

把堆积荷载作为影响斜坡稳定性因素，数值计算可以分为两种计算工况，工况1：自重荷载和地震荷载；工况2：堆积荷载、自重荷载和地震荷载。

3.2 数值模型

13#斜坡一共可以简化为4层介质，由于各个介质空间位置相互联系，边坡可以划分为4个区域。如图1所示，模型X方向为斜坡滑动方向，X正向指向SE方向，Y方向指向滑坡体内部，Z方向指向重力相反的方向。X方向范围为0～50 m，Y方向范围为0～450 m，Z方向范围为2 356～2 448 m。

图1 两种计算工况

3.3 边界条件

Flac3D可以施加包括应力、速度和位移等边界条件。本次数值模拟的边界条件采用位移边界条件，设置如下：底部边界(X、Y和Z)的位移为0，左右侧边界X方向位移为0，Y和Z方向自由，前后侧边界X和Z方向自由，Y方向位移为0。

3.4 计算结果分析

通过Flac3D软件计算边坡模型时，采用弹塑性模型，摩尔-库伦屈服准则，将整个数值模拟分为两步：第一步是施加自重，第二步是施加地震荷载。分别在两种工况下计算边坡的稳定性，分析位移和塑性区。计算结果见图2。

图2 位移计算结果

由图2所示的位移云图可知，在工况一条件下，边坡z向移分布均匀，从上到下逐渐减小，等值线平滑，最大位移发生在坡顶，其值为8.68 mm。总位移分布和Z向位移相似，最大位移也出现在坡顶，其值为8.68 mm。在工况二条件下，整个填筑体都发生了较大位移，其中Z向位移最大值集中在坡体右上角，大小为20.80 mm，坡脚有较小的隆起，大小为7.91 mm；总位移最大值发生在坡脚，大小为25.90 mm。

通过两种工况下位移云图分析可知，斜坡堆载后，在地震荷载作用下，坡体竖向位移最大值由坡顶转到堆积体，总位移最大值由坡顶转到坡脚，且坡脚发生隆起，最大值有明显增大，增大幅度为17.22 mm。

由图3所示塑性区云图可知，在工况一条件下，即坡体未发生塑性变形，处于稳定状态，软弱夹层中零星分布着一小部分塑性区，但并不影响整个滑坡的稳定性。在工况二条件下，坡体出现大面积塑性区，主要分布在软弱夹层中，说明该部分坡体曾达到屈服强度，已经发生剪切破坏，并影响了滑坡的稳定性。

图3 塑性区分布

通过两种工况下塑性区云图的分析，堆载前，坡体未发生塑性变形，处于稳定状态；堆载后，边坡出现贯通的塑性区，导致滑坡位移较大。

4 结论

对两种工况下斜坡稳定性计算分析结果表明：

(1)在地震荷载作用下，堆载后，边坡竖向位移增大且坡体向坡脚方向发生挤出变形，边坡填筑体有较大位移，位移最大值从坡顶转移到填筑体，最大值由8.68 mm变为25.90 mm，位移发生明显增大。

(2)在地震荷载作用下，堆载前坡体未发生塑性变形，堆载后边坡出现明显的贯通塑性区，发生剪切破坏，导致填筑体沉降变形大。