观音岩滑坡潜在失稳模式及成灾规模数值模拟
在岩土工程和地质灾害研究领域，边坡失稳问题一直是一个备受关注的重要课题。传统的基于网格的数值模拟方法在模拟边坡失稳临界点时具有一定的优势，但在分析边坡滑移距离及成灾范围等大变形过程时，由于单元大变形引起的网格畸变问题，会导致数值系统能量发散，从而限制了完整模拟边坡失稳滑移全过程的能力。
一、PCMM方法简介
基本设想传统的基于网格的数值模拟方法在模拟边坡失稳临界点时具有一定的优势，但当分析边坡滑移距离及成灾范围时，传统的网格类方法当单元大变形时容易出现网格畸变，从而导致数值系统能量发散。
因此无法完整模拟边坡失稳滑移的全过程。颗粒离散元由大量随机分布的球形颗粒组成，主要用于固体材料细观力学行为的模拟。
颗粒离散元中颗粒间的接触信息及拓扑关系复杂而有序，如以颗粒间的接触关系为基础，将彼此接触的颗粒构成连续介质单元，并通过颗粒的接触状态实现单元的删除及重建，即可实现固体材料大变形、大运动过程的模拟。
基于上述想法，提出了一种基于颗粒接触的二维无网格方法(PCMM)。
PCMM方法采用显式欧拉前差技术，每一时步的工作主要包括:
（1）删邻居:将不再与本颗粒接触的邻居颗粒索引从本颗粒接触链表中删除；
（2）找邻居:重新寻找与本颗粒接触的邻居颗粒，并将邻居颗粒的索引增加至接触链表；
（3）删单元:根据本颗粒的当前接触链表，删除本颗粒单元链表中不满足条件的单元；
（4）找单元:根据本颗粒的当前接触链表重新创建连续介质单元；
（5）计算单元力:循环本颗粒单元链表，根据Mohr-Coulomb模型计算单元应力，转化为颗粒节点力；
（6）计算颗粒接触力:基于点-棱接触模型，计算颗粒与刚性线单元间的接触力；
（7）计算颗粒运动:根据牛顿定律计算颗粒的加速度、速度、位移等。
三角形创建的必备条件包括三条:
（1）构成三角形单元的三个颗粒必须彼此接触；
（2）所构建三角形单元的任何一个内角角度应在30°～150°之间(保证单元不会太奇异)；
（3）构成的三角形单元三条边的任意一条边不应该小于三个颗粒平均半径的0．5倍。在随机排列的颗粒体系中，为了构建稳定的连续介质单元系统，接触容差δ必须足够大以防止空隙的发生。
二、三角形单元的删除条件
经过一段时间的迭代计算后，单元将出现大变形，如果不进行任何修正继续计算，将出现由于单元畸变导致的系统发散等问题。
考虑到所有场量(如密度、速度、位移、应力、应变等)均位于颗粒上，且单元畸变后某一个方向尺寸非常小(单元体积可忽略)，将畸变后的单元删除不会引起系统总质量及总能量的改变。
因此，提出了如下的三角形单元删除原则(满足以下三个条件的任何一个：
（1）组成该三角形单元的三个颗粒不再彼此接触；
（2）该三角形单元的任何一个内角小于30°或者大于150°；
（3）该三角形单元的任何一条边长小于三个颗粒平均半径的0．5倍。
单元应力及颗粒合力求解三角形单元创建完毕后，即可采用经典连续介质理论进行单元应力及颗粒合力的求解。
采用基于FVM的增量方式进行相关计算，计算步骤包括:颗粒加速度求解→颗粒速度求解→颗粒位移增量求解→单元应变增量求解→单元应力增量求解→单元应力全量求解→基于塑性准则的单元应力修正→颗粒变形力求解→颗粒阻尼力求解→颗粒合力求解→下一时步颗粒加速度求解。
为了模拟边坡的失稳滑移过程，引入了理想弹塑性模型(含最大拉应力模型)。
三、观音岩滑坡数值分析
滑坡概况观音岩滑坡位于黔江河左岸，起点于污水处理厂大桥位置，终点于峡口观音庙位置，地形上基本呈一单面折线坡型，总体上呈上陡中缓下陡的地形，上部地形坡角25°～30°，中部地形坡角为10°～20°，下部地形坡角15°～25°。
地形总体北东高，南西低；中部与上部地形陡缓交接的部位基本为滑坡的后缘地带，地形标高在687．0m左右；
坡脚为黔江河，为滑坡的前沿剪出口位置，标高536．73m～551．92m；整个滑坡最大高差150m左右。
该滑坡平面形态近似长方形，沿河方向上长480m左右，沿滑动方向上宽190m左右，面积约69736m2。
滑体厚5．0m～35．0m，滑体总体前沿和后缘薄，一般厚度5．0m～10．0m；中部厚，一般厚度15．0m～35．0m。
该滑坡平均厚度16．5m左右，总体积约115．0×104m3，主滑方向236°，属中深层横长式大型土质滑坡。
据调查访问，该滑坡近年来变形迹象明显，每年雨季后在滑坡后部均会出现裂缝，缝宽一般2cm～8cm。
