基于颗粒流的顺层岩质边坡稳定性分析

李丽君,龚 盛

(1. 新疆塔里木河流域干流管理局,新疆 库尔勒 841000; 2. 西安科技大学 建筑与土木工程学院,西安 710054)

0 前 言

边坡工程在水电工程建设中的地位较为特殊,其稳定性直接关系到工程安全和生产安全[1],水电工程中的边坡失稳不仅会影响电站的正常运行,也威胁着下游人民的生命财产安全。边坡稳定性计算中常用的方法有极限平衡分析法和数值分析法。随着计算机技术的不断发展,数值分析方法被广泛地运用于边坡稳定性分析当中,采用较多是有限单元法(FEM) 、有限差分法 ( FDM) 和离散元法(DEM),这3种方法中有限单元法和有限差分法均假设工程对象为连续介质,忽视了岩土体颗粒间的非连续性、非均质性和高度非线性,在计算上存在一些不足。

离散元法是由Cundal[2]基于动力学方程提出的一种显式求解算法,显式求解是对时间进行差分,在求解过程中,每个增量步内不需要进行迭代求解,无需形成切线刚度矩阵,对计算机的内存需求较小,每个增量步内计算量相对较小,但若求解的时间过长或计算步长过小均会导致计算时长大幅度增加。颗粒流方法( particle flow code,PFC) 是一种特殊的离散单元法,利用质量刚性体来模拟岩土体颗粒,适用于研究非连续介质力学问题,颗粒之间可彼此分离,并在计算过程中自动识别新接触,能够反映岩块或土体之间的滑动、平移和旋转等,在大变形计算方面有较大的优势[3]。

Wei[4]等基于数字高程模型构建了四川乐山马边市滑坡的颗粒流模型,从滑坡速度、位移和能量特征对边坡破坏过程进行了分析,模拟结果与真实的滑坡特征和滑坡运动过程非常吻合;王培涛[5]等通过定义安全系数和选择合适的失稳判断准则,将强度折减法引入离散元计算,实现了边坡安全系数的定量化判定。戴建建[6]等以某边坡工程为例,通过绘制特征颗粒位移与折减系数关系曲线求解边坡稳定系数,得到的边坡安全系数与有限元强度折减法和Bishop法结果一致;代远[7]等基于数值双轴试验标定细观参数,建立颗粒流边坡模型,采用局部强度折减法对滑带内土体强度进行折减,计算结果表明局部强度折减法所得边坡破坏面更接近极限平衡法;Tao[8]等利用颗粒流方法研究了提防中的管涌侵蚀破坏问题,较为系统地揭示了渗流方向上局部孔隙水压力和局部接触力分布规律;张帆[9]等通过刚体极限平衡法与离散元法的计算分析,得到了较为相近的边坡安全性系数,且后续监测反映的现场情况与计算较为相似。唐红梅[10]针对三峡龚家方2号斜坡采用PFC2D进行建模计算,将边坡滑动过程分为4个阶段,得到坡脚出现“底鼓”,其为边坡是否滑动的重要预兆,得出模拟结果与实际破坏情况相一致的结论。张志飞[11]等基于 PFC软件进行计算,并在考虑岩体结构面参数的基础上对反倾层状边坡进行了模拟,数值模拟结果表明岩层倾角对边坡变形破坏模式有很大影响,随边坡倾角增大边坡破坏方式由滑移型转为倾倒破坏型。

本文基于颗粒流方法,以西北地区某水电工程岩质边坡为例,建立颗粒流计算模型,通过双轴试验标定岩体细观力学参数,并运用强度折减法进行边坡稳定性计算,分析边坡变形失稳特征。

1 工程概况及岩体参数的确定

1.1 工程概况

某水电站位于青海省区域,库区边坡距离峡谷出口约10 km,距离梯度上游约50 km。水库正常蓄水位2 990.00 m,坝高254 m。边坡整体主要以破碎的板岩和砂岩为主,以散体结构为主;滑体中的板岩相对完整,边坡地质体整体呈层状～碎裂结构。

1.2 边坡岩体细观参数标定

基于颗粒流方法计算所需参数称为细观参数,细观参数与宏观参数存在相关性,但并非一一对应的关系[12],因此需要通过数值试验来进行标定。边坡相关岩体宏观力学参数通过现场试验测定,参数如表1。颗粒流计算细观参数主要指颗粒及黏结间的参数,颗粒的细观参数包括颗粒法向刚度、切向刚度和摩擦系数等;黏结间的细观参数包括法向黏结强度、切向黏结强度和黏结模量等。通过颗粒及其黏结间的细观参数可以反映岩土体的宏观特征。

表1 边坡宏观岩体力学参数

本文利用双轴模型试验来标定颗粒的细观参数,为贴近岩土体的真实物理特性,颗粒的半径越小越好,但由于计算机容量与计算速度的限制,可适当放大颗粒半径,双轴试验中所用粒径范围为0.1～0.15 m,,共生成颗粒1336个,同时需考虑双轴试验中尺寸效应的影响,即模型短边需大于40倍的平均粒径,因此数值试验选取的试样尺寸为6 m×12 m(宽×高)。双轴压缩试验模型如图1所示。

图1 双轴压缩试验模型

首先在计算域内根据颗粒大小和粒径比生成颗粒;然后颗粒赋予初始接触模型,由于没有施加黏结,所以是线性模型;然后缩小颗粒位置,使得边界颗粒与墙体有一定距离,可防止后续在计算平衡过程中颗粒逃离至墙体外,添加密度和局部阻尼,进行初始平衡计算;将初始平衡模型添加平行黏结模型,消除线性力,改为增量模式,清除颗粒位移和角速度;最后采取位移加载的方式,上侧的边界向下移动,下侧边界向上移动,直至试样失稳破坏,并在加载过程记录加载应力、应变。

图2为试样破坏时的黏结图。加载初期试样局部黏结首先产生张拉破坏裂纹,随着加载的进行,张拉破坏裂纹逐渐增多,并且开始出现剪切裂纹,当试样达到峰值强度后,岩样内产生的微裂纹数量迅速增加并出现宏观裂纹。这一过程与室内试验的岩石破坏形式较为一致[14]。

图2 颗粒黏结状态

为获得能表征岩石宏观物理力学特性的细观参数,采用不同围压对双轴试样进行多次加载,并绘制不同围压下的摩尔应力圆,画出所对应的强度包络线,并通过多次试算与宏观参数进行对比,表2为数值试验试算对比得到的力学参数值。

表2 颗粒流模拟的细观参数

以板岩为例,采用表2的细观参数,对板岩进行围压为100、200 kPa和300 kPa的双轴试验,可得试样峰值强度分别为431、625 kPa和823 kPa,绘制摩尔应力圆及摩尔库伦包络线如图3所示。计算包络线斜率及截距可得黏聚力为93 kPa,内摩擦角为17.96°,对比表1认为,可在后续边坡稳定性分析中运用该组细观参数。

图3 摩尔应力圆与摩尔库伦包络线

2 颗粒流模型构建

本文采用膨胀法建立边坡模型,并在边坡的四个部位布置了8个监测点,所建立的模型如图4所示,具体建模步骤如下

图4 边坡离散元模型

(1) 在CAD中画出模型图,以DXF格式将边坡模型导出,并记录下各个地层的面积。将DXF文件导入软件,建立模型的计算边界墙,并利用记录的区域面积,并通过颗粒平均半径和孔隙率计算出各个区域所对应的颗粒数量。

(2) 在计算域内生成小颗粒,然后对颗粒最小、最大粒径按等比放大,颗粒粒径与双轴试验保持一致,进行应力平衡计算直到充满模型区域。

(3) 初始试样有些颗粒处于悬浮状态,这些颗粒与周围颗粒的接触少,可以认为对力学模型不起作用,但会导致模型不收敛,直接删除悬浮颗粒会导致局部孔隙率变大[15],为了不删除悬浮颗粒,每次将悬浮颗粒的半径增加1.05倍后,再次遍历悬浮颗粒,当有悬浮颗粒时,再次增加半径。

(4) 对填充好的模型进行初始化应力场和消除漂浮颗粒运算后赋予模型相应的细观参数。

3 岩质边坡稳定性分析

3.1 基于强度折减法的边坡失稳分析

边坡稳定性计算主要有两种方式:一是通过降低岩土体强度参数达到极限平衡的强度折减法;另一种是通过增大岩土体重度或重力加速度实现超载以达到极限平衡的重力增加法[16]。

在有限元及有限差分计算中,可以直接折减岩体宏观强度参数即摩擦角和黏聚力,将边坡达到临界平衡状态时的折减系数作为稳定性系数;颗粒流计算时采用的是细观参数,颗粒流法在折减时主要是针对黏结强度与摩擦系数进行折减,折减时同样采用等比例折减方法,寻找坡体失稳破坏的临界值[16]。因此,颗粒流强度折减法稳定性系数可以定义为:折减前的细观参数平行黏结摩擦角(pb_fa)、平行黏结内聚力(pb_coh)和摩擦系数(fric)与折减后的临界强度参数的比值,即式(1)所示:

(1)

在初始边坡应力平衡计算完成后的边坡模型强度折减计算过程中,首先对初始模型计算结果中的 位移场和速度场进行清零;然后采用二分法通过对砂岩的岩体力学强度参数,即黏结强度与摩擦系数分别除以折减系数;最后对赋予新的强度参数的模型进行等同于初始场计算的运算步;当达到了预设的精度时停止计算,最大不平衡力达最大接触力的1×10-5为预设停止计算的精度。通过上述步骤以模拟边坡岩体力学参数劣化过程中所产生的变化,为揭示边坡潜在问题、控制因素及影响范围提供参考。

3.2 边坡变形失稳特征分析

根据式(1)进行滑坡强度折减模拟。图5为折减过程中边坡关键监测点位移随折减系数的变化图。由图5可知,随折减系数的增大监测点位移逐渐增大,当折减系数超过1.23时边坡位移突增,因此认为当折减系数为1.23 时,滑坡刚好处于临界位移状态,此时的折减系数即为滑坡的稳定性系数。

图5 监测点位移随折减系数变化

边坡破坏过程模拟如图6所示,显然,坡体的位移以砂岩层滑体为主,根据计算所得不同时步的边坡位移分析认为边坡的破坏模式表现为典型的牵引式滑坡特点。在折减到1.23时颗粒间黏结产生大量破坏,使砂岩层底的颗粒向外滑出,导致上部颗粒失去支撑而向下滑动。由于坡体前缘发生滑坡,因而失去支撑面而向后缘牵引,使得坡表发生整体滑动。但因为坡体下层板岩的强度较大,折减至1.23时不足以使得下层板岩间黏结发生破坏,因此边坡会出现砂岩层整体向下滑动出坡面。

图6 边坡破坏过程位移

模型运行至5000时步时,边坡坡脚在重力的作用下产生了0.1 m左右的位移,此时边坡还未形成明显的剪切带;运行至10 000时步时,边坡位移依然集中在坡脚部位,其位移增长至0.14 m左右,已在边坡下部初步形成剪切带;运行至30 000时步时,边坡最大位移增长至0.45 m,剪切带沿砂岩与板岩的分界层逐渐向边坡上部发展;运行至50 000时步时,边坡剪切带已经贯穿坡体,边坡的位移量达到1.39 m,由于边坡的位移巨大,滑带上侧岩体出现架空现象,滑带呈现出下小上大的塔式结构。

由于坡体呈现松散结构的特点,对于浸水不敏感,但坡体内易于产生拉裂变形,在防治时应注重前缘的压脚、挡护与防排水措施,且前缘变形易导致坡表出现张拉裂缝,因此在滑坡前易于监测。

图7所示的边坡矢量表明,坡脚处颗粒的滑动趋势远大于滑面其他部位,滑动方向整体朝下,且坡脚有显著的隆起趋势;图8为边坡平行黏结图,坡体黏结沿滑面附近的失效严重,从整体上看,由于坡体下部位移率先发生滑动且位移量大,下部坡体颗粒间黏结失效范围大,坡体上部颗粒间黏结完全失效,岩体脱离坡体形成裂隙,坡体滑面已经沿坡脚剪出,从侧面验证坡体已经进入失稳状态。

图7 边坡位移矢量

图8 边坡平行黏结状态

坡脚对边坡的滑动起着关键的制约作用,坡脚隆起的趋势随计算时步的增加越来越明显,并且由于坡脚的约束作用,致使此上部岩体产生向临空方向的弯曲变形。局部并出现微弱的剪裂、架空现象,进而促进了边坡向深部的变形与破坏,在上覆压力作用下,变形继续加剧,从而使岩体的破坏也进一步发展,当坡脚岩体被剪断,滑面上部岩体缺少下部的支撑,上部岩体完全脱开,滑面完全贯通。

4 结 论

本文基于某电站边坡实例开展颗粒离散元数值模拟研究,探讨了边坡的稳定性问题及边坡变形失稳特征,得到以下主要结论:

(1) 采用颗粒流法分析边坡稳定性,无需假定边坡滑面,并且可以得到边坡破坏的全过程,可以较好的得到边坡失稳破坏各部分的位移情况,对关键块体的分析及处理有着预警作用。

(2) 坡体呈现松散结构特征,对于浸水不敏感,但坡体内易于产生拉裂变形,运行期应注重坡体前缘的压脚档护以及防排水措施,且坡体前缘变形易引发后缘坡表出现张拉裂缝,因此,边坡发生滑坡前易于监测。

(3) 通过对坡体失稳后的位移及矢量的变化分析得到,坡脚对边坡的稳定起着重要作用,故而施工期应尽量避免切脚开挖。