基于局部强度折减NMM的岩质边坡稳定性分析

2019-05-13李志军朱承金周海霞王述红王艳飞胡署光

水利与建筑工程学报 2019年2期

李志军，朱承金，周海霞，王述红，王艳飞，胡署光

(1.中交隧道工程局有限公司，北京 100020; 2.东北大学资源与土木工程学院，辽宁沈阳 110819)

全世界范围内均有大型岩质边坡失稳灾害频发，阻碍了社会发展，同时也对人类生命和财产安全造成了不可挽回的损失。准确地对灾害发生位置进行预警预报，并及时对危险部位进行支护加固处理，保卫人类生存空间，一直以来都是工程界亟待解决的难题[1-4]。

强度折减法近年来在边坡稳定性分析中应用广泛。郑颖人等[5]发展了强度折减法对边坡安全系数的求解。后又衍生出局部强度折减法，其凭借独特的对局部岩体进行折减且不必事先确定失稳面的优势被广泛应用及推广在岩质边坡安全性分析中，该方法注重岩体本构关系、开挖及支护组件的相互影响，相较于常用的极限平衡法具有绝对优势，但其理论原理限制了岩体破裂过程的模拟实现。由于岩体破裂过程及其复杂，现有的主流数值分析软件，例如有限元在对裂缝扩展问题分析时，需要不断的更新网格，使得前处理工作复杂。由Shi[6]首创的数值流形方法，通过数学覆盖定义权函数并控制计算精度，通过物理覆盖定义系统自由度；NMM继承了与DDA相同的接触检索计算方法，能够很好地模拟非连续面之间的挤压、摩擦等相互作用以及块体的平移、转动等大变形特征，并在统一计算架构下灵活而有效地处理连续非连续问题。而恰恰岩质边坡失稳是一个连续-非连续的渐进破坏演化过程，其连续-非连续特质决定其可以较好地揭示岩质边坡失稳的渐进破坏过程，实现失稳过程的动态表征。Zheng等[7]运用无网格逼近函数，实现了无网格流形法的构建。Ning等[8]基于Mohr-Coulomb准则组建NMM的裂纹扩展准则，并利用改进的NMM方法对边坡失稳现象进行了计算分析。Wong等[9]采用NMM模型对岩石边坡的渐进破裂过程实现了高精度模拟与计算。

岩质边坡稳定性评价核心问题在于安全系数及相应失稳面的求解[10]。安全系数求解后，由于传统局部强度折减法无法很好地进行失稳面滑落模拟，对其位置确定的精确性及后续失稳的持续性还存在较大疑问，无法确定该安全系数是否为真实安全系数。本文为进一步完善岩质边坡安全系数及相应滑落面的确定手段，借鉴NMM模拟大变形、大位移的优势，考虑边坡实际的局部损伤性，提出LSR-NMM，并将其应用到重庆兴隆隧道进洞口边坡稳定性评价中，且与传统强度折减法、现场监测数据进行对比分析，验证方法的可行性及计算精度，并给出滑落面位置与相应安全系数，为现场下一步治理提供参考性意见。

1 LSR-NMM原理

1.1 NMM有限覆盖系统

NMM创新性地采用了数学覆盖(Mathematical Cover, MC)和物理网格(Physical Mesh, PM)两套覆盖系统，能够在同一计算架构中进行连续和非连续问题的处理和计算[11]。MC系统由人为进行定义，可以是规则网格或是多个任意形状分片的集合，并且该集合要求能够覆盖整个求解域。PC系统则是由物理网格(材质裂隙、边界等)对MC系统的切割形成的面域。图1给出了NMM覆盖系统与流形单元的生成过程。如图1所示，MC由三个数学覆盖分片(矩形分片M1，圆形分片M2和M3)重叠组合而成，经过物理网格(物理边界-α和内部非连续-β)的再剖分，形成了包含7个PCs(P11，P12，P13，P21，P22，P31和P32)的物理覆盖系统。流形单元(Manifold Element, ME)由PCs的重叠区域或者单独的PC构成。如图中所示，单元E3由P12和P22两个PCs的重叠区域构成，单元E7由P13，P22和P32三个PCs的重叠区域构成，而部分单元则由独立的PC构成，如E1由覆盖P12独立构成等。

图1有限覆盖系统与流形单元

在物理覆盖系统形成之后，NMM在每个PC上构建局部位移近似函数，通常采用多项式进行局部位移近似：

ukpc(x)=PTd

(1)

式中：x为位置矢量；PT为多项式的基函数矩阵；d为位移自由度列阵；k为对应PC的编号。在原始NMM程序中，局部位移近似函数一般采用0阶线性位移函数进行近似：

(2)

式中：dk1，dk2为位移函数的系数。通过权函数对式(1)进行加权平均，即得全部PC系统上的整体位移近似为[12]：

(3)

构成ME的PC数量在式(3)中用npc表示，综上得全局位移函数表达式：

(4)

Nk(x)为形函数，wk(x)为覆盖k对应的权函数，且满足：

(5)

1.2 系统控制方程

NMM基于最小位能原理构建系统的总体平衡方程[13]。通过对各项的计算获取系统的总位能：

(6)

基于平衡条件：

∂Π=0

(7)

推导得到系统的控制方程为：

KD=F

(8)

式中：K表示系统的整体；D和F分别为系统的位移和荷载。NMM中系统的控制方程为线性方程组，总体刚度矩阵为带状稀疏矩阵，原NMM计算程序中采用逐次超松弛迭代法(Successive Over-Relaxation, SOR)进行方程组的求解[14]。

1.3 LSR-NMM实现原理

LSR-NMM在全单元强度折减的基础上，获取剪切带位置，提取剪切带岩体单元，仅将剪切带位置的部分岩体作为计算模型，折减其黏聚力c及内摩擦角φ，持续减小该部分岩体参数，直至计算不收敛，此时岩质边坡进入失稳临界状态，安全系数F可由式(9)表示。

(9)

由公式(9)可知，F随τa的变化而改变。利用极值方法对函数F进行最小值的求解，结果如式(10)，并将其定义为岩体单元安全系数Fs，Ds(σs，τs)为相应在应力圆上的点，其坐标见式(11)[15]：

(10)

(11)

过Ds点作摩尔圆切线，使得切线方程为

τ=tanφ′·σ+c′

(12)

Rc和Rφ各自表示黏聚力与内摩擦角的折减系数，Rc与Rφ可表述如下：

(13)

联立式(9)—式(13)可导出：

(14)

由式(14)可知，岩体单元安全系数与使得岩体单元处于临界位置时的折减系数相一致，进而安全系数就可以通过LSR-NMM来求解，使岩体达到临界位置，进而求得岩体安全系数。

2 工程应用

2.1 重庆兴隆隧道进洞口边坡概况

兴隆隧道位于重庆市渝北区，进洞口边坡坡向295°～298°，坡角60°～80°大部分基岩裸露，主要不良地质现象为崩塌问题，边坡全貌如图2所示。

图2兴隆隧道进洞口边坡

据工程地质勘察报告，仰坡地层岩性分为3层，自上而下为强风化泥岩、中风化泥岩和强风化砂岩。连续施工及持续降雨的原因导致坡面已出现多条张拉裂缝，该边坡如果继续受到工程扰动或其他非可控因素影响，根据经验可初步断定其存在发生失稳可能。布置深部位移测斜孔，为准确判断出岩质边坡整体失稳的潜在滑面和安全系数，向下一步治理工作提供治理方案及建设性意见。

2.2 现场监测

测点间设置间隔是0.5 m，近两个月内各测斜孔不同深度测得的累计位移值见图3。

图3各测点累积总位移值

据图3分析可得,P1孔总位移值较小且随时间增长较慢，基本稳定在1 mm～3 mm；P2孔总位移值随时间增长较大且存在多处位移值突变，最为明显的一次突变发生在深度约为23 m，瞬时突变值约为1.5 mm～2.5 mm，23 m之上总位移值较大，维持在4 mm～9 mm，23 m之下总位移值较大，绝大部分维持在1.5 mm～4.0 mm，说明P2孔深度为23 m附近岩体存在较大失稳嫌疑；P3孔总位移值随时间增长较大且存在位移值突变，最为明显的一次突变发生在深度约为13 m，瞬时突变值约为1.5 mm～2.0 mm，说明P3孔深度为13 m附近岩体失稳可能性较大。

2.3 兴隆隧道进洞口边坡稳定性分析

对工程现场所采集岩体进行室内试验可得其岩体参数，具体数值如表1所示。

表1 兴隆隧道进洞口边坡岩体参数

基于LSR-NMM对边坡失稳过程进行动态模拟，结果如图4所示。

图4边坡失稳过程(单位：mm)

首先对全单元进行折减，获取剪切带位置，提取剪切带岩体单元，仅将剪切带位置的部分岩体作为折减计算对象，持续折减其黏聚力及内摩擦角，直至计算不收敛，此时岩质边坡进入失稳临界状态，计算得出1.176为该边坡安全系数,且滑落面与现场位移监测预测滑落面基本一致。利用LSR-NMM模拟动态失稳过程的先进性对此工程岩质边坡下一阶段变形破坏进行了预测，在判断出安全系数与滑落面实际位置的基础上，为减小边坡失稳造成的损失，根据模拟结果提出治理方案，依照模拟结果，对上部处于强风化泥岩区的预失稳岩体不同位置进行剥离，其已经不能为岩体稳定提供支持，故剥离之，减小即将失稳岩体体积，剥离后，对其下部分进行注浆处理，将其固定在中风化泥岩区，保证边坡稳定性。