基于RFPA3D不同加载方法的边坡渐进破坏分析

康佳辉 ，刘晓辉，张永彬，向柏宇

（1.西华大学流体及动力机械教育部重点试验室，四川 成都 610039；2.四川航天建筑工程有限公司，四川 成都 610100；3.力软科技（大连）股份有限公司，辽宁 大连 116600；4.中国水电工程顾问集团成都勘测设计研究院，四川 成都 610072）

随着我国高新技术产业的飞速发展，社会对电能的消耗与日俱增。特别是我国西南地区，水力资源丰富，是我国水电开发的重要基地[1]。近年来，我国在建及已建的大型水电工程有白鹤滩、锦屏一级、溪洛渡、二滩等。这些水电工程多分布在地质条件复杂的深山峡谷之中，坝顶两侧还有高陡的自然边坡。在长期地质构造运动的影响下，两岸边坡内部断层发育，高边坡稳定性问题十分突出。此外，这类边坡的平均坡度大多超过40°，边坡受断层等软弱结构面切割严重，边坡稳定性难以控制，边坡失稳必然对水电建筑物造成巨大破坏。表1 统计了部分我国在建及已建高陡边坡的开挖高度。由表1 中数据可知，边坡开挖高度最大已超过500 m，超高边坡及特高边坡不断出现，成为工程建设成败的关键，影响着水电站的安全运行和效益，对高边坡稳定性问题及潜在失稳模式进行模拟、研究、预警显得十分必要。

表1 高陡边坡最大开挖高度示例

数值模拟是研究高边坡问题最常用的手段，随着计算机技术的进步更是得到了迅猛的发展。有限元法[2]、边界元法[3]、离散元法[4]、快速拉格朗日法[5]、DDA[6]、RFPA[7]、SPH[8]、XFEM[9]等方法在高边坡问题的分析中得到充分运用和发展。郑颖人等[10]较早在强度折减法领域取得丰硕成果，随后赵尚毅等[11−12]将强度折减理论引入到有限元分析中，解决了有限元在边坡稳定性分析中的应用问题。陈国庆等[13]提出基于动态和整体强度折减的边坡动态稳定性评价方法，并通过算例进行了验证。王庚荪[14]针对边坡破坏过程的渐进性，提出了新的接触单元模型来模拟边坡滑面上的接触摩擦状态，取得不错的效果。上述方法主要基于有限元分析，但有限元方法在处理大变形、大位移问题时仍然存在缺陷，常常会因为计算不收敛而终止，所得结果也与实际情况不符。钟波波等[15]采用RFPA-Centrifuge 研究了不同岩桥倾角下的岩石裂纹扩展模式。刘欣欣等[16]将DDD 与PRFA 离心加载法结合，研究了黑山铁矿边坡的滑移机制。刘佳等[17]采用RFPA 软件模拟了露天高边坡的失稳破坏模式，研究了边坡高度与角度对破坏模式的影响。唐春安等[18]采用RFPA-SRM 分析了软硬互层边坡的变形破坏过程，对比研究了RFPA-SRM 在边坡稳定性分析中的可行性，认为RFPA 在处理连续介质力学问题时更符合实际。杨莹等[19]等采用RFPA-Centrifuge 对白鹤滩电站左岸边坡进行了稳定性分析，通过RFPA 独有的声发射监测数据与微震监测进行了对比分析，研究发现RFPA 模拟的边坡渐进破坏过程与微震监测结果基本一致，说明RFPA 能很好模拟边坡在渐进破坏过程中的变形破坏规律。虽然以上两种方法在分析边坡稳定性问题上都有所应用，但对其适用条件的研究则相对较少，特别是对于大型水电工程边坡而言，采用不同计算方法所得计算结果差距较大，选择合适的计算方法将在保证工程安全的前提下降低工程成本，使工程更加经济可靠。

基于此，本文以锦屏一级左岸边坡为例，采用RFPA3D离心加载法和强度折减法再现边坡渐进失稳过程，分析其变形破坏过程中的应力分布和能量释放规律，模拟边坡在自重应力条件下岩体微裂纹的萌生、扩展、贯通直至破坏的全过程，分析边坡潜在破坏模式。

1 计算方法

RFPA3D基于岩石损伤本构关系，同时考虑了岩石材料的缺陷分布随机性，以此来模拟岩石材料强度的离散性。表2 给出了RFPA3D强度折减法和离心加载法的计算原理。表2 中：Step为单元破坏数最大时对应的加载步数；Δg为 离心加载系数；γ材料容重；Fs为 折减系数；s tep为强度折减中边坡失稳时的加载步。

表2 RFPA3D强度折减和离心加载计算原理

可以看到，两种方法所采用的本构关系是一致的，均为弹脆性本构关系，单元在达到破坏准则前始终满足弹性力学特性。安全系数的定义根据加载方式的不同而有所区别，RFPA-SRM 安全系数与一般强度折减法的安全系数的物理意义相同，即以边坡刚好处于临界失稳时的折减系数作为其安全系数。RFPA-Centrifuge 则以模型失稳时的单元自重与初始单元自重的比值作为其安全系数。

Griffiths 等[20]认为在边坡稳定性分析中适当考虑材料的拉伸破坏是必要的，特别是在岩质边坡稳定性分析中，拉伸破坏有时是构成边坡失稳的最关键因素[20]。所以，本文引入了带拉伸-截断的Mohr-Coulomb 准则[21]，充分考虑岩石材料的拉伸破坏，使其与实际情况接近，带拉伸截断的摩尔库伦破坏准则如图1 所示。

图1 带拉伸截断的摩尔−库伦破坏准则

图中：Ⅰ区为压缩引起的剪切破坏，Ⅲ区为拉伸引起的拉伸破坏，Ⅳ区剪切和拉伸复合破坏，相关参数关系为：

式中：σc为 Ⅰ类相变阈值；σt为 Ⅱ类相变阈值；φ为内摩擦角；λ为拉压比；σ1为 最大主应力，σ3为最小主应力。

可以看到，在RFPA3D中软件会首先考察单元的拉伸破坏，如果单元的最大拉应力先达到材料的拉应力阈值，则单元发生拉伸破坏；否则在单元没有发生拉伸破坏的前提下，同时细观单元的应力状态满足Mohr-Coulomb 准则，则单元发生剪切破坏。同时，软件根据设定的残余强度系数，认为破坏后的单元仍具有一定的承载能力。

2 工程背景

锦屏一级坝址区两岸山体高陡，谷坡陡峭，基岩出露，相对高差千余米，为典型的“V”字型河谷，是天然的建坝区域。左岸边坡岩层走向与河流流向基本平行，为反向坡，如图2 所示。

图2 锦屏一级水电工程全貌

地质调查表明，左岸坝头地貌上为两山梁夹一楔形凹槽，地形坡度在40°～50°范围内。地质勘探表明，左岸上覆岩体为砂板岩，下覆基岩为大理岩，坡内发育有近SN(N20°W)向的深部裂缝条带，缓倾F42-9 断层(近EW，S∠40°～60°），陡倾煌斑岩脉X(N60°～80°E，SE∠70°～80°）和断层F5、F8，地质特征参数如表3、表4 所示。经初步判定，左岸F5、F8、F2、F42-9、煌斑岩脉X 可能为边坡主控断层，在边坡稳定性分析中应重点考察。据此，本文在以锦屏一级左岸边坡为计算模型的同时，重点分析边坡主要断层对边坡稳定性及破坏模式的影响。

表3 左岸边坡岩性分布

表4 左岸边坡主要断层及层间挤压带

3 模型及参数

本文以左岸边坡沿拱轴线剖面建模，所采用的有限元单元类型为Solid 45，即三维八节点单元，该单元类型在模拟塑性及大变形位移时有较好的适应性。

模型计算范围顺河向820 m，横河向906 m，垂向852 m。模拟范围包括4 类岩体（Ⅱ1、Ⅲ1、Ⅲ2、Ⅳ1）和3 条主要断层F2、F5、F8 以及煌斑岩脉X。参考相关资料[22]，确定岩体物理力学参数如表5所示。

表5 岩体物理力学参数

建模时先通过工程地质剖面图经适当简化后拉伸为三维实体模型，再采用参数化建模软件ANSYS APDL 进行网格划分、边界条件设置及材料参数设置，最后通过数据转化接口Ator，将ANSYS 中的节点、单元、边界条件等信息提取并导入RFPA3D中进行计算。

有限元模型坐标系X轴正向指向右岸坡内，Y轴正向指向河流上游，Z轴正向竖直向下。共划分六面体单元数1 390 500 个，节点数1 432 951个，计算域四周法向约束，底部完全固定，地表自由。有限元实体模型、网格划分及约束情况如图3 所示。

图3 计算模型及网格划分

两种计算方法所采用的加载系数均为0.01，即离心加载法按单步增量为0.01 进行加载至边坡失稳，强度折减法按单步增量为0.01 进行折减至边坡失稳。本文暂未考虑边坡在开挖和锚索加固条件下的工作状态，仅重点关注边坡主要断层、岩脉等软弱结构面对边坡稳定性及变形破坏情况的影响。

4 结果分析

4.1 渐进破坏过程分析

图4 为边坡渐进破坏过程，其中图4（a）—图4（d）为RFPA3D-Centrifuge 计算结果，图4（e）—图4（h）为RFPA3D-SRM 计算结果。可以看到，离心加载条件下，边坡的渐进破坏过程为：初期（图4（a））边坡在横河向和竖直向的位移变形值很小，未出现宏观破裂；加载中后期（图4（b））边坡内部进行应力调整重分布，内部软弱结构面被压密，但结构面未出现滑动；加载后期，边坡开始出现较密集破坏，当加载到76～3 步时，拉裂缝沿下游方向贯通，黄斑岩脉X 向下贯通，边坡出现整体沿煌斑岩脉X 向下滑移现象，F5 外部岩体出现碎裂破坏，F5、F8 断层失去抗剪断能力，出现斜向拉裂缝贯穿至坡顶；当加载到76～5 步时，黄斑岩脉X 以外岩体完全破坏，边坡以F5、F8、煌斑岩脉X 为滑裂面，向下滑移剪断F2 断层整体失稳破坏，此时计算得到的安全系数为1.75。

强度折减条件下，边坡的渐进破坏过程为：折减初期（图4（e）），边坡位移变形值较低，煌斑岩脉X 高高程区域出现部分破坏；当折减到第25～0 步时，黄斑岩脉X 沿横向、斜向扩展，边坡开始出现宏观破坏；当折减到34～7 步时，煌斑岩脉X、断层F8 向下完全贯通，断层F5 横向扩展，F8、F5 于交汇处贯通，向下剪切F2 断层，新生裂隙沿F8、F5 断层于交汇处沿斜向上扩展；最终边坡在折减至34～11 步时失稳破坏，此时计算得到安全系数为1.52。

图4 边坡渐进破坏过程（单位：MPa）

从计算揭示的边坡渐进破坏过程可以看出，边坡在水平方向上表现为沿指向河谷方向上的变形，竖直方向上表现为沿竖直向下的变形。进一步分析表明，边坡在煌斑岩脉X、F8、F5 交汇处最先破坏，出现应力集中，而后裂隙沿横向和斜向上扩展、贯通至边坡失稳破坏，说明此处是影响边坡潜在失稳的关键部位，同时也印证了软弱结构面是控制岩质边坡破坏模式的主要因素之一[23]。

综上所述，边坡主要破坏模式表现为以煌斑岩脉X 为后缘面，以岩层分界线为底滑面，剪断Ⅳ1和Ⅲ2类岩体的剪切滑移模式和以F5 为后缘面，以岩层分界线为底滑面，剪断Ⅳ1和Ⅲ2类岩体形成的剪切滑移模式。该结论与王如宾等[24]采用西原模型及有限元法对锦屏一级左岸边坡长期稳定性变形分析的结论是一致的。

表6 为不同计算方法得到的锦屏一级左岸边坡在天然状态下的安全系数。周钟等[25]采用三维刚体极限平衡法计算得到的安全系数为1.202；许强等[26]采用FLAC3D 强度折减法计算得到的安全系数为1.277；漆祖芳等[27]采用ANSYS 强度折减法计算得到的安全系数为1.45，文献[28]中采用PHASE2和FINAL 计算得到的安全系数分别为1.70 和1.75。可以看到，不同计算方法得到的安全系数差异较大，但报告中FINAL 软件和PHASE2计算得到的安全系数与本文采用RFPA3D-Centrifuge计算得到的安全系数较为一致。另外，从模拟得到的边坡渐进破坏过程来看，RFPA3D-Centrifuge 模拟揭示的边坡失稳模式与徐奴文等[29]采用微震监测对锦屏左岸边坡的分析结果及报告的结论较为一致，说明采用RFPA3D-Centrifuge 法模拟边坡破坏模式更为合适，但结果偏于安全。

表6 不同计算方法安全系数

4.2 声发射事件分析

岩石声发射（AE）可定义为岩石在变形过程中，由于微破裂（微损伤）使聚集的能量突然释放而产生的一种弹性波，显然那些微破裂是岩石内部微损伤的结果[30]。也就是说，声发射活动代表了岩石微损伤程度。岩质边坡的失稳破坏通常是由岩体内部微破裂、微损伤的产生和不断聚集并逐渐累积而形成的。RFPA 自带的声发射监测功能可对岩石从最初的裂纹萌生直至宏观破坏进行实时记录，进而可以分析岩质高边坡在整个失稳破坏过程中的能量释放规律。

图5 为采用RFPA3D-Centrifuge 计算方法得到的边坡渐进破坏过程中声发射图。需要指出的是，在声发射图中，蓝色表示拉伸破坏，红色表示剪切破坏。由图5 可知，边坡上部主要为拉破坏，坡脚主要为剪切破坏并伴随少量的拉破坏。在离心加载条件下边坡的声发射情况在加载初期十分微弱，仅坡顶部分岩体力学参数较低的区域出现破裂，边坡并未出现宏观失稳破坏。加载中后期（图5（b）），煌斑岩脉X 上部出现较为密集的声发射现象，断层F2 出现微破裂并扩展；至76～3 加载步，煌斑岩脉X 完全贯通，断层F5 出现大量剪切破坏，声发射现象强烈，边坡岩体出现大面积损伤；最终黄斑岩脉X、断层F5、F8 剪断F2 释放巨大能量，边坡失稳破坏。

图5 离心加载条件下边坡破坏过程声发射图

图6 为离心加载过程中声发射事件数变化过程。需要指出的是，由于边坡在最大破坏步时的声发射事件数十分剧烈（为1 332 010 个），较前一步增加了近1 200 倍，故仅取起始加载步至最大破坏前一步的声发射事件数以作统计分析。可以看到，边坡在整个失稳破坏过程中能量释放巨大，接近边坡失稳时的声发射事件数密集而强烈。在68 加载步后，声发射事件数开始剧烈增加，此时边坡已累积损伤89%。边坡在整个失稳破坏过程中的能量释放是一个逐渐累积的过程，能量释放主要集中在加载后期，这与图4 的渐进破坏过程表现相一致，再一次验证了边坡破坏过程的累进性。

图6 声发射事件数变化过程

图7 为强度折减过程中得到的声发射图，与离心加载初期未出现微破裂不同的是，折减初期，黄斑岩脉X 上部、F2 下部出现微破裂和微损伤并逐渐扩展；折减中期（图7（b））煌斑岩脉X 外部、F5 断层内部岩体区域完全弱化损伤；折减至34～7 步时，F5、F8、黄斑岩脉X 完全贯通，三者交叉处产生大量剪切破坏，声发射密集；折减至34～11 步时，浅层边坡岩体完全损伤，形成以黄斑岩脉X 和断层F5、F8 为上滑面、以F2 为下滑面的滑动块体，此时声发射最为强烈，边坡整体损伤严重。

图7 强度折减条件下边坡破坏过程岩体声发射图

图8 为强度折减过程中得到的声发射事件数变化过程。同样地，由于边坡在最大破坏步时的声发射事件数十分剧烈（为1 222 629 个），较前一步增加了近1 100 倍，故仅取起始加载步至最大破坏前一步的声发射事件数以作统计分析。与离心加载结果不同的是，边坡在整个渐进失稳过程中声发射现象较活跃。折减前期，声发射活动水平相对平稳但总体水平高于采用离心加载法所得结果，说明边坡在折减条件下的损伤更为严重，这也印证了采用强度折减法计算得到的安全系数比采用离心加载法得到的安全系数偏低的结论。至折减中期，声发射事件数开始出现较为剧烈增长，折减后期（第27 折减步）时，声发射事件数剧烈增长，此时边坡已累积损伤82%。

图8 声发射事件数变化过程

综上所述，通过声发射图可以看到，边坡在加载初期就已有声发射现象，但十分微弱；随着加载进行，微破裂微损伤逐渐扩展，边坡已有的软弱结构面开始出现滑动和拉破坏，裂隙面开始延展，边坡声发射数逐渐增加，声发射现象逐渐活跃；加载后期，边坡能量释放十分剧烈，软弱结构面彻底贯通，薄弱岩层被剪断，边坡失稳破坏。从声发射事件数变化过程可以看出，边坡从初始加载至失稳破坏过程中的能量释放经历了初始增长−稳定发展−加速发展−剧烈增长−失稳破坏五个阶段。

5 结论

本文采用同一模型，相同参数，不同计算方法来分析锦屏一级左岸高边坡的失稳破坏过程，并与多种计算方法进行对比，得到以下结论：

1）采用RFPA3D-Centrifuge 和RFPA3D-SRM 计算方法进行边坡稳定性分析是可行的。通过建立三维有限元模型，编制Ator 数据转化接口，采用RFPA3D软件可实现大型岩质边坡原位地质模型的模拟，再现边坡失稳破坏过程中岩体微裂纹、微破裂的萌生、发育、扩展、相互作用直至完全破坏的全过程。

2）计算结果表明，边坡在天然状态下具有一定安全裕度。边坡的潜在滑移模式可能为以煌斑岩脉X 为后缘面，以岩层分界线为底滑面，剪断Ⅳ1和Ⅲ2类岩体的剪切滑移模式和以F5 为后缘面，以岩层分界线为底滑面，剪断Ⅳ1和Ⅲ2类岩体形成的剪切滑移模式。

3）渐进破坏过程揭示出F8、F5、F2 以及黄斑岩脉X 等软弱结构面是控制边坡失稳模式的主要因素。边坡声发射事件分析表明，锦屏一级左岸边坡从初始加载至失稳破坏过程中的能量释放经历了初始增长—稳定发展—加速发展—剧烈增长—失稳破坏5 个阶段。