基于数值流形法的反倾层状岩质边坡倾倒破坏分析1

王欢欢 郭明珠

（北京工业大学，城市建设学部，北京 100124）

引言

反倾层状岩质边坡指岩层走向与边坡走向接近一致，而倾向与坡面倾向相反的边坡，倾倒变形破坏是岩质边坡典型破坏失稳形式（李明霞，2015）。倾倒破坏模式通常分为块体倾倒、块体弯曲倾倒和弯曲倾倒（Goodman 等，1976）。青藏高原地区是板块运动最活跃的地区，该地区地壳内、外动力条件强烈交织和转化，使大型滑坡灾害发育（黄润秋，2007），具有灾害类型多、影响范围广、破坏程度大等特点，其中滑坡、泥石流、崩塌产生的直接危害最严重。青藏高原地区反倾层状岩质边坡深层倾倒现象较常见，特别是金沙江、岷江等深切峡谷地区，倾倒深度可达200～300 m，对川藏铁路、水资源开发等重大工程建设产生影响。反倾层状岩质边坡倾倒破坏严重影响山区经济社会发展，危及群众生命财产安全，因此，需研究反倾层状岩质边坡倾倒破坏特征。

目前，国内外学者对反倾层状岩质边坡的数值模拟分析可分为以下两类：一类是将边坡看作连续体进行数值计算（如极限平衡法、有限元法等），是对真实情况的数值近似表征，分析结构应力与变形后，再进行边坡稳定性判断，但对裂缝等不连续问题和结构破坏后运动过程的处理较困难。韩贝传等（1999）、王宇等（2013）、代仲海等（2018）、王霄等（2018）构建了反倾层状岩质边坡有限元模型，探讨了坡体变形破坏机制及稳定性影响因素。另一类是将岩石看作离散体进行数值计算（如离散元法、非连续变形分析方法），以不连续力学为基础，将岩体结构看作离散的块体，块体之间采用接触连接，不仅可分析块体运动、变形、破坏等，还可对块体间的相互作用进行分析。然而，该方法采用全一阶位移近似和小角度假定，且在块体接触、参数取值、边界处理等方面存在不足，影响含裂缝的不连续结构计算结果准确性。刘云鹏等（2012）、王章琼等（2014）、陶志刚等（2017）等使用UDEC 分析了层状岩质边坡破坏模式。

1992 年，石根华首次提出数值流形法（NMM），充分发挥以上两类方法的优势，解决了有限元、DDA和解析法存在的计算问题，可同时处理连续和非连续问题，具有广泛应用前景。NMM 以数值流形为核心，在非连续变形方法的块体系统理论基础上，融入有限元法和解析法，是连续分析方法。因此，该方法既可计算块体内部变形，也可模拟块体系统力学行为（裴觉民，1997）。

针对反倾层状岩质边坡的研究多采用有限元、离散元或非连续变形分析方法等,采用NMM 进行地震力作用下的边坡破坏研究较少。本文采用NMM 建立反倾边坡模型，对地震作用下不同内摩擦角的反倾层状岩质边坡进行倾倒破坏模拟，并分析倾倒破坏过程。

1 边坡概况

滑坡区主要出露石英片岩、砂岩等硬岩及少量片麻岩和石灰岩等软岩，岩体结构面发育。石英片岩是坡内出露的主要岩性，成层性较好，表面风化强烈。堆积体主要成分为石英片岩夹少量片麻岩，区域岩性变化较小。滑坡堆积体方量4 500 万m³，水平向长度1 400 m，高度方向落差1 000 m，滑坡堆积体距坡趾最远达1 022 m。滑坡区主要受3 组节理控制，倾向分别为375°∠40°、78°∠82°和185°∠78°，岩层表现为反倾向。通过调查分析，绘制贡扎滑坡体平面图和主剖面图，分别如图1、2 所示。

图1 贡扎滑坡平面图Fig. 1 Plan of gong-za landslide

图2 贡扎滑坡主剖面图Fig. 2 Main section of gong-za landslide

2 室内力学试验

点荷载试验是获取岩石力学参数的主要方法之一，在进行室内点荷载试验的基础上，通过仪器油表获得岩石点荷载强度。为确定数值模拟主要物理力学参数，进行室内点荷载试验。由于对贡扎滑坡进行了现场实地勘察，并从现场带回了多块堆积体岩石（图3），因此开展室内点荷载试验较方便，试验数据具有可信度。

图3 试验所用硬岩和软岩Fig. 3 Hard and soft rocks used in the test

根据点荷载计算原理，进行径向测试时，等效岩芯直径De即为岩芯直径D。进行不规则岩块轴向测试时，等效岩芯直径De计算如式（1）所示。当De取50 mm 时，可换算出直径为50 mm 标准试件点荷载强度，由于原型岩石形状不规则，故采用式（2）、（3）进行计算。通过点荷载强度Is（50）及式（4）可推导出岩石单轴抗压强度值σc及单轴抗拉强度στ，进而推导出黏聚力c及内摩擦角φ。

式中：W为两加荷点最小截面宽度（mm）；F为修正系数；Is为未修正点荷载强度（MPa）；m为修正指标，一般取值为0.45；k为修正系数，当岩石为硬岩时k=0.8～0.9，当岩石为软岩时k=0.6～0.7。

计算得到贡扎滑坡地区硬岩黏聚力c=17.23 MPa，内摩擦角φ=47.39°；软岩点荷载强度Is（50）=0.74 MPa，单轴抗压强度σc=14.8 MPa，单轴抗拉强度στ=1.11 MPa，黏聚力c=2.02 MPa，内摩擦角φ=41.57°。

3 算例分析

3.1 计算模型

为简化数值模拟计算过程，采用图4 所示简化边坡模型，其中坡角为30°、岩层倾角为65°、边坡长1 400 m、高1 200 m。模型边界固定采用位移边界固定法，对左右边界及底部进行约束。各岩层均采用莫尔-库仑准则岩石本构模型，同时，在边坡上布置3 组编号为1～6 的监测点，以便分析内摩擦角对监测点位移的影响。

图4 简化边坡模型Fig. 4 simplified model of slope

3.2 模型参数选取

坡体主要包括石英片岩、砂岩及少量片麻岩和石灰岩等，石英片岩是坡内出露的主要岩体，成层性较好，表面风化强烈。堆积体主要成分为石英片岩夹少量片麻岩，通过室内试验得到相关力学参数如表1 所示。动力系数为1，单步允许最大位移率为0.01，总时间步选为10 000 步，泊松比为0.25，超松弛系数为1.35。

表1 岩石力学参数取值Table 1 Rock mechanical parameters

3.3 数值计算过程

主要研究内摩擦角φ对岩质边坡倾倒破坏的影响，取计算总时间步数为10 000 步，时间步长为0.008 s，最大位移比为0.001。计算模型如图5 所示，在边坡底部输入日本岩手地震记录到的地震波（图6）。

图5 边坡网格划分模型图Fig. 5 Grid division model of slope

图6 加速度时程曲线Fig. 6 Acceleration time history curves

NMM 计算步骤如下：（1）对含节理、裂缝的岩体进行流形法网格剖分，定义物理和数学网格；（2）进行第1 步加载，计算地震作用下的应力分布及变形；（3）计算各裂缝尖端应力强度因子；（4）读取每个流形单元数据并计算其接触方式；（5）判断接触方式并计算变形；（6）形成法向弹簧和剪切弹簧的加法矩阵，进行迭代求解；（7）计算变形模式。

模拟结果如图7 所示，由图7 可知，随着内摩擦角的增大，块体之间的摩擦系数增加，块体抗滑力增加。因此坡体从开始破坏到新的平衡状态和达到最大位移所需的时间越短。随着内摩擦角的逐渐增大，摩擦系数增大，当增大至一定程度时，阻止块体滑动的力将大于块体滑动力，这时斜坡上的岩石块体处于稳定状态，如图7(e)所示。

图7 内摩擦角不同时边坡破坏特征Fig. 7 Failure characteristics of slope with different friction angles

3.4 模拟结果分析

通过观察数值模拟结果可知，不同内摩擦角下边坡最终破坏状态不同，但边坡上部均开始出现裂缝，然后上部破坏，继而推动边坡下部岩体滑动，导致边坡倾倒破坏。使用Origin 软件进行数据处理，得到监测点位移如图8 所示。

图8 监测点位移图Fig. 8 Displacement diagram of monitoring points

内摩擦角<40°时，大部分岩块产生整体滑移，靠近坡顶处的岩块发生轻微转动，推动前面的岩块加速滑动，呈倾倒-滑移模式。内摩擦角＞40°时，坡顶岩块首先产生滑动并转动驱使前面的岩块滑动，推动坡脚处岩块产生滑动，呈渐进式倾倒破坏，产生整体性破坏的可能性较小。

由图8 可知，前15 s 滑坡体整体移动，说明在水平地震作用下整体下滑，位移变化差异小；然后下部坡体在水平力作用下抛出，中间的倾倒区域块体出现倾覆现象，最终坡体倾倒破坏。坡脚位置的监测点5、6水平位移远大于其他位置，说明在地震来临时，边坡下部的岩体对地震作用更敏感，移动距离较大。

4 结语

采用二维数值流形法分析了在其他岩体物理力学参数保持常量时，仅改变内摩擦角，地震作用下反倾层状岩质边坡倾倒破坏特性。结果表明，内摩擦角<40°时，坡体首先在水平地震作用下整体下滑，然后下部坡体在水平力作用下抛出，中间的倾倒区域块体出现倾覆现象，呈倾倒-滑移模式；内摩擦角>40°时，少量岩块发生倾倒破坏，边坡较稳定。