基于重力增加法的边坡失稳破坏全过程模拟

，，， ，

(武汉大学 a.水资源与水电工程科学国家重点实验室；b.水工岩石力学教育部重点实验室，武汉 430072)

1 研究背景

边坡失稳破坏[1-2]是一个渐进变化的过程。在外力作用下，坡体内部应力状态不断发生改变，岩体部分范围出现应力集中，产生局部滑裂面。随着局部裂缝延伸，形成一条贯通坡体的裂缝，最终裂缝上部的岩体破碎滑落堆积成稳定状态。传统刚体极限平衡法由于模型简易、计算高效而广泛运用于边坡稳定分析领域[3-7]，但需要人为预先假定滑裂面的位置和形状，并只能模拟边坡临界失稳状态。

近年来，随着计算机软件技术的不断成熟，有限元法、离散元法等在边坡失稳分析领域中迅速发展。重力增加法作为极限平衡分析的一种工具通常被用来结合有限元法、离散元法评估边坡的稳定性。康亚明等[8]将重力增加法运用于有限元法中，对边坡破坏过程中的位移场进行分析，得到了边坡失稳破坏的最危险滑动面；王述红等[9]通过流形元程序采用重力增加法模拟了块状岩体边坡的大变形破坏过程；Li等[10]将重力增加法应用于RFPA程序中，开发了RFPA-GIM程序，对实际工程中的边坡破坏进行了数值模拟；Utili等[11]将重力增加法应用于离散元法中，模拟了理想化岩质边坡的破坏过程。有限元法虽然能反映岩体在变形过程中的应力-应变关系，但无法模拟岩体在破坏过程中复杂的裂缝演化过程。而离散元法虽然能够模拟边坡岩体的大变形，呈现边坡破坏的动态演化过程，却无法反映边坡岩体破坏之前的连续状态。

为了模拟边坡岩体从连续状态向非连续状态转化的渐进破坏全过程，部分学者开始采用连续-离散耦合分析方法[12-13]。Intrieri 等[14]运用连续-离散耦合分析方法研究了Torgiovannetto di Assisi滑坡的诱发机制和演化过程；Elmo等[15]和Havaej等[16]采用连续-离散耦合程序ELFEN研究了不同岩质边坡的破坏过程；Mahabadi等[17]开发了程序Y-Geo，对比分析了海蚀作用下均质岩体和节理岩体的渐进破坏过程。

本文尝试将重力增加法(Gravity Increase Method，GIM)应用于连续-离散耦合分析方法(Combined Finite-Discrete Element Method，FDEM)中，得到基于重力增加法的连续-离散耦合分析方法(FDEM-GIM)，即在有限元法中引入断裂力学的内聚力模型，将界面单元插入边坡模型的表层岩体中，建立连续-离散耦合的边坡计算模型，采用重力增加法对边坡临界破坏状态进行分析，依据边坡系统总动能突增判断边坡的极限状态，计算得到边坡的安全系数和滑面的位置和形状。以红石岩堰塞体边坡工程为实例，对比了FDEM-GIM与刚体极限平衡法的边坡稳定分析结果，验证了基于重力增加法的连续-离散耦合分析方法进行边坡稳定分析的合理性；通过继续增加重力加速度，模拟了边坡达到临界失稳状态之后的后续破坏过程，得到了边坡的最终滑落堆积方量。

2 基本理论与计算方法

2.1 连续-离散耦合分析方法

本文将内聚力模型[18-19]用于连续-离散耦合分析方法中，模拟岩石材料的开裂过程。假设岩石为胶凝材料[20]，在岩石破坏过程中实体单元只发生弹性变形，损伤和断裂发生在界面单元上。在各实体单元之间插入无厚度四节点界面单元，实现短暂的“连续”效果，为了便于示意，界面单元被拉开,如图1所示。采用考虑单元刚度退化的应力-分离准则来模拟裂缝的产生和扩展。

图1 界面单元与实体单元的共节点方式Fig.1 Mode of common nodes between cohesive elements and elastic elements

考虑到岩石等准脆性材料的破坏大多是由拉断和剪断导致，因此本文采用二次应力准则定义裂缝的起裂，即

(1)

岩石材料抗剪强度的计算采用带拉断的Mohr-Coulomb准则，即

(2)

式中：c为材料的内聚力；φ为材料内摩擦角。

当界面单元开始出现损伤后，加载仍然存在，直至界面单元完全失效，产生裂缝。在裂缝不断扩展时，界面单元的本构关系为

(3)

式中：kn，ks分别为界面单元的法向刚度、切向刚度；δn，δs分别为界面单元的法向应变、切向应变；D为无量纲的损伤因子，当D=1时，界面单元失去承载能力，可将界面单元从岩石材料中剔除。

采用基于线性软化的Benzeggagh-Kenane准则，界面单元的应力-分离曲线如图2，对应公式为

(4)

图2 界面单元应力-分离曲线Fig.2 Stress-separation curves of cohesive elements

2.2 FDEM-GIM边坡失稳破坏分析方法

本文采用大型有限元软件ABAQUS进行边坡失稳破坏过程的数值模拟。为真实反映边坡模型的初始状态，对边坡进行地应力平衡计算。施加重力作用，计算时忽略上部岩体及节理裂隙的开裂行为。

重力增加法常被用于计算边坡安全系数，其基本原理为在控制岩体的抗剪强度参数c和φ不变的前提下，通过逐渐增加重力加速度(即增大自重力作用)，使得边坡处于临界失稳状态，此时对应的重力加速度(即临界重力加速度glimit)与实际重力加速度g0的比值即为边坡的安全系数，即

(5)

式中：glimit为边坡处于临界破坏状态时的重力加速度；g0为边坡实际重力加速度，通常取g0=9.8 m/s2；Fs为边坡的安全系数。

本文边坡临界失稳状态的判别标准是边坡的系统总动能发生突增，其计算公式为

(6)

式中：EK为边坡的系统总动能；ρ，v，V分别为积分点处的密度、速度和体积。

文中将连续-离散耦合分析方法与重力增加法相结合，采用FDEM-GIM分析方法模拟边坡失稳破坏的全过程，具体流程如图3所示。

图3 边坡失稳破坏模拟流程Fig.3 Flow chart of simulation of slope instability

3 工程实例分析

3.1 红石岩堰塞体概况

受云南省鲁甸县发生在2014年8月3日的地震影响，在鲁甸县火德红乡李家山村和巧家县包谷垴乡红石岩村交界的牛栏江干流上，两岸山体大规模塌方形成堰塞湖，震后地形地貌如图4。

图4 地震后地形地貌Fig.4 Topographic features after earthquake

堰塞湖集水面积11 832 km2，堰塞体位于红石岩水电站取水坝下游600 m处，在取水坝与厂房之间。堰塞体顶部高程1 222 m，估算堰塞体总方量约1 200万m3，两岸边坡坡高在800～1 100 m之间。由于震后边坡岩性强度较为薄弱，在外力作用下，可能导致两岸边坡继续失稳滑塌，对坡脚堆积体和工程安全产生不利影响。因此,有必要对两岸边坡的稳定性进行分析，研究边坡失稳破坏的全过程。

图5 剖面g-g边坡计算模型和边坡材料分区Fig.5 Computational model and material partition of section g-g

3.2 计算模型与材料参数

计算选取红石岩堰塞体右岸边坡的一个典型剖面g-g断面，剖面g-g的右岸边坡坡高为1 020 m。进行二维FDEM-GIM边坡失稳分析，计算模型如图5(a)所示。在滑动部分的各实体单元之间插入厚度为0的四节点界面单元。为了提高模拟的精度并减少计算量，仅对上部滑动部分岩体的网格进行了加密。上部岩体与卸荷裂隙网格尺寸为2 m，下部岩体网格尺寸为20 m。剖面g-g右岸边坡模型的实体单元数量为26 376个，界面单元数量为32 775个。根据相应工程资料，材料分区如图5(b)所示，计算模型参数见表1。结合实测资料并通过实验室单轴压缩数值试验反复调整[21-22]，获得表层岩体和节理裂隙的界面单元力学特性参数，见表2。

表1 计算模型参数Table 1 Parameters of computational model

表2 表层岩体和节理裂隙的界面单元力学参数Table 2 Mechanical parameters of cohesive elements of surface rock mass and joint fractures

3.3 边坡临界失稳状态分析

图6 剖面g-g系统总动能历时曲线Fig.6 Curve of time duration of systematic total kinetic energy of section g-g

根据FDEM-GIM计算结果提取得到了剖面g-g边坡模型的系统总动能历时曲线，如图6所示。当t=16.8 s时，剖面g-g的系统总动能发生突增，边坡处于临界失稳状态，对应临界重力加速度glimit为1.680g0(文中取g0=9.8 m/s2)，边坡安全系数Fs=1.680。

为了验证计算结果的可靠性，采用刚体极限平衡法对剖面g-g进行边坡稳定分析，将分析结果与FDEM-GIM的计算结果进行对比。表3列出了FDEM-GIM与刚体极限平衡法求得的边坡安全系数。

表3 剖面g-g右岸边坡模型安全系数Table 3 Factor of safety of slope

由表3知，通过刚体极限平衡法计算得到的边坡安全系数为1.616，与FDEM-GIM所得到的边坡安全系数基本一致。图7为刚体极限平衡法与FDEM-GIM模拟得到的边坡滑裂面对比，当边坡处于临界失稳状态，刚体极限平衡法与FDEM-GIM搜索得到的最危险滑裂面的位置以及形状基本吻合，因此采用基于重力增加法的连续-离散耦合分析方法分析边坡的安全系数是可行的。

图7 刚体极限平衡法与FDEM-GIM边坡滑裂面对比Fig.7 Comparison of the slope slip surface obtained respectively from rigid body limit equilibrium method and FDEM-GIM

图8 FDEM-GIM模拟边坡开裂-滑落-堆积全过程Fig.8 Whole process of fracture, sliding, and deposition of slope simulated by FDEM-GIM

3.4 边坡破坏全过程模拟

图8为FDEM-GIM模拟得到的剖面g-g右岸边坡的失稳破坏动态演化全过程。当t=10.0 s时，边坡在g0下能够维持其自身稳定。随着重力加速度的增大，表层岩体首先从其底部逐渐开裂，当t=16.8 s时，底部出现一条局部贯通裂缝，边坡处于临界失稳状态，此贯通裂缝为坡体的第一滑裂面，是边坡潜在最危险滑裂面。当继续增加重力加速度，使达到2g0之时，底部岩体不断相互碰撞并破碎成若干小块，边坡开始产生一条沿卸荷裂隙延伸的裂缝。当t=32.0 s时，所加重力加速度仍然为2g0，边坡逐渐产生一条贯通坡体的裂缝，边坡岩体在向下运动的过程中，不断撞击破碎成若干岩块，呈现出明显的坍塌破坏模式，最终堆积成稳定状态。

由于边坡失稳破坏是一个渐进破坏过程，当边坡达到极限状态时，继续增加重力加速度会导致上部岩体随之开裂。边坡表层岩体的底部滑落使得上部失去了承托，在后续破坏阶段，边坡沿着第一滑裂面自下而上延伸出一条贯通坡体的裂缝。从图8中边坡滑落堆积的情况可以看出，边坡在出现第一次滑落堆积之后仍然有后续滑落过程。因此FDEM-GIM能够模拟边坡在达到临界失稳状态之后的开裂行为，边坡最终滑落方量为贯通裂缝上部体积，而刚体极限平衡法仅能用来分析边坡临界破坏状态。

4 结 论

本文基于重力增加法，采用连续-离散耦合分析方法对边坡岩体进行失稳破坏的全过程模拟，以红石岩堰塞体边坡工程为实例，将刚体极限平衡法得到的边坡安全系数与FDEM-GIM得到的边坡安全系数进行对比，验证了基于重力增加法的连续-离散耦合分析方法运用于边坡失稳破坏问题的合理性，并模拟了边坡岩体从未发生破坏到裂缝产生、扩展、破坏直至滑落堆积的全过程,结论如下：

(1)将FDEM-GIM计算结果与传统刚体极限平衡法的结果进行对比，结果显示2种方法得到的边坡安全系数基本相同，最危险滑裂面位置与形状基本吻合，验证了基于重力增加法的连续-离散耦合分析方法进行边坡失稳破坏分析的合理性。

(2)边坡加固措施必须依赖于潜在滑动面位置的确定，基于重力增加法的连续-离散耦合分析方法能自动搜索获得边坡临界破坏滑面，克服传统刚体极限平衡法需要预先人为假定滑裂面的缺陷。刚体极限平衡法仅能用于分析边坡临界失稳状态，不能预测后续滑块的形成。而基于重力增加法的连续-离散耦合分析方法能模拟边坡破坏的全过程，并且可以显示模拟后续滑坡体的形成，为边坡防治处理提供可靠依据。