月牙山坠落式危岩稳定性离散元数值分析

张 研，吴康丽，王 伟，高 宇

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

随着经济的发展，我国西南山区建设速度呈逐年增长趋势，但由于自然地理条件、地质构造和地层岩性等复杂因素的影响，使得地质灾害成为威胁西南山区居民生命财产、城镇建设和交通运输安全的首要因素。危岩崩塌属于三个地质灾害种类中重要的一种，其破坏性虽比地震、火山要小，但由于其发生的频率和广度远超其它灾害，故其危害性不容小觑[1]。岩体的稳定性常常由岩块之间的裂隙、节理、断层等不连续结构面所决定。纵横交错的结构面将岩体进行切割，形成我们日常所见的不完全连续体，并可逐渐形成为危岩体[2]。

目前，针对危岩稳定性分析及演化所采用的研究方法主要有实验法、解析法和数值模拟法。其中，实验法最为直接，但是具有周期长、成本较高和偶然性大等诸多缺点。又因危岩稳定性受多种因素影响，解析法无法将所有影响因素进行系统分析，故其计算结果代表性不佳。随着计算机技术和理论的发展，使得数值模拟技术已经广泛应用于危岩稳定性问题研究[3]。目前，主要有有限元法、有限差分法、离散元法等。有限元法和有限差分法主要以连续介质为基础，而危岩的崩塌过程大多是不连续的，这使得在实际工程中无法求解有初始间隙的相互接触体；并且以单元内应力作为屈服依据，可能会导致在计算位移和应力时出现计算失真[4]。离散法允许岩体之间出现滑动、平移、转动和岩体断裂等复杂过程,可用于模拟岩块产生大变形、节理面发生滑动甚至滑落等危岩稳定性问题时，不会对模型产生较大影响而产生误差，此外还具有识别新的接触面、重构计算模型进行分析等诸多优点[5-6]。

本文通过对桂林月牙山进行现场地质勘察，共发现10处群体危岩，其中W1危岩为坠落式危岩。从岩体规模和所处位置的高度来看，其体积较山体北侧的其他危岩体都要大，一旦发生危岩崩塌，W1坠落式危岩体所造成的影响将为最大。本文选取W1坠落式危岩为研究对象，运用UDEC离散元软件对其稳定性进行定性评价，揭示了危岩体的变形发展趋势，为合理选取工程防治措施提供参考。

1 月牙山区域地质概况

1.1 工程概况

桂林月牙山位于七星景区内，为典型的喀斯特岩溶地貌，山体主要由石灰岩组成，岩体长期裸露，受降雨等外界环境侵蚀，岩体表面裂隙发育丰富(见图1)。随着时间推移，风化、水侵蚀等作用使得裂隙逐渐扩大、加深，在危岩体自身重力或外力影响下，随时会发生崩塌落石灾害。

图1月牙山岩体裂隙

1.2 危岩分布概况

月牙山主体由硬质的石灰岩岩体、岩块组成，经过初步的野外地质勘察，共发现危岩体10处(见图2)，分别有倾倒式危岩、滑移式危岩、坠落式危岩，其中倾倒式危岩2处，滑移式危岩5处，坠落式危岩3处，月牙山危岩体具体分布情况见表1。

图2月牙山初步勘察危岩分布概况

表1 月牙山野外初步勘察危岩体资料表

2 月牙山危岩成因分析

2.1 地形地貌特征

桂林地区为典型的喀斯特地貌，受岩溶地质作用十分明显，山体岩性多为上泥盆统融县组(D3r)石灰岩，石灰岩属碳酸盐岩，坚硬、质纯、厚层，由于燕山运动的原因，境内山体多高耸直立、陡峭挺拔，山体上植被发育茂密，亦有石灰岩大量裸露，容易发生危岩崩塌灾害[7-8]。

2.2 水文气象条件

桂林所处纬度较低，接近北回归线，属亚热带季风性气候。全年气候温和，降雨量大且集中于夏季，夏长冬短且雨热相随，年平均相对湿度为73%～79%。该地区内地表河径流丰富，集雨面积较大且达到100 km2及以上的河流有六十余条，地下河流的分布亦相当丰富，平均每年的水分蒸发量约为1 490 mm～1 905 mm。因此，该地区大部分石灰岩山体长期受到降雨、雾气和地下水活动的影响作用，产生危岩崩塌的可能性更大。

2.3 人为因素影响

月牙山位于著名的七星景区内，游客数目巨大，频繁的人类活动对自然生态和地理环境造成一定程度的影响。工业和城市化进程的不断发展也给空气带来了越来越多的工业废气和汽车尾气，多种物质经过一系列复杂的化学反应后会产生SO2-和NO3-等酸根离子，这些离子与空气中的水蒸气结合后形成酸雨。酸雨的出现，使得岩石中的一些矿物在酸性条件下发生化学分解，导致原有的岩石结构破坏，从而降低岩石的力学性能，加大危岩崩塌、滑落的风险[9]。

3 离散元法基本原理

目前在工程领域广为应用的数值模拟方法主要有有限单元法、有限差分法和离散单元法。有限元法(如ANSYS、ABAQUS)虽然应用最广，但不适用于大变形、不连续体等问题的模拟；有限差分法(如FLAC3D)存在计算边界和单元网格划分的随意性；不连续变形分析法(如DDA)亦有参数选取的局限性[10-11]。离散元法除具有将岩体视为非连续结构体、允许发生接触位移、模拟过程更加真实等优点之外，还包含丰富的材料本构关系和接触本构关系[12]。是研究和分析非连续岩体稳定性问题的优选方法。

基于牛顿第二定律建立起来的离散元法与其他数值模拟方法相比，不需要复杂的数学推理过程，只需满足物理方程和运动方程：

(1) 物理方程。假定块体之间的法向作用力Fn和法向位移Wn成正比，则表达式为：

Fn=Kn·Wn

(1)

式中：Kn为界面的法向刚度系数，GPa/m。

块体与块体之间剪切力用剪切增量ΔFs表示，则表达式为：

ΔFs=Ks·Ws

(2)

式中：Ks为界面的切向刚度系数，GPa/m；Ws为块体之间的相对位移，m。

(2) 运动方程。由于块体之间是相互接触的，因此每个块体上受到不同角度的力。依据公式分别计算出作用在块体X方向、Y方向上的合力FX、FY和合力矩M。

(3)

(4)

μ=F合/m

(5)

式中：X0为质点横坐标；Y0为质心纵坐标；μ为加速度，m/s2；F合为合力，N；m为块体质量，kg。

美国ITASCA软件公司开发的二维离散元法UDEC软件利用显式差分法基本原理，可以为不稳定物理过程提供稳定解和为分析精度提供基本技术保障[13]。应用离散元法模拟分析时，岩体被当成仅由块体和接触面组成的模型，岩体整体模型设置有块体模型和节理模型两种，可以根据实际情形进行相应设置。UDEC离散元法程序软件中有两个较为常用的裂隙命令：Jset和Crack。Jset命令可以根据现场勘查所得的岩体节理参数资料赋值建立岩体裂隙，一般用于大量节理组的建立，并且可根据实际模型情况设置节理相关参数(如起始坐标、倾角、间距、岩桥等)，准确建立节理在模型中分布的位置和情况。Crack命令通过对任意两点坐标进行设置，两点坐标的直线连线即为岩体中的裂隙，这一命令常适用于建立单一节理面岩体的裂隙[14]。本文UDEC软件进行危岩体稳定性分析，具体实现步骤如图3所示。

图3 UDEC模拟计算流程图

4 工程实例分析

4.1 模型的建立

本文块体模型选用岩石力学界广为使用的摩尔-库仑本构模型[15]，掌握模型岩体实际的密度、体积模量、剪切模量、内摩擦角、黏聚力、剪胀角以及抗拉强度，并且根据这几个参数对模型岩体进行赋值建立初步块体模型。同时采用Jset命令建立模型的节理。选取W1坠落式危岩体为对象进行研究，危岩算例高约20 m，见图4，岩体本构模型采用摩尔-库仑模型，岩体顶部由于当地气候多雨的关系出现溶痕，渐渐发展成裂隙、裂缝(图4中黑粗线段表示裂隙)，岩体右侧下部受到差异风化的作用或是下部受到地下水侵蚀的作用逐渐被掏空，形成了此坠落式危岩模型。

图4山体模型南北向横截面图

在模型中，模型底部边界X轴和Y轴需要固定X、Y轴速度位移都为0，即：Xvel=0，Yvel=0；模型的左侧边界X方向速度需要边界约束，即：Xvel=0。

4.2 模拟参数确定

本文结合现场实地勘察和已有的桂林地区岩层研究成果获得了相关的岩石力学性质指标，结合地质勘察经验和监测分析，依据文献合理的选取参数数值[7-8,16]，见表2、表3。

表2 UDEC岩体模型石灰岩参数取值

表3 UDEC岩体模型节理参数取值

4.3 数值模拟结果分析

运算前将UDEC模型初始应力场调整为稳定的应力状态，山体的横截面右侧和山体底部设置为固定边界，在自然条件下仅仅受到重力作用，UDEC模型中重力的作用通过设置重力加速度达到，重力加速度数值取10 m/s2。在初始应力场中，最大主应力与最小主应力随着深度的变化而变化，在地势陡峭的地方主应力数值最大，而与岩体平行的地方主应力数值最小，几乎为零。靠近模型陡崖底部最大主应力和最小主应力会逐渐发生偏转，由竖直转为水平或由平行转为垂直。

模拟运算开始(见图5、图6、图7)，随着时间的推移，迭代步数的增加，在受重力的作用下，W1危岩体与山体之间的裂缝、裂隙逐渐变宽。由于此时危岩体与山体间出现摩擦和拉扯等力的作用，从初期应力场图可以看出应力变大且集中于裂隙下部或是危岩体与山体连接部，由初期位移矢量图可以看出裂隙、裂缝变宽后危岩体整体位移明显，岩体上部位移矢量值最大，岩体中部次之，位移矢量值大小由上至下依次减小。再由初期速度矢量图可以看出危岩体整体下落趋势趋于明显，尤其岩体上部加速度数值最大，岩体中部加速度数值次之，此时危岩体的潜在威胁变大。

图5 山体初期应力场图

图6 山体初期位移矢量图

图7山体初期速度矢量图

随着计算的进行(见图8、图9、图10)，迭代步数逐渐变大，后期山体与危岩体的裂隙、裂缝变得越来越宽，此时岩体之间的摩擦、拉扯现象更为剧烈，力的作用也更为明显。通过山体后期应力场图可以看出裂隙下部或者是山体与危岩体连接部的应力增大和集中现象也变得越来越明显，山体后期的位移矢量图中可以清晰的看到位移矢量方向较山体初期位移矢量图趋于向下的态势明显，而且危岩体上、中、下部位移矢量方向凌乱，与之前的初期位移矢量图相比位移矢量方向不均匀。山体后期的速度矢量图可以看出危岩体整体下落趋势已经十分明显，岩体上、中、下部的速度矢量方向凌乱但整体方向向下，危岩体此时的潜在威胁变得更大。

图8 山体后期应力场图

图9山体后期位移矢量图

图10山体后期速度矢量图

根据极限平衡法进行危岩体稳定性计算的结果为：W1岩体的稳定系数为1.46且小于1.5，说明该岩体处于欠稳定状态。当条件具备时，W1岩体极易发生坠落式崩塌破坏。由此可见，上述UDEC所模拟岩体变化发展趋势结果与该结论较为一致，且与实际相吻合，故可将计算结果应用于该危岩体稳定性的判定。同时，针对月牙山的其它危岩体稳定性也应及时做出判断和治理。

5 结 论

本文采用UDEC软件对月牙山坠落式危岩体(W1)稳定性进行数值分析，得出如下结论：

(1) 岩体的稳定性主要取决于结构面，危岩体在重力的影响下会使裂隙、裂缝下部以及山体与危岩体的连接处产生应力增大和应力不断集中的现象。

(2) W1危岩体已处于欠稳定状态，受危岩体所处的高度(重力势能)、裂隙深度的发展和外界因素等方面影响，随着时间的推移，该岩体终将会发生破坏。

(3) 将UDEC离散元法数值模拟分析软件应用于危岩崩塌模拟，结合工程实例经验和先前学者的研究，科学合理地选取岩体力学参数和调整软件相应设置。通过与极限平衡法计算结果进行对比，表明该方法对危岩稳定性的分析结果准确，可为日后危岩的防治提供有价值的参考数据。