基于DEM-LEM考虑裂隙扩展的岩质边坡稳定性分析

曾海艳，吴玉龙，张建海

(1.四川省水利水电勘测设计研究院， 四川 成都 610072； 2.成都市市政工程设计研究院， 四川 成都 610023；3.四川大学 水利水电学院 水力学及山区河流开发与保护国家重点实验室， 四川 成都 610065)

岩质边坡内部通常存在许多复杂的节理裂隙，在外荷载作用时，内部的节理裂隙将有可能出现应力集中、裂隙扩展、裂隙面拉裂以及边坡滑移破坏等现象。因此，分析岩质边坡稳定性时考虑其内部潜在的滑移块体十分必要。目前岩质边坡裂隙扩展问题常用的方法主要是地质力学模型试验法和有限元法、无单元法、扩展有限元法、离散元法等数值计算方法。地质力学模型试验法[1]可以比较准确的模拟岩质边坡主要结构以及其扩展对边坡稳定性的影响，但模型材料组合复杂，加工成本较高，且弹性模量、强度参数的测试比较困难。有限元法一般采用裂隙尖端网格加密或者设置奇异单元进行数值模拟，该方法易出现单元扭曲变形导致求解误差加大、计算终止等问题[2]。无单元法[3]在模拟裂隙扩展时，只需引入少量节点并在局部进行必要的人工介入即可自动完成裂隙路径的确定与扩展而得以广泛应用[4]。但该方法随积分域的选取产生不稳定的计算结果，同时由于近似函数积分复杂从而严重限制其应用。扩展有限元法[5-6]除了拥有广义有限元法的优点，还考虑了引用其他方法确定裂隙的实际位置，跟踪裂隙生长，但难以处理复杂裂隙分割的单元，且附加函数难以确定，对于大变形失稳问题的求解存在一定的局限性。离散元法[7]在模拟计算模型时能反映块体之间接触、滑移、分离和倾翻等大位移，但经典的离散元程序在模型前处理时会自动删除无法形成封闭回路的裂隙，难以模拟裂隙起裂、扩展及失效过程。

刚体极限平衡法是经典的力学分析方法，作为较早的岩质边坡计算方法之一，发展相对成熟，且同时又有一套与之相应的规范，因此得到了广泛的应用。但由于刚体极限平衡法需事先假定危险滑动面，故难以应用于复杂裂隙岩质边坡的计算。为解决内部存在复杂节理裂隙的岩质边坡稳定判定问题，作者综合离散元法和刚体极限平衡法的优势，提出一套基于离散元法及刚体极限平衡法(DEM-LEM)并能够考虑裂隙扩展的岩质边坡稳定分析方法，同时研制开发了岩质边坡稳定分析系统SAFESLOPE。

1 计算原理

图1表示的是SAFESLOPE分析系统的基本原理，该方法首先需建立岩质边坡的原始计算模型，利用虚拟裂隙生成算法获取边坡节理裂隙的扩展模型文件。将扩展模型导入离散元程序UDEC中，计算裂隙扩展模型的应力场并获取裂隙尖端的应力值。根据格里菲斯准则判定尖端是否沿虚拟裂隙起裂扩展，若满足扩展条件则计算裂隙扩展长度，然后更新计算模型并改变虚拟裂隙面力学参数，用UDEC重新计算扩展模型的应力场，进而进行滑块重组，同时计算该时步各滑移块体的稳定安全系数。重复上述步骤，直至裂隙扩展完成。最后根据具有最小安全系数的组合块体，确定岩质边坡关键滑移块体及关键滑移路径。

图1 SAFESLOPE稳定分析系统基本流程图

1.1 虚拟裂隙生成算法

如图2所示，虚拟裂隙的形成是以悬枝裂隙末端点为圆心，以固定增量的半径(r=r+Δr)画圆，直至找到最近的边界。增设的虚拟节理裂隙主要分为悬枝裂隙末端点向临近悬枝裂隙末端点扩展的虚拟裂隙(Ⅰ类虚拟裂隙)、悬枝裂隙末端点向计算边界扩展的虚拟裂隙(Ⅱ类虚拟裂隙)、悬枝裂隙末端点向临近真实裂隙扩展的虚拟裂隙(Ⅲ类虚拟裂隙)。最终所形成的三类虚拟裂隙其力学参数不发生改变。

1.2 虚拟裂隙扩展及长度计算

对已经确定虚拟节理裂隙的扩展计算模型施加动荷载作用，获取相应时刻岩质边坡的应力场，同时利用不连续位移法(DDM)[8]计算悬枝裂隙尖端扩展点的应力值。

图2虚拟裂隙类型及生成方式

根据Griffith理论[9]，节理裂隙尖端的应力集中现象在受到荷载作用时效应放大，真实裂隙随着荷载的改变沿着虚拟裂隙扩展，判断裂隙尖端是否扩展的依据为公式(1)[10]。

K1≥KIC

(1)

式中：K1为应力强度因子，它反映了裂纹尖端奇异性的强度；KIC为断裂韧度，是材料对裂隙扩展阻力的度量，是一个通过实验确定的常数；某一岩块中悬枝裂隙的长度为α；悬枝裂隙尖端受σt的有效应力；Y为几何修正因子，本文取为1.0。

裂纹扩展速度利用公式(2)进行计算[11]。如果悬枝裂隙尖端满足裂隙扩展判据，则根据动力计算时，每个时步的步长和裂纹扩展速度可计算得该时步虚拟裂隙的扩展长度。

lgV=b+dK1

(2)

式中：V为裂纹扩展速度；b、d为常数(通过实验确定)。

若整个计算过程中悬枝裂隙尖端应力值始终达不到扩展起裂值，则扩展累积位移为零，虚拟裂隙的力学参数不作改变，亦不必更新计算模型。若计算过程中存在裂隙扩展，则裂隙的力学计算参数为真实扩展部分与虚拟扩展部分力学参数的加权平均，更新计算模型并计算边坡稳定安全系数。

1.3 滑块组合算法

由于潜在滑块的不稳定可能导致附近块体滑动，从而致使滑移区域扩大形成组合块体，因此搜索潜在滑块周围的可能滑动块体并组合成新的滑移块体十分必要。图3给出了滑移块体组合算法的流程图，邻近块体为至少一条边与潜在关键块体边界重合的块体，只要邻近块体的所有边不处于约束边界上，即认为该块体可以滑动。

图3滑移块体组合算法流程图

1.4 安全系数的计算

对岩质边坡进行动力稳定分析时，每一个计算时步均可获得块体的应力场，通过对块体积分获取其总阻滑力FZ(t)和总滑动力FH(t)[12]，利用刚体极限平衡法[12-13]计算边坡在动荷载作用下的稳定安全系数(如公式3)。

(3)

式中：FZ(t)、FH(t)为t时刻块体的总阻滑力和总下滑力；r(t)为t时刻块体所受的合力；S为滑动方向(滑动方向的确定参考文献[14-15])；Kd为t时刻滑块瞬时安全系数。

2 案例分析

为了说明裂隙扩展对岩质边坡稳定性的影响，此处假设了图4(a)所示存在悬枝裂隙的岩质边坡简化模型。当两个未贯通的节理裂隙末端较为靠近时，悬枝裂隙尖端将形成Ⅰ类虚拟裂隙，即生成图4(b)所示含潜在关键滑块的计算模型。如图4所示，裂隙扩展前计算模型仅存在一个可移动滑移块体，裂隙扩展后可移动块体增加到两个，其中块体2为潜在关键滑移块体。按照块体组合算法，扩展模型的稳定性评价涉及到块体1、块体2及其组合块体的稳定安全系数计算。表1给出了岩质边坡物理力学参数[16]，其中岩体的剪切模量G为0.3 GPa，泊松比μ为0.25，断裂韧度KIC为1.0 MPa·m1/2，常数b=-38.5，d=2.46×10-5。另外，虚拟裂隙物理力学参数取为已扩展为真实裂隙部分力学参数与未扩展部分力学参数的加权平均值。

图4 含未贯通节理边坡模型(m)

根据计算时考虑裂隙扩展与否，将采用两个工况进行计算分析(如表2)。计算荷载除施加岩体自重外，还将引入以“5·12”汶川地震震源机制解和已记录到的强震加速度为基础的人工合成地震波。两个方向的地震波加速度随时间变化曲线如图5所示，其x轴、y轴方向的峰值加速度分别为2 m/s2，-1.866 m/s2与1.608 m/s2，-2 m/s2。取计算步长为0.01 s。

岩质边坡稳定计算时，工况1仅考虑块体1的稳定性，而工况2则为分析块体1、块体2以及两者组合块体的稳定性。结果表明(见图6)：地震前5 s，两种工况条件下关键块体及组合块体的最小动力安全系数曲线基本重合。工况2条件下虚拟裂隙在地震作用5 s时，开始向真实裂隙扩展，滑动面力学参数随之发生变化，随着时间的推移，节理裂隙扩展进一步加剧直到17 s时，真实节理裂隙贯通整个虚拟节理面，关键滑移块体滑动面力学参数不再发生改变，17 s以后工况2条件下对应的最小安全系数曲线变化规律与工况1条件完全一致。不考虑节理裂隙扩展时，整个动力计算中的最大安全系数为3.912，最小安全系数为2.753。考虑裂隙扩展时，整个动力计算中的最大安全系数为3.395，最小安全系数为1.463。由此可见不考虑虚拟节理裂隙扩展会使得计算结果偏危险，而工程实际中反复外荷载的作用将使得岩体内部节理裂隙发生扩展，因此动力计算时考虑关键滑移面上的节理裂隙扩展十分必要。

图5 地震波加速度时程曲线

图6最小动安全系数时程曲线图

3 结 论

综合离散元法与刚体极限平衡法的优势，提出一离散元法及刚体极限平衡法(DEM-LEM)考虑裂隙扩展的岩质边坡稳定分析方法。研究表明，同一岩质边坡考虑裂隙扩展能够更加准确的模拟其内部破坏机理，对准确判断岩质边坡的稳定性十分必要。本次研究详细的给出了岩质边坡裂隙扩展平面问题的解决方法，对三维边坡裂隙扩展研究将是下一步研究重点。除此之外，本文仅考虑了岩质边坡滑移破坏模式，对于坠落、倾倒等破坏模式也是未来需要解决的问题。

参考文献：

[1] 杜应吉.地质力学模型试验的研究现状与发展趋势[J].水资源与水工程学报,1996,3(2):64-67.

[2] 戚靖骅,张振南,葛修润,等.无厚度三节点节理单元在裂纹扩展模拟中的应用[J].岩石力学与工程学报,2010,29(9):1799-1806.

[3] 王国庆.岩体水力劈裂试验及裂纹扩展的无单元法计算[D].南京：河海大学,2004.

[4] 寇晓东,周维垣.应用无单元法追踪裂纹扩展[J].岩石力学与工程学报,2000,19(1):18-23.

[5] 石路杨,余天堂.多裂纹扩展的扩展有限元法分析[J].岩土力学,2014，35(1):263-272.

[6] 张晓东.基于扩展有限元法及虚拟裂缝模型的混凝土断裂过程区分析[J].计算机辅助工程,2016,25(1):61-67.

[7] 蒋明镜,张 宁,申志福,等.含裂隙岩体单轴压缩裂纹扩展机制离散元分析[J].岩土力学,2015，36(11):3293-3300.

[8] Xiao H T, Yue Z Q. A three-dimensional displacement discontinuity method for crack problems in layered rocks[J]. International Journal of Rock Mechanics & Mining Sciences, 2011,48(3):412-420.

[9] 龚声武,陈才贤.急倾斜煤层放顶煤开采顶煤破碎的Griffith理论分析[J].湖南科技大学学报(自然科学版),2007,22(2):9-12.

[10] 邵 勇,王向东.基于断裂力学的高拱坝底缝稳定性分析[J].水电能源科学,2008,26(5):80-82.

[11] 罗 礼,孙宗颀.双扭法研究岩石亚临界裂纹扩展速度和断裂韧度[J].岩土工程学报,1992，14(3):40-48.

[12] 刘喜康,张建海,赵文光.基于动力有限元分析的任意滑块稳定分析方法[J].岩土工程学报,2013，35(S2):363-368.

[13] 冯 欣，熊子正，滕志强，等.基于积分法的均质边坡稳定性区间分析[J].水资源与水工程学报，2017,28(1):192-196.

[14] 陈革强.刚体极限平衡法浅析[J].海河水利,1999(2):16-18.

[15] 吴玉龙,张建海,李仁鸿,等.断层对孟底沟高拱坝坝肩变位及稳定性的影响研究[J].水利与建筑工程学报,2016,14(4):141-146.

[16] 雷远见,王水林.基于离散元的强度折减法分析岩质边坡稳定性[J].岩土力学,2006,27(10):1693-1698.