内部裂隙对节理岩质边坡失稳影响的数值分析

马鹏飞，郭德龙，许文年，陈 新*，夏 栋

(1.河海大学岩土力学与堤坝工程教育部重点实验室，江苏 南京 210024；2.中铁十四局集团有限公司，山东 济南 250014；3.三峡大学水泥基生态修复技术湖北省工程研究中心，湖北宜昌 443002；4.三峡大学水利与环境学院，湖北 宜昌 443002)

在形成现实岩体的过程中,自然界岩石综合体经受风化作用、冷热交替、水的破坏及改造,且在构造变动、工程建设、卸荷作用等条件状态下,促使了岩石内部各种宏细观裂纹裂隙的形成和发育[1]。这些损伤降低了岩层的强度并且不易被观察,是地下硐室群开挖、地基工程、岩质边坡等岩石工程的潜在威胁。地基工程和地下硐室工程基本上在有限变形理论领域里实施,其岩体或围岩变形量是受限的；而人工开挖边坡为一个有侧向临空面的半开放系统[2],在环境应力改变(如降雨、坡顶荷载)后既有瞬时变形又有与时间有关的流变属性,边坡一旦发生类似无限变形的显著位移即可造成塌方、滑坡和堵江等严重的地质灾害。因此,赋存于一定地质物理环境且作为工程对象的岩质边坡的稳定安全性问题,不仅自身属于自然科学问题,更是整个岩石力学界及工程地质学科最为关心的应用领域之一。

发育着不连续裂隙缺陷体的岩质边坡,在内外界因素诱发下,内部裂隙源将沿一定的路径发生拓展、汇合、互相衔接贯通直至最后宏观裂缝(或破碎带)的产生。工程实践及大量反复的试验已证实,岩体工程的失事大都与岩石节理的发育、原有结构面的衍生及其不稳定扩展演化过程相关联[1,3-4]。也正是因为实际岩石中含有具备不同特性且以不同组合方式组成的结构面,岩质边坡在失稳时,受裂隙的宏观结构特征及其力学性质的影响严重,造成其破坏模式繁杂多样,而节理面常被看作边坡关键性破坏面的组成部分。所以，对含内部裂隙岩质边坡滑移破坏机制与稳定性的进一步阐明具有重要而鲜明的指示意义。

由于岩体工程各种相互作用、材料性质及初始条件等问题的复杂性,以固体力学和经典数学为理论基础的理论分析无法解决现有的许多工程难题[5],故而要想给出服务于工程需求的结果，采用数值方法在不少情况下已成为一种不可缺少或者不得已的手段。岩土数值计算通过对工程实体的计算机模拟,能从整体上获得一些原本可能不了解的边坡变形失效细节,甚至在计算中也可以了解如何去分析边坡产生过度变形或破坏的原因,从而节约科研时间和经费。近些年新兴的基于断裂力学的动力扩展有限元法(XFEM)在处理物体内部不连续介质力学问题时特别有效[6-7],该方法不仅较为精准,并且还秉承了常规有限元法的优点,尤其在模拟剪切带不连续场、裂隙开裂以及裂纹扩展方面具有明显的优势。

回顾先前的研究,主要存在着以下不足：①针对边坡稳定性问题,岩土力学工作者采用强度折减有限元法基本对各种类型的边坡的稳定都进行了广泛的计算,但评价内部裂隙对边坡滑动面的影响及其演化规律的研究仍偏少,目前的边坡稳定性分析多是结合匀质边坡模型展开的,亦或将边坡模型建立为均质材料层的“堆砌体”,忽略了岩土体的裂隙性；②有些研究虽考虑了裂隙的影响,却直接将带裂隙服役的岩体通过设置很弱的强度参数来等效裂隙对边坡整体稳定性的不利作用,导致数值仿真结果失真[8-9],即无法模拟出预期理想的边坡介质材料开裂/破断效果。实质上,上述做法视边坡内部基岩或基土整体为连续体,未把裂隙等不连续界面夹杂其中,这也恰恰暴露了传统连续方法如有限单元法(FEM)、有限差分法(FDM)等的局限性。为此,本文以Abaqus程序软件环境下的XFEM二维模式为基础,结合温控参数的强度折减方法原理,通过几个算例阐介了将XFEM运用于裂隙岩质边坡相关问题研究的切实可行性，并简要探讨了不同的裂隙面排布状况及含量对岩质边坡的稳定系数、滑裂面演化路径以及滑坡最终破坏形态的控制,揭示工程岩体失稳的渐进破坏现象,以期为指导水利水电、公路、采矿工程中边坡的设计、稳定建设与科学修缮提供理论支撑。

1 基于扩展有限单元的强度折减法求解

1.1 动力扩展有限元法(XFEM)基本理论

有限单元法(FEM)在模拟微裂隙时,要预先确定裂隙扩展方向,并且划分的单元边界须遵从结构内部几何或物理界面,特别是在裂隙端部,计算网格的布置要求苛刻[6],形成了裂尖附近区域极其密集其余区域稀疏的不均匀网格分布,从计算机的运行成本方面考虑是不合适的。XFEM由力学家Belytschko和Black[10]于1999年首次提出,该方法采用水平集法描述裂隙位置,是在标准有限元框架内的重要延伸。随后,Mo⊇s等[11]进一步将阶跃函数作为结点富集项引入,弥补了传统有限元每个裂隙扩展步中都需要对裂尖区域重构网格以适应位移强、弱间断等不连续边界演化问题的缺点。

基于单元分解思想,引入反映裂隙局部特性的附加函数对有限元的逼近位移模式进行改进,将含裂隙物体的位移u(x)拆分为连续和非连续两项之和,得到XFEM的近似位移场如下[7,12-13]：

(1)

在此位移模式基础上,仿照常规有限元的虚功原理就可写出XFEM的支配方程(有时又叫作求解方程)[14-15],然后求解总刚矩阵、整体等效节点荷载矩阵,继而求解节点位移列向量。将XFEM应用于实际问题的分析中,裂隙的张开、闭合、生长和偏折行为完全与计算网格相互独立。

1. 2 温度相关参数强度折减理论及其方法实现

自Dawson等[16]提出强度折减法(SRM)以来,该方法就被普遍运用在工程稳定性系数(FoS，简称稳定系数)的自动搜寻中,其优点之一是不用事先知道破坏面。遗憾的是,Abaqus/CAE早年版本仅支持通过手动修改inp文件的方式来实现强度折减功能,不甚方便[17]。这里将等同于强度折减系数FV概念的温度Temp置于边坡体中(见图1),通过相对简洁的温度相关参数强度折减计算得出FoS。应当说明的是：①温度相关参数强度折减法是指由温度变量关联操纵整个强度折减过程的方法,其在使用时比较便捷,节省了在关键字中插入关于FV变化的描述语句这部分工作量,因此作者团队将该方法命名为“convenient strength reduction method”，简称C-SRM；②Abaqus新版本是否内置了强度折减法,笔者不敢断言；③下文各处出现的NT11就为图1中的Temp。

图1 抗剪断强度折减定义窗口(示例)

按量变到质变理论,由下面公式(2)对边坡岩体或土体所能提供的宏观抗剪强度黏聚力c和内摩擦角φ实施连续的线性折减[17]：

ck=c/NT11,tanφk=tanφ/NT11

(2)

式中:ck、φk分别为维持边坡平衡所需要或实际发挥的材料黏聚力(kPa)和内摩擦角(°);NT11表示坡体温度(量纲为1)。

为了防止边坡可能发生初始破坏,令NT11以固定步长从某一比1小的数值开始增加(温度初始值一般习惯设定为0.25、0.5或0.7,折减步长的选择则是依据对计算结果准确性的期望、对边坡安全性的预估及一些书籍上的通常做法等),将不同的ck、φk代入计算模型进行计算,判断边坡系统是否处于平衡状态，直至折减到边坡到达破坏或临界静力失稳状态为止,此时坡体的温度NT11即为平常人们所说的边坡安全系数(严格来讲应该叫作边坡稳定系数)[18]，且此时SRM自动识别出的最危险滑面所在的位置也即实际滑动面[19]。

C-SRM在ABAQUS有限元软件中的实施流程要点如下：在Property模块中定义与温度相关的c、φ；进入Load模块,先选择分析步initial,对集合rock创建温度分布场(宜事先对rock和crack-i部件各建立一个“Set”),温度值填入温度初始值,再选择分析步reduce,将温度进一步修正为某一足够大的终值。

强度折减数值计算中边坡失稳状态的判定标准主要有：①边坡域内最先形成连续的贯通区或塑性点呈大面积分布；②边坡特殊关键部位的位移量出现突增现象；③有限元应力-应变方程在有限迭代次数内不收敛。

在进行各类裂隙发育边坡的稳定性求解时,裂隙的存在诱致了所研究问题的极端非线性。数值计算结果不收敛情况有两种：第一种是在边坡达到极限平衡状态之前单元就畸变,并引发模型过程过早停止,事实上,因模拟方法不够完善这种现象较易出现,此时上述边坡失稳的判定依据①～③均无用武之地,所以务必让强度折减的整个过程持续到正常不收敛才可停止[措施包括：修改模型计算进程的默认分析参数和控制参数(分析步增量步长、允许的不连续增量数目I0、IR以及迭代不收敛次数IA等)；尝试重置网格单元边长]；另一种是在塑性应变区贯通后有限元分析不收敛,但边坡内特征点位置处此刻已产生了无穷大的形变(或位移-折减倍数曲线末端出现几近竖直的线段),说明边坡已或者早已破坏,故据此得到的边坡稳定系数偏大。因此，本文认为将边坡失稳状态的判定标准③作为裂隙岩质边坡失稳破坏标准是不可取、不明智的。

1. 3 屈服准则

岩土介质的非线性变形规律比钢材等其他材料要复杂得多,但边坡岩土体的屈服和破坏符合Mohr应力空间中的经典M-C准则,同时M-C准则因其形式和使用简单，在诸岩土失效准则中至今仍居于统领地位[20-22]，故文本在分析岩质边坡岩块体的强度问题时选用M-C条件,并认为材料屈服只取决于大、小主应力σ1、σ3,不考虑中间主应力σ2,其表达式为

(3)

理想弹塑性模型的库仑-莫尔屈服条件可以表示为[20]

(4)

式中：I1和J2分别为应力张量的第一不变量和第二偏应力不变量，这里I1的单位即应力的单位，J2的单位则是应力单位的平方;θσ为应力罗德角(-30°≤θσ≤30°)；F为屈服面函数。

1. 4 数值模拟的基本假定条件

本次数值模拟的基本假定条件如下：

(1) 边坡纵向视为足够长,在确保计算效率的前提下进行平面应变状态下的分析,对一些外轮廓形状沿纵向变化缓慢的边坡,多数情况下二维处理能近似满足工程要求[23-24]。

(2) 对节理岩体而言,其包含着节理裂隙面,在模型结构上类似于复合材料结构,节理岩体裂隙界面的力学性质弱于基岩部分的力学性质。为了降低研究难度,将裂隙简化为不占体积的空材料,并忽略其厚度,基质体按介质分布连续、均匀的理想弹塑性体考虑。

(3) 实际的工程岩体中存在规模尺度不一的大小裂隙,在现有技术条件下,很难把小微裂隙网络也纳入模型坡体,故暂只考虑内部长大裂隙,并将各种形态、特点的非连续结构面概化为坡顶裂缝和坡内斜直裂隙。

(4) 岩石材料的损伤起始判据应用最大主应力(Maxps)准则,其损伤演化规律则基于能量、线性软化、刚度衰减方面的规则。Maxps损伤准则认为裂纹/裂隙沿着最大周向拉应力达到最大值的方向发展,待满足界面裂纹起裂扩展准则时裂纹扩展。

(5) 滑体采用小位移计算模式,Abaqus-CAE前处理的Step模块中几何大变形开关仍拨置off,岩质边坡某处断裂的图形表示效果可通过变形放大倍数调控。

(6) 斜坡面为单一直线。

2 算例分析

相比同规模土质边坡,岩质边坡滑坡能量更大,毁坏性更强,所以岩质边坡稳定性评价成为广大同行共同关注的热点。本次数值模拟方案是建立多种人工设计节理岩质边坡算例背景并对算例背景中的边坡岩体在相互作用裂隙端部塑性区扩展影响下的稳定性响应(包括等效塑性应变、变形位移、折断过程等)逐一进行初探,以此研究XFEM在含岩体结构面边坡工程及在类似工程(特别如膨胀土边坡[25]、崩岗治理区崩壁土坡[26-27]等)中的应用。

2. 1 对比算例分析

2.1.1 含4条非共线裂隙(背景一)的岩质边坡算例

图2简示了某一假想的岩质边坡地质断面几何模型,设该边坡内部存在原始裂隙1～4(表面张裂隙1和埋藏型裂隙2、3、4)，数值分析模型中裂隙空间展布信息，见表1。

图2 含4条非共线裂隙(背景一)岩质边坡算例稳 定性计算的概化模型及裂隙布置(单位：m)

表1 含4条非共线裂隙岩质边坡算例裂隙模型几何参数

模型的单元库类型使用CPE4R(平面四节点四边形元,缩减积分),与完全积分单元相比,有意识地少取积分点是对高斯数值积分的一种处理,反而会使精度更高。计算网格大小以在模型分区边线上布结点数的方式设置,由于研究侧重于滑面附近及滑面左侧岩块的变形,故围绕含裂隙重点计算域hbcg的边缘上结点较稠密，以提高计算精度,而在扩展控制域外的gdeh区域分网较粗。模型有限元网格的结点数和单元数分别为1 395、1 317,但值得注意的是,不要对裂隙部件进行网格划分。地基的底边在X、Y方向上皆固支,在地基的左右两侧和岩质边坡的右侧安装辊轴支撑,坡面自由约束。

在建模分析时不需要输入任何关于内部结构面的力学参数,而输入的岩体(结构块)参数有：泊松比ν为0.29,弹性模量E为7.22 GPa,密度d为3.08 g/cm3,黏聚力c为100 kPa,内摩擦角φ为20°,塑性势角ψ为0°(忽视体积变形),能量释放率GC为500 J/m2,岩石所能提供的断裂能KIC为151.1 N/m,岩石损伤稳定黏度系数为1×10-5，重力加速度g按10 N/kg计。为了使含裂隙岩质边坡在自重一定的条件下产生变形破坏,借助材料降强的方式将边坡带入极限状态[21],得到如图3所示的该岩质边坡失稳演化中塑性区的变化情况。

图3 Abaqus塑性区云图显示的含4条非共线裂隙岩 质边坡算例滑裂面延展过程

由图3分析可知：①强度折减初期坡体的温度NT11(亦称降强系数)较低,也就是前期的边坡岩石抗剪强度指标值远高于实际黏聚力和内摩擦角,单元应力没达到屈服面,因此就未形成塑性区；②塑性区最先在裂隙3的上端部位被隐约发现[见图3(a)],紧接着是在裂隙3、4和2下端局部薄弱地带出现[见图3(b)]；③随强度折减程度的加强,裂隙扩展前端塑性区继续向前扩展,当强度折减到第15次时,裂隙尖端塑性区加速向对方延伸,即裂隙3、4之间塑性区明显接触,同时由于裂隙3靠近易遭破坏的坡脚,其下尖端与坡趾之间形成连通的滑弧通道，当强度折减到第17次时,裂隙4上端与裂隙2之间塑性区也互相连通；④伴随坡体温度NT11继续递增,裂隙2上部尖端塑性区逐渐向上发展,最终在后缘被张拉裂隙撕裂后,与裂缝1下部尖端相连接,坡内发育形成了从坡脚逐个连通裂缝,直至贯通坡顶的塑性区，此时的温度值约为1.916、强度折减次数为21,按照边坡工程失稳判据第①条,确定本算例边坡的稳定系数FoS为1.916。

图4 终极破坏时的含4条非共线裂隙岩质边坡 算例状态(两级破坏形态)

由图4可见,岩质边坡内含有多条裂隙的条件下,初始破坏面(指最先形成的贯穿滑面)形态并不呈圆弧状,而是相当不规则。分析其原因认为：岩体中的层面、节理、片理面、裂隙、裂缝、断层等地质界面的力学强度一般比岩石材料弱得多,这些优势结构面具控稳机制,体现在其对边坡的滑坡型式、变形规模和速率、位移场展布以及附近岩石中应力传递都具有控制作用[28]。待岩质边坡达到真正意义上的破坏状态时,更大规模的次级滑动趋势面已经酝酿完成。

边坡特征点处的位移值和坡体整体运动具有相同的趋势,为了结合位移突变理论得出边坡符合失稳判据②时的稳定系数FoS,图5中记录了该岩质边坡临空面端点c水平位移(U1)随降强系数(NT11)变化的发展趋势。

图5 非共线4条裂隙岩质边坡算例坡顶监测位点c的 水平位移(U1)随降强系数(NT11)的变化曲线

由图5可见：该岩质边坡临空面端点c的水平位移(U1)在坡体温度NT11由初始值增加到1.837的过程中只有微微的变化(当前阶段特点是没有滑坡迹象[29])；随岩石材料塑性参数c、φ值的持续下降,临空面端点c的水平位移U1产生了数量级的变化,位移切线角愈来愈大直到变形曲线上位移突变点显现后,-U1出现大幅度激增,表示滑坡过程由渐变阶段转化到了红色预警的急剧加速阶段[29],实际上此时边坡退化成了一个机构。

图5中位移拐点难以辨识(边坡稳定系数不易于合理地推求),但依旧可借助“双切线法”来辅助确定：过U1-NT11关系曲线末端及曲线上降强系数NT11为1.837的点分别做2条切线,交于点R,再过点R做铅直线,其与位移曲线的交点就是位移量拐点,对应的横坐标在2.083附近。所以在边坡失稳判据②下求得的边坡稳定系数FoS≈2.083,与采用边坡失稳判据①求得的FoS之间相差8.72%。考虑到边坡失稳判据①可操性强、直观性强、获得的FoS值偏于保守,且模型中裂隙之间应变突变区动态连通过程即为边坡渐进失稳的过程,故后文中统一依据塑性区分布情况来判别边坡的失稳状况。

2.1.2 含3条非共线裂隙(背景二)的岩质边坡算例

为了论证内部裂隙数量对岩质边坡稳定性与安全性的影响,对背景一岩质边坡算例模型稍作改变，即剔除裂隙4后重新计算。图6清晰地展示了6个典型强度折减分析时步时的模型等效塑性区发展云图。

图6 含3条非共线裂隙岩质边坡算例的失稳过程

由图6可以看出：含3条非共线裂隙岩质边坡算例的塑性区首先在裂隙3的两尖端处形成,然后裂隙2端部也出现了不明显的塑性区；随着强度折减过程的推进,边坡坡脚处出现了局部剪坏面,并与裂隙3下端连通,且按弧线状逐渐向上开展，但是由于受裂隙的影响迫使失稳滑面沿两个方向扩展,一条是与裂隙2尖端相连接的主滑面(一级滑裂面),另一条为按原路线延伸的二级滑裂面[图6(d)、6(e)]；在第20次强度折减完成时,边坡域内塑性应变已贯通,但该岩质边坡还未进入整体破坏,且在此后边坡塑性应变继续增大、二级滑裂面继续拓展；当塑性区发展到一定程度(第450次强度折减),有限元计算不收敛,该岩质边坡也彻底破坏,次级破坏面发展为与裂隙2交汇的第2条贯通塑性面[见图6(f)和图7]。

图7 等比例放大1 717倍的含3条非共线裂隙岩质 边坡算例网格变形图

总体言之,若裂隙之间距离较远,则两裂隙中间部分可视作无裂隙扩展，换言之，随着坡体温度的升高(在现实中表现为作用于边坡的自然力及其作用时间的增加,或人类应力扰动的日益加强),域内已形成的广义塑性应变区向外传播,直至塑性面延展到与边坡内某个与之最靠近的预制裂隙(“nearest-fracture,简记N-F”)之间距离足够小时,塑性面之后的扩展才会显然受到N-F的影响(对应于本文,当强度折减到图6(c)的状态时从坡脚引出的滑面才开始分叉)；相反,若内部裂隙之间距离较近,两裂隙端部塑性应变则易发生连通(如裂缝1、2相隔较近,两者之间塑性区先于裂隙2与3间的塑性区连通)。

2.1.3 含2条非共线裂隙(背景三)的岩质边坡算例

设上述岩质边坡内仅有裂隙1和2,其余计算条件均不变,重新计算。为了研究分析，图8展示了施加加速度g引起模型在第13、17、20、48次强度折减时的积分等效塑性应变云图。

由图8可以看出：该岩质边坡破坏区在裂隙两侧端面最早出现,这一点在前述诸背景算例中也得到体现，这是由于事实上裂纹/裂隙是应力释放和应力产生奇异的关键部位,其局部高应力必会致使材料屈服,直接影响着岩体的变形和稳定，此外,众所周知脆性岩石抗压不抗拉,且边坡滑动时岩石的抗拉强度未能充分发挥,故塑性区首先出现在围绕裂隙尖端一个小范围拉应力集中区域；而后,坡趾附近受力较大也渐渐被征服,当强度折减到第20次时求得了最弱的滑动面[见图8(c)],它是刚刚贯通的临界破坏区,此时边坡转动破坏机制为坡面破坏；主滑移面贯通后如强度折减过程继续下去,第二条塑性面也将贯通,直至边坡发生最终的坡脚破坏[见图8(d)]。

图9给出的是该岩质边坡滑动位移矢量图。

图9 含2条非共线裂隙岩质边坡算例的最终破坏 形态(变形放大595倍)

由图9可大概判断出该岩质边坡稳定失效时的滑体范围，但因裂隙2端部距离坡脚的位置太远,或者说裂隙2边缘与坡脚之间的部分相当于无裂隙区,坡内预制裂隙几乎影响不到二次贯通塑性区的发展路径,所形成的过坡趾的整体性潜在失稳滑面(与之对应的降强系数为2.331)与裂隙1、2都不相交。

2.1.4 不含裂隙(背景四)的岩质边坡算例

在背景一岩质边坡算例的基础上,删除计算模型中的所有裂隙。因均质岩质边坡不含裂隙,滑动面剪入口自然就不会承受后缘裂隙的影响,为了降低边界效应,根据圣维南原理(仅仅对周围一定区域有影响,对更远的地方影响可忽略不计),现将图2中的cd边、of边往右加长了437.5 cm。不难看出,边界范围已给图8(d)示出的计算结果带来了一定影响。

采用传统有限元法对此模型的稳定性进行分析,利用最后2个强度折减步(Frame)得到该均质岩质边坡增量位移分布云图[见图10(a)],可判定其潜在滑带大致呈弧形,符合完整边坡滑动面的形态特点,这与考虑裂隙发育岩质边坡的破坏机制存在较大的差异。温控强度折减过程中该均质岩质边坡最早明显贯通的破坏面，见图10(b)。

图10 均质岩质边坡内部大致的滑动面形状及位置

2.1.5 含2条共线非贯通裂隙(背景五)的岩质边坡算例

图11 含2条共线非贯通裂隙的岩质边坡示意图

设定的岩质边坡岩石物理力学性质同第2.1.1节,图12大致展示了模型的运行结果。

图12 利用XFEM模拟的含2条共线非贯通裂隙(背景 五)岩质边坡连续塑性破坏面和裂缝的演化

由图12(b)可以看到,计算域内出现一条通过非连续弱面且贯穿的潜在滑面,在裂隙1′起裂以后边坡将启动平面滑动型失稳模式。

结合软件中输出参数PHILSM(裂纹面水平集函数)，可查看计算不收敛时裂隙扩展的结果,结合响应量STATUSXFEM(裂隙扩展状态参量)可得两裂隙之间的损伤值(1代表该单元完全开裂；0代表该单元无裂隙),见图12(c)和12(d)。即,在边坡系统到达最后失稳破坏状态时,裂隙5、6在各自两尖端处沿着与岩桥呈某一钝角方向有所扩展(长度约6～12 cm),但萌生出的次生或翼裂隙较为短小。

综上,表2列出了不同背景算例下计算得到的岩质边坡最小稳定系数FoS。

表2 不同背景算例下计算得到的岩质边坡最小稳定系数FoS对比

由表2可知：一方面，只要岩质边坡内存在裂隙,含裂隙岩质边坡的FoS就明显低于均质岩坡的FoS；另一方面,当裂隙数目一样时,即使裂隙6比裂隙3长了1 m,但背景二下岩质边坡的FoS比背景五下岩质边坡的FoS还小,这可能的原因是裂隙尖端位置和两条平行埋藏裂隙的倾角不同。

假设有这样一种背景算例(命名为**)：将背景二中裂隙2和裂隙3绕各自的下部尖端逆时针旋至与bc面平行,则旋转后的裂隙2、3也接近共线,其他计算条件均不改变。可计算出背景**下岩质边坡的FoS为2.290(>2.086)。可见,节理边坡稳定系数对滑动层面(结构面)倾角的变化比较灵敏,因此在滑坡治理工程中,针对内部含不同倾角、产状裂隙的岩质边坡应给予区别对待。

该岩质边坡(指图11模型)的最终破坏形式及相应的总位移分布特征，见图13。

图13 含2条共线非贯通裂隙(背景五)岩质边坡的 最终破坏形式及相应的总位移分布特征

由图13可见，沿着被该岩质边坡裂隙切割出的复合滑动面上方单元岩块位移明显大于其下方滑床位移,滑裂面的上部和中部处均有脱离母岩现象,该模拟研究结果与预期变形情况(即岩质边坡已滑出失稳并向坡底滑落)的符合度较高。

2.2 独立算例(滑移路径呈阶梯状的岩质边坡)分析

在中国西部山区,由于岩体结构经历了多期表生演化作用,阶梯式贯通变形破坏常见于受浅表部陡、缓两组结构面控制的卸荷岩质高边坡的破坏模式中,其破裂面形态整体上为陡-缓相间的阶梯状[30]。为了进一步考察XFEM的可用性与推广性,本文以文献[31]中的平行断续节理岩质边坡模型为实例(与第2.1节无关),就更为复杂难解的裂隙性边坡变形稳定问题尝试开展额外的仿真和分析。该岩质边坡剖面地形结构如图14所示,内含有一组与边坡倾向一致的平行等间距主控节理裂隙,岩石主要物理力学参数的取值见表3。

图14 平行断续节理岩质边坡模型(90°岩桥 )[31]

表3 岩石主要物理力学特性

对图14模型进行扩展有限元计算,得到该岩质边坡在程序正常不收敛时的裂隙扩展情况、运动矢量及相关等值线图(涉及塑性应变、剪应力S12、水平位移U1、最大主应力),分别列于图15。经数次反复的温控强度折减试算可知：当塑性区刚刚从坡脚到坡顶贯通时对应的岩质边坡模型温度(即FoS)约为5.255；如果清除该岩质边坡内部的节理裂隙,由FEM求出该岩质边坡的FoS约为8.17。

图15 独立算例的动力扩展有限元模拟结果汇总

由图15可得出几点有参考价值的信息：①观察图15(a),在原始裂隙尖端处,裂隙发生了扩展,同时伴有翼裂隙,岩桥失效,断续结构面间在宏观层面形成连通；②坡体整体失稳后的塑性应变带由4个相互连通的近椭圆状的塑性区域组成,见图15(b),并且椭圆长轴的两端点基本位于原始主控裂隙面的中心,椭圆短轴的两端与原始裂隙尖端重合；③播放图15(d)的动画发现,坡体先沿临近坡脚的节理产生局部滑移,而后再沿着由断续节理系以及扩展萌生的裂隙所组成的滑动面发生整体朝下的阶梯状贯通滑裂；④以贯通面为界线,边坡位移场分布在模型内产生明显分异(与图13相一致)，且对于图15(e),X轴方向位移值在可滑动块体的脚部位置处达到最大,坡体中部位移较大,坡肩部位位移次之,这说明阶梯状变形破坏机制具有“牵引渐进式”破坏性质；⑤如图15(f)、(c)所示,节理裂隙的存在将引起岩石应力重分布，裂隙尖端部位应力较为集中,说明开裂扩展还可能会继续下去，而相对于边坡其他地带,节理周围处剪应力大小也较大,这与周子涵等[32]利用FLAC3D得到的数值计算结果很类似。

3 动力扩展有限元法(XFEM)的优劣分析

以上建立的各背景算例岩质边坡数值模型大体上达到了预期效果,有限元计算所得的一系列成果图与现有文献中的研究结果可良好吻合(如边坡主滑动区集中出现在岩体结构面上部坡体区域)。基于XFEM的温度相关参数强度折减法(C-SRM)意义明确,不仅可以一次性地让Abaqus算出断续裂隙岩质边坡的滑裂面位置及最小稳定系数值,而且有限元计算中无需裂隙力学参数的帮助就可正确地反映岩质边坡的变形及滑动面孕育过程(若在以往,要通过假定、室内测试、经验公式等方式事先获得节理裂隙的强度参数),同时本文提出的方法能针对折线滑面以及同时包含折线和曲线的滑面进行分析,并能较好地描述裂隙对边坡稳定状态的削弱作用,间接验证了基于XFEM的C-SRM具有更加高效、合理与简洁等优越性。

但是，本文方法也存在一些局限性：无法预测内部含完全贯通节理(包括节理从坡面延伸到坡顶和从坡趾延伸到坡顶两种情况)岩质边坡的稳定系数；无法对边坡坡内的弯折裂隙及交叉型裂隙进行模拟；此外,XFEM计算的收敛性尚需提高。

4 结 论

不是使用过去习惯的传统有限元强度折减法,而是在提出基于温度场变量的便捷降强计算方法的基础上,本文改为采用XFEM(从而使内部含非连续弱面这一类型边坡的稳定问题能够被有效地处理和评价)研讨和开拓了滑坡全过程中多裂隙岩质边坡内部裂尖塑性区的萌生、发展规律,以及几种简单的节理岩质边坡模型在静力荷载下可能的失稳模式。通过研究获得了如下结论：

(1) Abaqus中的XFEM提供了直观的裂隙形式,有别于条分法等经典的极限平衡法,采用XFEM得到的岩质边坡主滑面更符合滑坡形成特点,由于日后可考虑某具体边坡工程后缘张裂隙的存在,故不会高估岩质或土质边坡的稳定系数；而且将XFEM与强度折减法结合还能寻找出两个潜在的失稳滑面。因此,本文所做工作为XFEM在裂隙边坡稳定性计算领域的应用提供了有益参考，建议将XFEM在相近工程实践中加以推广应用。例如,可考虑利用XFEM模拟南方7省区崩岗区花岗岩残积土裂隙的动态行为,这无疑对崩壁崩塌机制的深层次和多角度发掘具有重要的意义。

(2) 裂隙的存在缩短塑性区贯通的路径长度(就好比裂隙扩展必须要沿着指定的薄弱面),促使边坡的抗力下降,因此一般情况下,裂隙数目越多,边坡的稳定系数就越小。

(3) 完整岩石坡体在外部荷载驱动下的失稳破坏面是穿过坡顶的典型圆弧形或对数螺线形滑裂面；含裂隙节理岩质边坡的破坏始于裂隙尖端形成的塑性区,其失稳破坏面是通过断续裂隙的不规则滑移面,且滑面形态与裂隙的发育程度及排布息息相关；阶梯状滑移方式形成的连通塑性区是由若干离心率较小(丰满)的椭圆塑性带与原有节理面共同组合而成的。

(4) 受内部裂隙的影响,岩体边坡处于极限状态时的位移场特征形成了以贯通滑面为界线的明显分区；在顺着裂隙长度方向上,边坡岩石应力分布规律也呈现出很强的统一性。

限于篇幅,下一步的研究工作拟开展含连续岩体节理(完全贯通)、内部裂隙倾角变化范围更广(如优势结构面逆向坡、反倾层状边坡)、坡顶存在多条竖直裂隙等工况作用下岩质边坡稳定性及其滑移规律的数值模拟研究,并与已有解析解(含有解析解的直线滑面的边坡算例)做对照分析。

致谢:谨向参加本论文审稿并提出建设性修改建议的专家、老师致以最真诚的谢意！