滑坡堵河冲击大坝区域建筑物过程的模拟预测方法及系统

技术领域

本发明属于滑坡模拟预测技术领域，具体涉及滑坡堵河冲击大坝区域建筑物过程的模拟预测方法及系统。

背景技术

水库库区滑坡是水利水电工程中常见的地质灾害，目前由于气候变化影响，强震、暴雨等极端事件发生概率增加，极易诱发滑坡灾害，冲击和掩埋建筑物，威胁群众生命财产安全。滑坡的有效防治取决于对其触发及传播动力过程的正确认识，包括滑坡演进过程及最终堆积形态等，前者主要影响滑坡对工程结构的冲击动力响应，后者对滑坡影响范围分析、危害性评估具有重要意义。

对于滑坡冲击研究，现有研究多集中在模型试验与数值模拟两方面，模型试验能较好地揭示滑坡传播侵蚀及冲击特征，但由于其成本相对较高，且滑坡运动属动力过程，试验条件复杂且难控制。因此数值模拟逐渐成为开展此类研究的重要手段。现有研究方法主要可分为三类，基于网格的方法如有限单元法，基于粒子的方法如离散单元法、光滑粒子流体动力学方法等，以及混合方法如物质点法、粒子有限元法等。传统的基于网格的方法存在计算滑坡大变形时网格易畸变导致计算失真，耦合计算流体力学的CFD-DEM方法存在DEM参数标定困难和计算尺度限制等缺点。光滑粒子流体动力学方法目前在不断发展中，逐渐应用与滑坡大变形模拟中，但目前仍存在三维精细化建模难度大、边界处理困难等问题，仍在不断发展中。此外，现有学者和工程界多关注单灾害的研究，对灾害链研究较少。目前研究多关注滑坡触发过程或者演进过程，两者是割裂的。目前对于滑坡堵河冲击大坝附属建筑物灾害链的模拟仍存在难度，其过程涉及滑坡触发、流固耦合堵河过程、侵蚀、冲击建筑物等复杂过程，采用一种数值模拟方法在统一框架下实现仍存在较大难度。由于灾害链演化过程的复杂性，目前对于灾害链过程较为完善的数值模拟方法及手段仍有待进一步发展。

发明内容

本发明是为了解决上述问题而进行的，目的在于提供滑坡堵河冲击大坝区域建筑物(大坝区域建筑物包含大坝及其附属建筑物)过程的模拟预测方法及系统，能够用于模拟预测强震或暴雨等作用下潜在滑坡体触发、演进、侵蚀沿途土体及建筑物、冲击下游建筑物等全过程，动态计算滑坡体的演变情况，量化滑坡灾害。

为了实现以上目的，本发明采用了以下方案：

<方法>

本发明提供滑坡堵河冲击大坝区域建筑物过程的模拟预测方法，其特征在于，包括以下步骤：

步骤1，根据待模拟的情况，通过建模软件建立相应的滑坡土、地形边界、大坝区域建筑物、岩石及水流三维实体模型，并离散为粒子模型，分别对滑坡土粒子、地形边界粒子、建筑物粒子、岩石(块石)粒子及水流粒子这五类粒子赋予初始参数值和用于区分各类型粒子的粒子属性信息；

步骤2，在步骤1基础上，构建滑坡触发及传播灾害链全过程模拟框架；

步骤2-1，滑坡触发前的模拟；

在滑坡触发前，对于未发生滑坡的滑坡土粒子i，采用固体弹塑性本构模拟其变形行为，并根据如下公式计算特性信息；特性信息至少包含应力、应变、加速度、速度、位移、密度信息；

式中：ρ为粒子的密度；t为计算时间；m为粒子的质量；v为粒子的速度；x是基点的位置坐标；W是光滑核函数，α与β表示应力的向量方向，应力为三阶张量，α与β均可对应取x、y、z三个方向；S为滑坡土粒子区域；B为边界粒子区域；在影响域内，以滑坡土粒子i为中心，周围的粒子j有的属于S，为滑坡土粒子，有的属于B，为边界粒子；σ

步骤2-2，根据滑坡触发临界转化的判断条件确定是否触发滑坡；

步骤2-3，滑坡触发后冲击沿途岩石及建筑物过程的模拟；

滑坡触发后，采用如下公式求解发生滑坡的滑坡土体粒子i的特性信息，并将当前滑坡土体粒子作为失稳滑坡土粒子用于后续侵蚀计算；

式中：S为滑坡土粒子区域，W为水流粒子区域，B为边界粒子区域，R为岩石或者建筑物粒子区域，p

步骤2-4，岩石或建筑物的接触力计算；

岩石或建筑物a由一系列相应的粒子k组成，基于如下公式计算粒子速度及接触力：

式中，岩石或建筑物a具有质量M

步骤2-5，滑坡土及岩石建筑物侵蚀启动判断；

滑坡触发后失稳滑坡土会对沿途的滑坡土、岩石及建筑物造成冲击，若达到侵蚀启动条件，则会带走滑坡土、岩石、建筑物一起移动；侵蚀包括对沿途土体的侵蚀和对沿途岩石及建筑物的侵蚀；对待判断的滑坡土粒子，通过临界剪切强度确定是否触发侵蚀启动，若触发侵蚀启动，则该滑坡土粒子被侵蚀，成为被侵蚀的滑坡土粒子；对于待判断的岩石或建筑物粒子，通过冲击碰撞力确定是否触发侵蚀启动，若触发侵蚀启动，则该岩石或建筑物粒子被侵蚀，成为被侵蚀的岩石或建筑物粒子；被侵蚀的滑坡土粒子作为失稳滑坡土粒子在下一个时间步计算过程中，采用步骤2-3的公式进行计算模拟；

步骤2-6，基于步骤2-1～步骤2-5计算得到当前时间步内各粒子特性信息和侵蚀情况数据，然后根据当前时间步计算得到的数据，确定下一个时间步计算用的初始参数值；重复前述过程，直至达到预设时间步数；

步骤3，根据步骤2构建的滑坡触发及传播灾害链全过程模拟框架，预测潜在滑坡土触发、演进、侵蚀沿途土体及建筑物、冲击下游建筑物全过程，预演滑坡触发及传播过程，量化滑坡灾害。

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测方法，在步骤2-2中，滑坡触发前，采用强度折减法计算滑坡土变形，得到塑性区贯通时候的强度折减参数，记录此时刻的折减系数RFS

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测方法，在步骤2-5中，当失稳滑坡土冲击碰撞力小于岩石或建筑物整体启动力时，则岩石或建筑物不发生整体侵蚀启动；当失稳滑坡土冲击碰撞力大于岩石或建筑物整体启动力时，则岩石或建筑物整体侵蚀启动发生，开始随着失稳滑坡土一起运动；此处冲击碰撞力通过冲击力响应计算算法得到：首先确定需要计算的目标结构(岩石或建筑物)的组成粒子信息，包括单个粒子k的质量、组成该建筑物的粒子的数量N；然后，通过公式2-3-2对影响域内位于粒子k周围的其余粒子进行插值计算，求得组成目标结构的所有粒子的加速度，其余粒子类型包括：滑坡土粒子、水流粒子、地形边界粒子、岩石粒子、建筑物粒子；对于每个粒子k，均根据x、y、z三个方向求得的分力进行

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测方法，在步骤2-5中，在失稳滑坡土冲击碰撞力小于岩石或建筑物整体启动力的情况下，岩石或建筑物还可能发生局部侵蚀，发生局部侵蚀时，被冲击碰撞的岩石或建筑物会发生碎裂解体效应，碎裂后的岩石粒子或建筑物粒子会被失稳滑坡土侵蚀带走；

①采用如下方法判断岩石或建筑物粒子是否发生局部侵蚀：

分别计算粒子的拉伸应力和剪切应力，粒子的拉伸应力计算公式为：

剪切应力计算公式为：

式中，F

②计算得到粒子的拉伸应力和剪切应力后，进行断裂损伤破坏判断，满足断裂破坏准则，则认为该粒子发生了局部侵蚀导致的损伤破坏并成为损伤粒子；

③破坏后粒子应力更新及损伤统计：

统计损伤粒子的数量，用于损伤进行实时量化，并将损伤粒子的抗拉强度σ

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测方法，在步骤2-5的②中，断裂破坏准则采用拉伸和剪切两种破坏模型进行判断，公式如下：

σ

t

式中，σ

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测方法，在步骤2-5中，当沿途的滑坡土粒子所受到的冲击剪切力大于临界剪切强度值时，当前滑坡土粒子被侵蚀铲刮带走，成为被侵蚀的滑坡土粒子，否则当前滑坡土粒子仍处于静止状态，采用Drucker-Prager屈服准则判断土体是否屈服。

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测方法，在步骤3中，根据滑坡土粒子、岩石粒子的特征信息和岩石及建筑物的侵蚀情况量化滑坡灾害链增容效应和灾害放大效应；增容效应量化包括：量化被侵蚀的滑坡土粒子和发生整体或局部侵蚀的岩石与建筑物的相应组成粒子使失稳滑坡土的体积增加的部分；灾害放大效应量化包括：量化发生侵蚀后滑坡的冲击速度和对下游建筑物的冲击力的变化情况。

<系统>

进一步，本发明还提供了能够自动实现以上<方法>的滑坡堵河冲击大坝区域建筑物过程的模拟预测系统，其特征在于，包括：

建模部，根据待模拟的情况，通过建模软件建立相应的滑坡土、地形边界、大坝区域建筑物、岩石及水流三维实体模型，并离散为粒子模型，分别对滑坡土粒子、地形边界粒子、建筑物粒子、岩石粒子及水流粒子这五类粒子赋予初始参数值和用于区分各类型粒子的粒子属性信息；

模拟框架构建部，在建模部的基础上，构建滑坡触发及传播灾害链全过程模拟框架；包括：

滑坡前模拟单元，滑坡触发前的模拟；在滑坡触发前，对于未发生滑坡的滑坡土粒子i，采用固体弹塑性本构模拟其变形行为，并根据如下公式计算特性信息；特性信息至少包含应力、应变、加速度、速度、位移、密度信息；

式中：ρ为粒子的密度；t为计算时间；m为粒子的质量；v为粒子的速度；x是基点的位置坐标；W是光滑核函数，α与β表示应力的向量方向，应力为三阶张量，α与β均可对应取x、y、z三个方向；S为滑坡土粒子区域；B为边界粒子区域；在影响域内，以滑坡土粒子i为中心，周围的粒子j有的属于S，为滑坡土粒子，有的属于B，为边界粒子；σ

触发确定单元，根据滑坡触发临界转化的判断条件确定是否触发滑坡；

滑坡后模拟单元，滑坡触发后冲击沿途岩石及建筑物过程的模拟；滑坡触发后，采用如下公式求解发生滑坡的滑坡土体粒子i的特性信息，并将当前滑坡土体粒子作为失稳滑坡土粒子用于后续侵蚀计算；

式中：S为滑坡土粒子区域，W为水流粒子区域，B为边界粒子区域，R为岩石或者建筑物粒子区域，p

接触力计算单元，计算岩石或建筑物的接触力；岩石或建筑物a由一系列相应的粒子k组成，基于如下公式计算粒子速度及接触力：

式中，岩石或建筑物a具有质量M

侵蚀判断单元，滑坡土及岩石建筑物侵蚀启动判断；滑坡触发后失稳滑坡土会对沿途的滑坡土、岩石及建筑物造成冲击，若达到侵蚀启动条件，则会带走滑坡土、岩石、建筑物一起移动；侵蚀包括对沿途土体的侵蚀和对沿途岩石及建筑物的侵蚀；对待判断的滑坡土粒子，通过临界剪切强度确定是否触发侵蚀启动，若触发侵蚀启动，则该滑坡土粒子被侵蚀，成为被侵蚀的滑坡土粒子；对于待判断的岩石或建筑物粒子，通过冲击碰撞力确定是否触发侵蚀启动，若触发侵蚀启动，则该岩石或建筑物粒子被侵蚀，成为被侵蚀的岩石或建筑物粒子；被侵蚀的滑坡土粒子作为失稳滑坡土粒子在下一个时间步计算过程中，采用公式式2-3-1～2-3-2进行计算模拟；

迭代单元，基于滑坡前模拟单元、触发确定单元、滑坡后模拟单元、接触力计算单元、侵蚀判断单元计算得到当前时间步内各粒子特性信息和侵蚀情况数据，然后根据当前时间步计算得到的数据，确定下一个时间步计算用的初始参数值；重复前述过程，直至达到预设时间步数；

模拟预测部，采用滑坡触发及传播灾害链全过程模拟框架模拟预测潜在滑坡土触发、演进、侵蚀沿途土体及建筑物、冲击下游建筑物全过程，预演滑坡触发及传播过程，并量化滑坡灾害；

控制部，与建模部、模拟框架构建部、模拟预测部均通信相连，控制它们的运行。

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测系统，还可以包括：输入显示部，与控制部通信相连，让操作员输入操作指令，并进行相应显示。

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测系统，在侵蚀判断单元中，当失稳滑坡土冲击碰撞力小于岩石或建筑物整体启动力时，则岩石或建筑物不发生整体侵蚀启动；当失稳滑坡土冲击碰撞力大于岩石或建筑物整体启动力时，则岩石或建筑物整体侵蚀启动发生，开始随着失稳滑坡土一起运动；此处冲击碰撞力通过冲击力响应计算算法得到：首先确定需要计算的目标结构的组成粒子信息，包括单个粒子k的质量、组成该建筑物的粒子的数量N；然后，通过公式2-3-2对影响域内位于粒子k周围的其余粒子进行插值计算，求得组成目标结构的所有粒子的加速度，其余粒子类型包括：滑坡土粒子、水流粒子、地形边界粒子、岩石粒子、建筑物粒子；对于每个粒子k，均根据x、y、z三个方向求得的分力进行

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测系统，在触发确定单元中，滑坡触发前，采用强度折减法计算滑坡土变形，得到塑性区贯通时候的强度折减参数，记录此时刻的折减系数RFS

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测系统，在侵蚀判断单元中，在失稳滑坡土冲击碰撞力小于岩石或建筑物整体启动力的情况下，岩石或建筑物还可能发生局部侵蚀，发生局部侵蚀时，被冲击碰撞的岩石或建筑物会发生碎裂解体效应，碎裂后的岩石粒子或建筑物粒子会被失稳滑坡土侵蚀带走；

①采用如下方法判断岩石或建筑物粒子是否发生局部侵蚀：

分别计算粒子的拉伸应力和剪切应力，粒子的拉伸应力计算公式为：

剪切应力计算公式为：

式中，F

②计算得到粒子的拉伸应力和剪切应力后，进行断裂损伤破坏判断，满足断裂破坏准则，则认为该粒子发生了局部侵蚀导致的损伤破坏并成为损伤粒子；

③破坏后粒子应力更新及损伤统计：

统计损伤粒子的数量，用于损伤进行实时量化，并将损伤粒子的抗拉强度σ

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测系统，在侵蚀判断单元的②中，断裂破坏准则采用拉伸和剪切两种破坏模型进行判断，公式如下：

σ

τ

式中，σ

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测系统，在侵蚀判断单元中，当沿途的滑坡土粒子所受到的冲击剪切力大于临界剪切强度值时，当前滑坡土粒子被侵蚀铲刮带走，成为被侵蚀的滑坡土粒子，否则当前滑坡土粒子仍处于静止状态，采用Drucker-Prager屈服准则判断土体是否屈服。

优选地，本发明提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测方法，在模拟预测部中，根据滑坡土粒子、岩石粒子的特征信息和岩石及建筑物的侵蚀情况量化滑坡灾害链增容效应和灾害放大效应；增容效应量化包括：量化被侵蚀的滑坡土粒子和发生整体或局部侵蚀的岩石与建筑物的相应组成粒子使失稳滑坡土的体积增加的部分；灾害放大效应量化包括：量化发生侵蚀后滑坡的冲击速度和对下游建筑物的冲击力的变化情况。

发明的作用与效果

1、本发明建立了考虑“滑坡-河流-地形-大坝区域建筑物”进行灾害链动力耦合计算的滑坡触发及传播灾害链全过程模拟框架，可用于计算强震或暴雨作用下潜在滑坡体发生大变形，堵塞河流并进一步冲击大坝附属建筑物的复杂灾害链动力过程，为潜在滑坡体风险评估和早期预警提供技术支持。

2、本发明不仅考虑了“滑坡土-河流-岩石-大坝区域建筑物”相互冲击作用，而且考虑了滑坡土对沿途滑坡土、岩石(块石)、大坝区域建筑物侵蚀作用，可用于预测强震或暴雨作用下潜在滑坡体触发、演进、侵蚀沿途土体及建筑物、冲击下游建筑物等全过程，动态计算滑坡体运动速度、滑移带涡度、堆积形态、对建筑物的冲击响应等，量化滑坡灾害链增容效应和灾害放大效应，可预演滑坡触发及传播过程。

3、本发明不仅考虑了滑坡土对岩石或建筑物的整体侵蚀作用，还考虑了滑坡土对岩石或建筑物的部分侵蚀作用，能够使得模拟预测过程更加符合实际情况，从而使得模拟预测结果更为准确。

4、本发明提出了滑坡触发及发生大变形冲击全过程的动态数值模拟手段，解决了实际工程中滑坡体演进过程难计算、堆积影响范围难预测、冲击响应难表征的难题。可用于预测和计算滑坡体运动速度、滑移带涡度、堆积形态、对建筑物的冲击响应等，预演滑坡触发及传播过程，对滑坡体潜在危险性进行预测和风险评估。

5、本发明可用于事先对滑坡体潜在危险性进行预测，在滑坡事故发生之前预先计算滑坡体触发之后达到下游建筑物的时间、到达速度及冲击力大小，为沿途及下游建筑物进行早期预警和预先采取防护加固措施提供指导。

附图说明

图1为本发明实施例涉及的实体模型转粒子模型算法示意图；

图2为本发明实施例涉及的滑坡建模图；

图3为本发明实施例涉及的接触计算中各粒子耦合机理图；

图4为本发明实施例涉及的土体侵蚀启动机制图；

图5为本发明实施例涉及的岩石及建筑物整体侵蚀启动机制图；

图6为本发明实施例涉及的建筑物受冲击力计算示意图；

图7为本发明实施例涉及的岩石及建筑物在滑坡冲击作用下相互碰撞、发生裂缝扩展及破裂部分侵蚀启动机制图，其中，(a)为岩石及建筑物还未发生侵蚀启动，(b)为岩石发生了整体侵蚀启动，(c)为侵蚀后的岩石碰撞建筑物产生裂缝，(d)为岩石何建筑物产生裂解形成碎块；

图8为本发明实施例涉及的不同折减系数下滑坡体塑性区分布情况图；

图9为本发明实施例涉及的不同折减系数下滑坡体特征点(顶点)处位移及增量变化；

图10为本发明实施例涉及的本发明方法计算模拟数据与现有试验实测数据冲击力定量对比图，其中，(a)为现有实验现场图，(b)为对比图；

图11为本发明实施例涉及的滑坡侵蚀沿途滑坡土的过程图；

图12为本发明实施例涉及的模拟滑坡侵蚀沿途碎石及启闭房的过程图；

图13为本发明实施例涉及的滑坡冲击过程中岩石及启闭房受力情况示意图；

图14为本发明实施例涉及的滑坡体(多层滑坡土形成)-地形-河道-大坝及其附属建筑物的三维精细化模型图；

图15为本发明实施例涉及的模拟滑坡体堵河冲击大坝区域建筑物的过程图；

图16为本发明实施例涉及的滑坡体最终堆积过程评估及预测情况示意图；

图17为本发明实施例涉及的大坝区域建筑物冲击力时程变化图。

具体实施方式

以下结合附图对本发明涉及的滑坡堵河冲击大坝区域建筑物过程的模拟预测方法及系统的具体实施方案进行详细地说明。

<实施例一>

本实施例所提供的滑坡堵河冲击大坝区域建筑物过程的模拟预测方法包括如下步骤：

步骤1，根据待模拟的情况，通过建模软件建立相应的滑坡土、地形边界、大坝区域建筑物、岩石及水流三维实体模型，并离散为粒子模型，分别对滑坡土粒子、地形边界粒子、建筑物粒子(形成各个大坝区域建筑物的粒子)、岩石粒子及水流粒子这五类粒子赋予初始参数值和用于区分各类型粒子的粒子属性信息。大坝区域建筑物包括大坝及其附属建筑物，大坝附属建筑物如启壁房、输水管道、闸门、溢洪道等。

首先可在基本建模软件如CAD或犀牛等建立滑坡体、地形三维实体模型，后使用采样算法将实体模型离散为等距离粒子模型，采样算法思路为：如图1所示，在实体模型内部设置一点作为种子生成点，并设置一范围包括了实体模型的计算域M，在此计算域内生成粒子，粒子生成使用判断标准如下：当粒子在实体模型内部边界N时则生成此粒子，当粒子不在实体模型内部边界N时则舍弃此粒子，其中M区域包含了N区域。通过设置粒子半径大小，滑坡体建模的精度可达厘米级，实现三维滑坡模型高保真建模，滑坡体建模过程如图2所示。本发明无需剖分网格，突破了传统网格法计算滑坡大变形发生网格畸变导致计算终止的局限，实现了三维滑坡体、建筑物及地形的精细化高保真建模，建模精度达厘米级别，滑坡体建模误差在5％以内，满足工程计算要求。

滑坡土粒子赋值信息主要包括固体本构参数信息(失稳前的那些粒子采用固体本构)及流体本构参数信息(失稳后或已知已成为流体状态采用流体本构)，具体包括密度、单个粒子体积、弹性模量、泊松比、粘聚力、内摩擦角、抗剪强度、强度折减参数以及流变参数、粒子类型(用于接触判断)。滑坡触发前，滑坡体基本保持静止，呈现固态，此时采用固体本构表征其运动特性。而在地震或暴雨等作用下，发生流态化，滑坡触发后运动，考虑到此时滑坡体已经表现出类似流体的运动行为特性，此时采用流体本构表征其运动特性。

地形边界粒子赋值信息主要包括摩擦系数、体积、半径、密度、弹性模量、泊松比、粒子类型(用于接触判断)。

建筑物粒子(包含形成大坝及其附属建筑物的粒子)赋值信息主要包括摩擦系数、体积、半径、密度、弹性模量、泊松比、接触刚度、粒子类型(用于接触判断)；对不同建筑物可以采用建筑物序号或标签进行区分。

岩石粒子赋值信息主要包括摩擦系数、体积、半径、密度、弹性模量、泊松比、接触刚度、粒子类型(用于接触判断)；对不同岩石也可以采用岩石序号或标签进行区分。

水流粒子赋值信息主要包括密度、单个粒子体积、半径、粘度、粒子类型(用于接触判断)。

步骤2，在步骤1基础上，构建滑坡触发及传播灾害链全过程模拟框架。

步骤2-1，滑坡触发前的模拟：

在滑坡触发前，滑坡土体仅小变形，对于未发生滑坡的滑坡土粒子i，仅需考虑与地形边界及同类型滑坡土粒子之间接触作用，采用固体本构模拟其变形行为，并根据如下公式计算特性信息；特性信息至少包含应力、应变、加速度、速度、位移、密度信息。

式中：ρ为粒子的密度；t为计算时间；m为粒子的质量；v为粒子的速度；x是基点的位置坐标；W是光滑核函数，α与β表示应力的向量方向，应力为三阶张量，α与β均可对应取x、y、z三个方向；S为滑坡土粒子区域；B为边界粒子区域；在影响域内，以滑坡土粒子i为中心，周围的粒子j有的属于S，为滑坡土粒子，有的属于B，为边界粒子；σ

步骤2-2，根据滑坡触发临界转化的判断条件确定是否触发滑坡：

采用强度折减法计算滑坡变形，折减过程如下，通过对土体强度参数进行不同程度折减，采用固体本构计算滑坡土体等效塑性区分布，通过计算不同折减系数下滑坡塑性区的分布，得到塑性区贯通时候的强度折减参数，此时为滑坡触发临界，记录此时刻折减系数RFS

滑坡触发临界转化条件，考虑地震以及暴雨条件下滑坡触发，此时土体强度降低，采用强度折减法对滑坡土体强度参数进行折减，公式如下：

式中，c

使用折减后的强度参数带入滑坡土体本构方程计算求得折减后的滑坡土体应力、应变、速度及位移信息之后，带入公式2-1-1和2-1-2求解下一计算步的应力、应变、速度及位移信息。以滑坡特征点处位移突变和塑性区贯通两个条件作为滑坡触发临界条件，必须同时满足滑坡才算触发。

步骤2-3，滑坡触发后冲击沿途岩石及建筑物过程的模拟：

本实施例中，关注滑坡触发后，强震或暴雨作用下潜在滑坡体发生大变形，堵塞河流并进一步冲击大坝附属建筑物的复杂灾害链演化过程。考虑到强震或暴雨作用下滑坡体已经表现出类似流体的运动行为特性，呈现出类似泥石流的非牛顿流体运动特征，可称之为流态型滑坡。因此采用流变本构模型模拟滑坡体触发后的运动过程。流变本构模型及其等效粘度系数公式如下：

τ

式中：γ为剪切应变速率；μ为表观动态粘度；τ

此时失稳后的滑坡土粒子在向下游冲击过程中会与滑坡自身同类型粒子、边界粒子、水流粒子、岩石及建筑物粒子进行接触。首先进行接触判断，以粒子影响域作为判断区域，以当前作为中心粒子的发生滑坡的滑坡土粒子i为研究对象，当周围粒子在滑坡土粒子i的影响域范围内，则存在相互接触作用，求解粒子特性信息，每一步计算处粒子的加速度、速度、位移、等信息，作为下一计算步的初始条件。

滑坡触发后，采用如下公式求解发生滑坡的滑坡土体粒子i的特性信息，并将当前滑坡土体粒子作为失稳滑坡土粒子用于后续侵蚀计算。

式中：S为滑坡土粒子区域，W为水流粒子区域，B为边界粒子区域，R为岩石或者建筑物粒子区域，p

步骤2-4，岩石或建筑物的接触力计算；

对于岩石及建筑物，与其接触的粒子类型包括：同类型的有岩石及建筑物粒子，异类型的有滑坡土粒子、地形边界粒子、水流粒子，滑坡冲击岩石及建筑物过程中，岩石或建筑物受到来自滑坡土粒子、水流粒子、地形边界粒子以及同类型岩石或建筑物粒子的接触作用。

设岩石或各个建筑物a由一系列相应的粒子k组成，基于如下公式计算粒子速度及接触力：

式中，组成岩石或建筑物a具有质量M

如图3所示，图中显示了滑坡土粒子与滑坡土粒子接触①情况、建筑物粒子与建筑物粒子之间接触②情况、滑坡土粒子与建筑物粒子接触③情况、岩石粒子与水流粒子接触④情况、滑坡土粒子与水流粒子接触⑤情况、滑坡土粒子与边界粒子接触⑥情况、水流粒子与边界粒子之间接触⑦情况、建筑物粒子与边界粒子之间接触⑧情况。对于水流粒子和边界粒子，其同类型粒子计算方式与滑坡土粒子相似，即接触①，在图中不再重复标记。

步骤2-5，滑坡土及岩石建筑物侵蚀启动判断；

滑坡触发后失稳滑坡土会对沿途的滑坡土、岩石及建筑物造成冲击，若达到侵蚀启动条件，则会带走滑坡土、岩石、建筑物一起移动；侵蚀包括对沿途土体的侵蚀和对沿途岩石及建筑物的侵蚀；对待判断的滑坡土粒子，通过临界剪切强度确定是否触发侵蚀启动，若触发侵蚀启动，则该滑坡土粒子被侵蚀，成为被侵蚀的滑坡土粒子；对于待判断的岩石或建筑物粒子，通过冲击碰撞力确定是否触发侵蚀启动，若触发侵蚀启动，则该岩石或建筑物粒子被侵蚀，成为被侵蚀的岩石或建筑物粒子；被侵蚀的滑坡土粒子作为失稳滑坡土粒子在下一个时间步计算过程中，采用步骤2-3的公式进行计算模拟。

(1)失稳后滑坡土对下游滑坡土粒子的剪切力来判断被侵蚀的滑坡土，当沿途的滑坡土粒子所受到的冲击剪切力大于临界剪切强度值时，当前滑坡土粒子被侵蚀铲刮带走，成为被侵蚀的滑坡土粒子，否则当前滑坡土粒子仍处于静止状态，采用Drucker-Prager屈服准则判断土体是否屈服：-αp+κ<2μγ，式中参数α和κ是与粘聚力和内摩擦角有关的参数，γ为剪切应变率。

如图4所示，失稳后滑坡土向下游冲击，侵蚀沿途土体(沿途滑坡土)，可将沿途土体分为：屈服液化土体、临界屈服区域土体、沿途未屈服土体。屈服液化土体此时已经失去强度，土体满足Drucker-Prager屈服准则，其所受剪切应力大于屈服应力，后续将会被侵蚀带走。临界屈服区域土体此时处于临界状态，其所受剪切应力等于屈服应力。沿途未屈服土体表示土体仍然静止，其所受剪切应力小于屈服应力。

(2)对于岩石或建筑物的侵蚀启动，主要通过冲击碰撞力计算。

如图5所示，当滑失稳滑坡土冲击碰撞力小于岩石(块石)或建筑物的整体启动力时，则岩石或建筑物不发生整体侵蚀启动；当失稳滑坡土冲击碰撞力大于岩石或建筑物整体启动力时，则岩石或建筑物整体侵蚀启动发生，开始随着失稳滑坡土一起运动。

如图6所示，此处冲击碰撞力通过冲击力响应计算算法得到：首先确定需要计算的目标结构(岩石或建筑物)的组成粒子具体信息，包括单个粒子k的质量、组成该目标结构的粒子的数量N；然后，通过公式2-3-2对影响域内位于粒子k周围的其余粒子进行插值计算，求得组成目标结构的所有粒子的加速度，其余粒子类型包括：滑坡土粒子、水流粒子、地形边界粒子、岩石粒子、建筑物粒子；对于每个粒子k，均根据x、y、z三个方向求得的分力进行

如图7所示，在失稳滑坡土冲击碰撞力小于岩石或建筑物整体启动力的情况下，岩石或建筑物还可能发生局部侵蚀，发生局部侵蚀时，被冲击碰撞的岩石或建筑物会发生碎裂解体效应，碎裂后的岩石粒子或建筑物粒子会被失稳滑坡土侵蚀带走。滑坡触发后，向下游运动，侵蚀沿途土体及岩石，当岩石被侵蚀带走后，运动过程中，岩石会与沿途岩石及建筑物相互碰撞，此时被碰撞的岩石或建筑物会发生碎裂解体效应，碎裂后的岩石被滑坡侵蚀带走，进一步增加了失稳滑坡体的体积。同时，岩石的冲击也会造成沿途建筑物发生损伤甚至断裂，不考虑此裂解效应则会低估岩石对建筑物的冲击力以及建筑物失事风险。

①采用如下方法判断岩石或建筑物粒子是否发生局部侵蚀：

分别计算粒子的拉伸应力和剪切应力，粒子的拉伸应力计算公式为：

剪切应力计算公式为：

式中，F

②计算得到粒子的拉伸应力和剪切应力后，进行断裂损伤破坏判断，满足断裂破坏准则，则认为该粒子发生了局部侵蚀导致的损伤破坏并成为损伤粒子。

本实施例中，断裂破坏准则采用拉伸和剪切两种破坏模型进行判断，公式如下：

σ

τ

式中，σ

③破坏后粒子应力更新及损伤统计：

统计损伤粒子的数量，用于损伤进行实时量化，并将损伤粒子的抗拉强度σ

步骤2-6，基于步骤2-1～步骤2-5计算得到当前时间步内各粒子特性信息和侵蚀情况数据，然后根据当前时间步计算得到的数据，确定下一个时间步计算用的初始参数值；重复前述过程，直至达到预设时间步数。

步骤3，根据步骤2构建的滑坡触发及传播灾害链全过程模拟框架，预测潜在滑坡土触发、演进、侵蚀沿途土体及建筑物、冲击下游建筑物全过程，预演滑坡触发及传播过程，量化滑坡灾害。

通过滑坡触发及传播灾害链全过程模拟框架每一计算步循环，最终预测强震或暴雨作用下潜在滑坡体触发、演进、侵蚀沿途土体及建筑物、冲击下游建筑物等全过程，动态计算滑坡体运动速度、到达时间、滑坡体各时刻运动形态及最终堆积形态、滑坡体对坝属建筑物的时变冲击力等，量化滑坡灾害链增容效应和灾害放大效应，可预演滑坡触发及传播过程。

根据滑坡土粒子、岩石粒子的特征信息和岩石及建筑物的侵蚀情况量化滑坡灾害链增容效应和灾害放大效应。增容效应量化包括：量化被侵蚀的滑坡土粒子和发生整体或局部侵蚀的岩石与建筑物的相应组成粒子使失稳滑坡土的体积增加的部分。对于多层滑坡土体，上下层粒子之间存在侵蚀作用，上一层滑坡触发后被侵蚀的下一层滑坡土粒子也算在上一层滑坡土体中作为增加的一部分。灾害放大效应量化包括：量化发生侵蚀后滑坡的冲击速度和对下游建筑物的冲击力的变化情况。

下面对本发明方法进行合理性验证，主要包括冲击力试验对比、侵蚀过程表征、土-岩-结构冲击过程表征，最后给出一个滑坡-堵河-冲击大坝附属建筑物工程案例进行佐证，包括滑坡体运动过程、堆积范围、大坝附属建筑物冲击力响应等。

滑坡触发前采用固体本构表征其小变形特性，采用强度折减法计算滑坡变形，计算滑坡触发临界强度折减参数，计算此时滑坡塑性区分布，并以此时刻作为滑坡触发时间点。如图8给出了以塑性区贯通和位移突变作为综合判断滑坡触发的条件。

图9表明折减系数在1.10和1.15之间塑性区基本贯通，滑坡体位移在1.05和1.10之间发生突变，说明此时滑坡体基本滑动，取两者最小值1.05作为最终触发时刻。图10表明本发明数值模拟得到的冲击力结果与试验实测冲击力数值较为接近，一致性较好，说明了本发明方法的合理性和可靠性。图11表明本发明可较好地表征土体与土体之间的侵蚀行为。图12～13的滑坡案例为某次水库滑坡，该滑坡在强降雨作用下触发。图12中第一幅图为滑坡初始示意图，模型包括：地形、岩石(岩石1、岩石2)、启闭房、滑坡土体(此案例认为水流与土体混合，不对水流进行单独建模)，第二至四幅图为模拟的动态滑坡过程不同阶段。图13为不同结构的受力情况。图12～13表明，本方法可较好地表征滑坡运动过程中土体、岩石及建筑物相互冲击碰撞动态过程及受力情况。由图12～13可知，本发明可较好地表征滑坡运动过程中滑坡土体、岩石及建筑物相互冲击碰撞动态过程，并且本发明在实际滑坡工程风险防控中，可根据预先计算得到的滑坡堆积形态判断建筑物是否会被掩埋，如发生掩埋，可根据进一步根据滑坡冲击范围以及建筑物受力情况，根据建筑物的抗冲性能及所受冲击力判断建筑物是否会在此类滑坡冲击荷载下被冲走或者损毁。因此，设计人员以及大坝管理人员可根据此算法结果评估滑坡触发后其潜在影响范围内的建筑物的安全状态。

最后，给出工程案例，证明本发明在实际工程应用的可行性，首先建立滑坡体-地形-河道-大坝附属建筑物三维精细化模型，如图14所示。采用本算法计算得到的滑坡体运动过程、堆积范围、建筑物冲击力响应等如图15～17所示。图15表明，本发明能较好地模拟滑坡体堵河冲击大坝附属建筑物过程。图16及图17给出了滑坡体最终堆积形态及建筑物所受冲击力的预测结果，表明本发明可用于滑坡堵河冲击坝属建筑物全过程模拟和危险性预测。

<实施例二>

本实施例二提供能够自动实现上述本发明方法的滑坡堵河冲击大坝区域建筑物过程的模拟预测系统，该系统包括建模部、模拟框架构建部、模拟预测部、输入显示部以及控制部。

建模部用于执行上文步骤1所描述的内容，根据待模拟的情况，通过建模软件建立相应的滑坡土、地形边界、大坝区域建筑物、岩石及水流三维实体模型，并离散为粒子模型，分别对滑坡土粒子、地形边界粒子、建筑物粒子、岩石粒子及水流粒子这五类粒子赋予初始参数值和用于区分各类型粒子的粒子属性信息。模拟框架构建部用于执行上文步骤2所描述的内容，在建模部的基础上，构建滑坡触发及传播灾害链全过程模拟框架。模拟框架构建部包括滑坡前模拟单元、触发确定单元、滑坡后模拟单元、接触力计算单元、侵蚀判断单元、迭代单元。滑坡前模拟单元用于执行上文步骤2-1所描述的内容，进行滑坡触发前的模拟。触发确定单元用于执行上文步骤2-2所描述的内容，根据滑坡触发临界转化的判断条件确定是否触发滑坡。滑坡后模拟单元用于执行上文步骤2-3所描述的内容，进行滑坡触发后冲击沿途岩石及建筑物过程的模拟。接触力计算单元用于执行上文步骤2-4所描述的内容，计算岩石或建筑物的接触力。侵蚀判断单元用于执行上文步骤2-5所描述的内容，进行滑坡土及岩石建筑物侵蚀启动判断。迭代单元用于执行上文步骤2-6所描述的内容，基于滑坡前模拟单元、触发确定单元、滑坡后模拟单元、接触力计算单元、侵蚀判断单元计算得到当前时间步内各粒子特性信息和侵蚀情况数据，然后根据当前时间步计算得到的数据，确定下一个时间步计算用的初始参数值；重复前述过程，直至达到预设时间步数。模拟预测部用于执行上文步骤3所描述的内容，采用滑坡触发及传播灾害链全过程模拟框架模拟预测潜在滑坡土触发、演进、侵蚀沿途土体及建筑物、冲击下游建筑物全过程，预演滑坡触发及传播过程，并量化滑坡灾害。输入显示部让操作员输入操作指令，并能够根据操作指令对各个部的输入、输出和中间处理信息进行相应显示。控制部与建模部、模拟框架构建部、模拟预测部、输入显示部均通信相连，控制它们的运行。

以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的滑坡堵河冲击大坝区域建筑物过程的模拟预测方法及系统并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。