桂林翻山危岩稳定性评价的离散元分析

刘宝臣， 王良玉， 曾榕*， 张研， 吴康丽

(1.广西岩土力学与工程实验室, 桂林 541004; 2.桂林理工大学土木与建筑工程学院, 桂林 541004)

危岩崩塌是广西发生率最高(占45%)、分布最广、危害最大的地质灾害，长期对该地区的旅游景区、房屋建筑、公路铁路、矿山水库等正常使用构成严重威胁[1-4]。危岩崩塌灾害已成为岩溶区人民生命财产安全的重要威胁之一，但其稳定性评价水平远未能满足社会发展需求[5-6]。

朱大鹏等[7]通过理论分析和离散元模拟法(3DEC)模拟相结合的方法对边坡滑移运动进行了分析,从运动学角度，给出楔形体的防治建议；张光福等[8]基于3DEC离散元软件,完成了对煤层井壁稳定的仿真模拟分析，验证了3DEC离散元法比弹性力学连续介质法更适用于煤岩井壁稳定研究；刘其琛等[9]运用3DEC建立关岭大寨滑坡三维模型,结合强度折减法对滑坡进行模拟,得出暴雨发生前和发生后的滑坡安全系数并对滑坡失稳过程进行分析；Wu等[10]利用3DEC评估2009年莫拉克台风引发的仙都山岩石崩塌破坏，模拟块体表现出局部岩石的弹性行为，成功地模拟了滑坡起始时的楔体破坏;Wang等[11]采用3DEC对某岩质边坡开挖进行了数值模拟，通过在不同条件下的对比，得出不同条件下边坡可能的破坏模式和安全系数的变化规律；刘波等[12]采用3DEC程序中离散裂隙网络技术,建立了能反映节理裂隙分布特征的离散裂隙网络模型，为研究工程岩体的稳定性提供了数据基础；康志强等[13]通过高精度不接触边坡扫描仪对边坡中节理面信息进行采集，应用持平投影法搜索出边坡临空面中存在的关键块体并使用3DEC模拟出了关键块体的失稳破坏模式；段启兵[14]采用3DEC数值模拟软件,模拟分析了不同工作面推进速度下上覆岩层垂直位移、顶板的运动矢量、推进速度和上覆岩层应力与位移关系以及不同工作面推进速度下基本顶失稳运动情况。

基于此，现以桂林翻山为研究背景，选取翻山Ⅱ区内危岩作为研究对象，对山体北侧进行勘查和统计，共发现存有危岩21处。采用离散元数值模拟法，依据研究区内危岩发育特征和破坏情况，选取滑移式危岩WY1#和坠落式危岩WY14#为目标岩体，对目标岩体在重力和降雨作用下的破坏过程进行模拟并对模拟结果进行分析，最后提出合理的治理方案。以期为此类危岩的稳定性分析和防治提供理论支持,具有一定的参考价值和实际意义。

1 离散元基本原理

离散元是以牛顿运动定律为基础的新的数值分析方法[15-17]。节理和裂隙将完整的岩石切割成若干小块体，如图1(a)所示，取其中一个块体进行研究，周围的块体对该块体形成约束，每一组力为Fxi、Fyi,如图1(b)所示。假设块体无法达到平衡，那么，块体B将发生转动，周围的块体定会阻碍该块体的运动，则物理方程满足位移和力的作用规律，每一个块体都满足牛顿第二运动定律。

1～5为岩块编号图1 块体集合及作用于个别块体的力Fig.1 The collection of blocks and the forces actingon individual blocks

1.1 物理方程

通过物理方程可以求出单元之间的作用力，表达式为

FN=KSSN

(1)

式(1)中：FN为块体的法向作用力；KN为几面的法相刚度系数；SN为法向位移。

块体与块体之间的剪切力用剪切增量ΔFS来表示，表达式为

ΔFS=KSSS

(2)

式(2)中：KS为界面的切向刚度系数；SS为块体之间的相对位移。

1.2 运动方程

由于块体之间是相互接触的，因此每个块体上受到不同角度的力。通过运动方程可以确定单元体之间的相互受力状态。依据公式分别计算出作用在块体x方向、y方向上的合力Fx、Fy和合力矩M。

(3)

(4)

(5)

式中：x0为质点横坐标；y0为质心纵坐标；μ为加速度；F合为合力；m为块体质量。

求解块体质心在j方向上的速度与位移，方程为

(6)

式(6)中：t0为初始时间；Δt为计算时步；t1=t0+Δt。

离散元法采用显示差分方法求解，实现对物理非稳定问题的稳定求解，可以追踪记录破坏过程和模拟结构的大范围破坏，包括块体间的完全脱离。还能模拟岩体静态或者动态荷载下的受力和位移，从而正确的分析岩体的变形、破坏模式以及稳定性[18-19]。

2 滑移式危岩体建模及分析

2.1 模型建立

采用Mohr-Coulomb模型，选取WY1#为例进行滑移式危岩模拟分析。危岩体位于山体上部，分布于205～210 m高程，主崩方向260°，危岩高5 m、宽4 m、厚2.5 m，体积50 m3。岩体后缘贯通主控结构面长约4 m，与水平方向倾角为70°，未贯通裂隙长约1 m。WY1#危岩体现场图片如图2所示。

2.2 边界条件的确定

3DEC力学边界条件分为应力边界和位移边界。本例依据实际情况，采取接近真实的边界位移约束，用bound命令分别对底面、左侧和后侧进行位移限制。即对左侧和右侧进行x和y方向固定位移，底面x、y、z3个方向均进行位移固定。在进行模拟时，只考虑危岩的重力和降雨作用下的位移和应力变化，忽略山体应力历史构造应力、温差等外界因素的影响。建立山体模型，尺寸为11 m×4 m×18 m的山体模型，采用genedage命令进行单元划分，共生成25 676个四面体单元，选择con2(摩尔-库尔塑性模型)。Gravity命令施加岩体自重。山体模型图及边界固定情况如图3所示。表1～表3为模型参数取值。

图2 WY1#危岩体现场图Fig.2 Site map of WY1# dangerous rock mass

图3 山体模型图Fig.3 Mountain model diagram

表1 岩块物理力学参数表

表2 岩石力学参数值

表3 节理力学参数

2.3 参数强度折减

桂林降水丰富，雨水进入岩体的裂缝中，不仅会对岩体产生裂隙水压力，还会对裂隙面起到润滑作用，加快裂隙开裂。因此，在模拟计算中要考虑长期暴露在外界环境下岩体后缘主控裂隙面在雨水作用下发生的物理参数变化。取折减系数k=1.25为此次降雨模拟的折减系数，经过迭代计算可以得出危岩在考虑降水折减条件下不同时期的整体位移矢量图和变形情况云图，如图4～图6所示。

图4反映的是滑移式危岩WY1#在降雨和重力共同作用下，岩体从前期开始到后期整体散落的全过程。当处于稳定状态的危岩体转为滑动状态时，总位移量值将不再是一个固定值，而是处于不断变化的状态，这是危岩体发生破坏的一个特征。从图4(a)和图5(a)可以看出，山体在初期重力场作用下的位移情况，由于模型底部进行x、y、z三轴方向位移的固定，因此，底部位移值为0，山体上部则均出现不同量值的位移。运行至4 100步时，山体位移基本达到稳定状况，主要的位移出现在为岩体上。由图5(b)可知，危岩体的矢量箭头方向与山体坡度滑动面平行且指向斜下方。迭代至8 000步时，通过图6(a)可以进行危岩体剪切滑动面的确定，同时可以看出岩体出现开裂并产生滑移，滑移方向沿切面向斜下方滑动。运行至40 000步时，可以观察出危岩在下滑过程中出现岩体裂隙开裂，并发生岩体散落现象，此现象在后期运动中更加明显，通过图5可以观察出。图6(b)和图6(c)中可以得出山体最大主应力值为372.63 kPa,最小主应力值为22.169 kPa，并且由于岩石的各项异性，山体应力呈不均匀分布。

为了凸显降雨作用对岩体的影响，选取8 000步的降雨强度折减作用下的位移和矢量云图与只考虑重力作用下的图进行对比分析，如图7所示。

图4 不同步长位移云图Fig.4 Asynchronous long displacement nephogram

图5 4 100步位移云图和矢量图Fig.5 4 100 step displacement nephogram and vector diagram

对比分析可知，经过强度折减的危岩体在运行至同样迭代8 000步的情况下最大的位移值为1.950 7 cm，只在重力作用下的最大位移值为1.904 4 cm。随着迭代步数的增加而运行至岩体破坏后期，位移量差值愈加明显。由此可见，经折减的危岩体整体位移量比未经折减的位移量大，强度降低同时使得位移矢量随之增大，岩体下滑速度加快。所以，降雨会影响危岩的发育和破坏。

3 坠落式危岩体建模及分析

3.1 模型建立

选取坠落式WY14#危岩体作为研究对象进行坠落式离散元模拟。经测量，危岩体高18 m，宽8 m，厚6 m，体积864 m3，属于体积较大危岩体。依据实测尺寸进行山体模型建立，整体模型建立，整体模型为45 m×63 m×8 m。岩体后缘贯通主控结构面长约10 m，与水平方向倾角为85°，未贯通裂隙长约8 m。山体图片与初期建模图如图8和图9所示。

图6 8 000步节理剪切位移和大小主应力云图Fig.6 Nephogram of joint shear displacement and magnitude principal stress in 8 000 steps

图7 8 000步位移矢量云图Fig.7 8 000-step displacement vector cloud diagram

图8 山体勘察图Fig.8 Mountain survey map

图9 初期建模图Fig.9 Initial modeling diagram

3.2 模型计算及结果分析

图10为坠落式危岩体的破坏过程云图，随着时间的推移，岩体实现了从潜在滑裂面的开裂—滑移坠落全过程。

图10 位移矢量云图Fig.10 Displacement vector nephogram

通过图10可以看出在模型迭代初期，由于山体受自重应力的作用使山体出现正体下移，由于滑裂面的存在，危岩体步位的位移下移，由于滑裂面的存在，危岩体部位的位移相对山体整体位移展现为较大值。运算至9 500步时，山体自身位移达到平衡稳定值，而危岩体继续保持下降趋势，出现了较为明显的错动，错动量值为0.4 m，此时危岩体尚处于一个黏结良好且相对完整的状态。当继续运行至40 500步时，危岩体上的层理之间则出现了分离和错位，最下部块体由于受上步块体的压力，相对于上部块体位移速度稍快。到达90 500步时，最初相对完整的危岩变为碎散状态并脱离主控结构面，完成了11.8 m的坠落位移，并会沿着山体边坡滚落至山脚，造成严重的危害。在整个坠落过程中，岩体失稳方向向下，沿着主控结构面滑移剪切破坏。

图11为运行至中期的山体最大主应力图和最小主应力图，从图11中可以看出，此刻的最大主应力为0.4 MPa，最小主应力为0.196 MPa。图11中应力分布不均匀，多集中在层理较为发育的部位。在山体与危岩体最底部临空面连接部位则出现了压应力集中区。

综上，可以分析出WY14#危岩体处于不稳定状态，随着时间的变化和外界荷载等因素的影响会促使此危岩发生坠落现象。

图11 21 500步主应力云图Fig.11 21 500 step principal stress nephogram

4 治理措施

危岩的主动防治措施主要包括：支撑、锚固、清除、嵌补等，被动防治措施包括拦石墙、柔性防护网和森林防护等。翻山Ⅱ区内共21处危岩体，6处坠落式，15处滑移式，如表4所示。

危岩的治理方法有3种：危岩加固、危岩避让和危岩清除。对于WY1#-WY13#、WY15#-WY21#危岩体，通过分析工程量的大小和成本等因素的影响，选择主动支护中的清除法，清除法选择静态爆破。静态爆破震动小，对周围的岩体影响小，治理效果好，由于施工场地范围内现为拟建公园园区，附近居民及行人较少，山体下方为大片空地，施工对周围环境影响较小，因此采用静态爆破清除最为适宜。对实施静态爆破厚的周边岩体采用灌浆法进行加固，减少雨水进入裂隙而对危岩产生孔隙水压力，防止危岩裂隙开裂速度加剧，从而形成新的危岩。

静态爆破主要施工工序：设置防护措施→搭设施工平台→危岩体钻孔→安装静态爆破剂→爆破→处理小块危岩→现场破碎岩体运出→验收→清理场地、退场。

由于危岩体体积较大，被动支护方式防护能力有限，因此采用主动支护方式。根据实际工程地质情况和支护环境，对WY14#设计最终方案时采用锚固设计。

表4 翻山危岩体分布详情统计表

经计算，选用10根锚固分5排、2列进行支护，排间距3 m，列间距4 m。每根锚固杆由3根直径32 mm的HRB400钢筋焊接组成，采用全场黏结支护方式进行锚固，锚杆长15 m。3DEC中设有钢索(CABLE)结构单元用于实现岩体的锚固单元施加。经计算和查取相关资料进行锚固工程相关参数的取值，如表5所示。

为了更好地监测和之后的支护对比，在此模型中设置12个监测点进行位移监测，监测点及模型初期如图12和图13所示。

通过模拟计算，可以分别得出危岩在支护前和采取支护后的位移变化折线图，折线图可明显地将各监测点的位移前后变化信息直观地展现出来，结果如图14所示。

采用记录数据(Hist)命令对危岩后缘裂隙面及潜在裂隙面设置1～10号监测点，同时又在岩体最外测边缘处最上端和最底端设置11号监测点和12号监测点，对危岩体坠落过程中产生的位移进行记录。1～10号监测点分别对应图14中的Histn=62～72，11～12对应Histn=72～73。其中1、2号监测点位于裂隙顶端，由于在计算过程中1、2号监测点受裂隙左侧的山体位移限制，位移值为0，故将其舍弃不进行分析。因此，在图中只显示记录Histn=64～73变化趋势。

表5 锚固参数表

图12 模型计算图Fig.12 Model calculation diagram

图13 锚杆支护分布图Fig.13 Distribution diagram of bolt support

图14 支护前后监测点位移Fig.14 Displacement of monitoring points before and after support

通过图14(a)可知，未经过支护处理的危岩，经过迭代计算，各点位移量值变化呈递增趋势。前期由于重力作用下，潜在破裂面被剪断出现滑移现象，监测点的位移变化趋势较为平缓；随着滑移坠落位置的改变，岩体重心位置不断发生变化，同时山体层理面也在力的作用下发生开裂，对监测点监测的量值产生影响，出现波动现象。

通过图14(b)可知，在采取支护措施后，山体只在前期出现0.4 m的位移量，运行至9 000步时，锚杆锚固作用被触发。由图中曲线可知，后期的Histn=64～69和n=72的位移值增长幅度较小，变化基本趋于直线。由此说明，锚杆对危岩体的锚固起到很好的作用，如7号监测点(Histn=68)在step 15 000未进行支护时位移量约为1.2 m，在进行支护后位移量被控制在0.4 m左右，位移量降低了66.7%；在25 000步未进行支护时位移量为1.35 m，支护后的位移值仍为0.4 m，位移量降低70.3%，支护效果愈加明显。

图15 锚杆图Fig.15 Anchor bolt diagram

由于n=71、73记录的是10、12号监测点的位移，而在进行锚杆支护时，最底部体积较小的岩块11 697块体未能进行锚杆支护，因此在自身重力作用下仍处于坠落状态。因此，对于这种情况下的岩体，除了进行锚杆支护外，还应对该碎石进行人工清除或在下部设置防护网，对碎石进行防护。

锚杆受力示意图如图 15(a)所示，图15(b)可以看出锚杆的最大位移出现在危岩体内，中间节理部位处的锚杆出现了明显的方位变化，这是由于锚杆起作用时，因主控结构面右侧的危岩体与 山体形成错动，使得在主控结构面部位处的锚杆承受较大的剪应力，从而在主控结构面 附近产生较大的几何错位和锚杆的径缩现象。危岩的位移值除了钢筋受拉弹性阶段产生的位移值外，还包括灌浆体与岩石和钢筋之间的相对滑移。当荷载逐渐增大到锚杆一定承载力，灌浆体与周围介质间的滑移将停止，否则锚固 杆体就会发生破坏。因此，经过支护的位移图后期变化较为平缓。

锚杆使得山体与危岩重新组合成为一个整体并对正在滑移坠落的岩体提供一定的抗滑力(支护力)，抗滑力将下滑力抵消且阻止了岩体持续破坏的趋势，提高了岩体的安全系数，使得危岩体稳定性增强。综合图15(a)和图15(c)可以看出，锚杆所受轴向拉力范围和产生位移部位集中分布在主控结构面处和结构面两侧1倍的锚杆直径范围内。而锚杆产生的力会作为附加应力作用在结构面上，提高结构面的刚度。危岩完成加锚支护施工后，沿主控结构面布设的监测点的位移值都比无锚支护时小。锚杆的支护起到了很好的控制岩体发生破坏的作用，限制了岩体的变 形。因此，采用此支护方式是可行的。

5 结论

对翻山Ⅱ区危岩体进行了实地勘察，并依据危岩破坏模式进行稳定性分析，综合判定危岩的稳定状态。在3DEC离散元模拟软件中建立翻山WY1#和WY14#危岩体的三维计算模型，模拟降雨对危岩稳定性的影响。通过模拟结果可以得出，岩体处于不稳定状态，对危岩稳定性做出评价后，对岩块体积小的危岩采用静态爆破的方法进行清理，对于体积较大的危岩选用锚杆支护的方法进行加固。危岩治理不是一劳永逸的工程，还需要长期的监测，做到早发现早治里，保障周边环境的安全。