CN109359391A - 一种基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法 - Google Patents

Abstract

本发明属于地质灾害机理研究领域，公开了一种基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法，利用离散元数值方法对覆盖型岩溶塌陷机理及灾变演化过程进行分析，从微观角度上揭露岩溶塌陷全过程，包括：裂纹形成、颗粒剥落、土洞形成、土洞扩张、覆盖层塌陷等；分析塌陷过程中覆盖层土体颗粒的位移、裂纹的发展趋势、系统不平衡力的变化情况等，以真实反映出外界营力作用下灾变主体(覆盖层土体)的多场变化特征。本发明运用离散元数值模拟方法得出覆盖型岩溶塌陷的临界土洞平衡高度和土洞开始形成的临界地下水流速度，并与理论结果进行对比，其模拟结果可为岩溶塌陷地质灾害的勘察处治、防灾减灾工作提供可靠依据。

Description

一种基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法

技术领域

本发明属于地质灾害机理研究领域，尤其涉及一种基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法。

背景技术

目前，业内常用的现有技术是这样的：

覆盖型岩溶塌陷灾害，因其具有隐蔽性、突发性、不确定性、机理复杂性等特点，导致预防、治理岩溶塌陷灾害十分困难。为了减少岩溶塌陷灾害的发生，对覆盖型岩溶塌陷进行分析十分必要。

岩溶塌陷数值模拟方法主要有限单元法、有限差分法、积分方程法和分界单元法。其中采用最多的是有限差分法——FLAC3D。FLAC3D是依据于连续介质力学理论数值模拟方法，往往是把分析体当成一个整体来分析，在分析土洞周围土体应力、应变问题上具有较好的模拟效果，但灾害体为非连续介质时，有限差分法则无法模拟。岩溶塌陷过程中，土体大位移、颗粒迁移、裂纹发展及土体崩落等现象及问题是有限差分数值模拟无法解决的，且无法得到覆盖层土体从裂纹发展——土洞形成——覆盖层塌陷等一系列离散过程和塌陷的临界阈值。现阶段，前人已用离散单元数值模拟方法在山体滑坡、工程爆破等领域进行过诸多探索；但鲜有文献运用离散元方法分析讨论岩溶塌陷的演化微观过程及灾变阀值等问题。

综上所述，现有技术存在的问题是：

对于覆盖型岩溶塌陷问题而言，之前的分析及模拟方法(如有限单元法、有限差分法等)几乎全部是基于连续介质理论的数值模拟方法。该类方法以“连续介质力学”为理论基础，仅仅能够模拟覆盖型岩溶区盖层土体(土洞)的整体大变形、挠曲和塑性流动等力学行为；但土洞的盖层土体由破碎的固体颗粒组成，土的宏观变形主要不是由于颗粒本身变形，而是由于颗粒间位置即排列方式的变化造成的，以上传统的数值模拟方法无法体现土颗粒的动态迁移、潜蚀以及崩解过程。

现有技术中，对粒状集合体的破裂和裂纹发展问题、以及颗粒的流动(大位移)问题的微观力学程序颗粒流PFC(Particle Flow Code)，还没有建立覆盖型岩溶塌陷数值模型，无法模拟土洞的扩张、颗粒的剥落、地表的塌陷，从而导致不能从微观角度上分析岩溶塌陷过程中土体颗粒的位移、裂纹的扩展等问题；不能定量化地揭示覆盖型岩溶塌陷的影响因素(地质条件、外界营力)，求解初始土洞形成、地表塌陷的临界阀值，也没有与理论值进行对比分析，无法为治理岩溶塌陷灾害提供理论依据。

覆盖型岩溶系统内的土洞，在其破坏过程中，不仅存在挠曲、塑性流动等基于连续介质假设的大变形力学行为；同时也存在土洞下表面土颗粒潜蚀、迁移、崩落等基于非连续介质假设的力学行为，但传统的基于“连续介质力学理论”数值模拟方法无法对这些力学行为进行刻画。

解决上述技术问题的意义：

解决上述技术问题后，带来的意义为：基于“非连续介质力学理论”，利用离散单元法进行模拟，从微观角度通过颗粒结构的微观参数研究来分析材料的宏观力学行为，不仅可以模拟塌陷体连续状态下及非连续状态下的变形特性，更可以揭示塌陷体由连续体到非连续体的渐进破坏过程，对丰富岩溶塌陷的理论研究体系有着重要理论意义，同时对于覆盖型岩溶区塌陷地质灾害的勘察处治、防灾减灾也有着重要的社会价值。

发明内容

针对现有技术存在的问题，本发明提供了一种基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法。对于覆盖型岩溶塌陷问题而言，之前的分析几乎全部是基于“连续介质理论”的数值模拟方法；而本发明的方法，是利用基于“非连续介质”的离散单元法理论，利用PFC软件进行岩溶塌陷的过程模拟。

本发明以典型覆盖型岩溶塌陷为分析对象，在基本明确致塌要素的基础上，利用离散元数值(PFC)方法对覆盖型岩溶塌陷机理进行分析，从微观角度展示岩溶塌陷全过程，包括：裂纹形成——颗粒剥落——土洞形成——土洞扩张——覆盖层塌陷等；分析塌陷过程中覆盖层土体颗粒的位移、裂纹的发展趋势、系统不平衡力的变化情况等，以真实反映出外界营力作用下灾变主体(覆盖层土体)的多场变化特征。本发明运用离散元数值模拟方法得出覆盖型岩溶塌陷的临界土洞平衡高度和土洞开始形成的临界地下水流速度，并与理论结果进行对比，其模拟结果为覆盖型岩溶区塌陷灾害的勘察处治及防灾减灾工作提供依据。

本发明是这样实现的，一种基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法，所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法包括：

在覆盖型岩溶区典型塌陷的地质结构概化模型基础上，利用颗粒流理论数值模拟方法对土洞的形成过程、覆盖型岩溶的塌陷过程进行模拟，获取初始土洞形成的临界地下水流速度和覆盖型岩溶塌陷的土洞临界平衡高度。

进一步，所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法具体包括：

第一步，运用PFC模拟渗流模型，结合土体颗粒位移曲线、裂纹增长曲线、不平衡力变化曲线判断是否发生渗流破坏；

第二步，利用PFC数值软件模拟土洞初始发育阶段，得到土洞形成的临界地下水流速度，并与理论值作对比；

第三步，运用PFC软件对覆盖型岩溶塌陷发展过程进行模拟：

对覆盖型岩溶塌陷过程中的应力变化情况进行模拟；

对覆盖型岩溶塌陷进行数值模拟，将塌陷过程分为四个阶段：第一阶段，土洞向上扩展，形成第一级土洞；第二阶段，土洞向两侧扩展，形成第二级土洞；第三阶段，土洞继续缓慢发展，覆盖层中多次发生小型坍塌，第三级土洞形成；第四阶段，土洞规模不在扩大，裂纹数目急速上升且形成垂连通面，地表发生塌陷；

覆盖型岩溶塌陷过程中的微观现象模拟；

利用PFC模拟软件计算土洞最大临界高度，与理论值对比。

进一步，第一步包括：应力平衡→应力集中→裂纹形成→裂纹形成连通面→颗粒剥落→应力平衡。

进一步，覆盖型岩溶塌陷过程中的微观现象模拟中，包括裂纹扩张、土体颗粒位移、土洞发育、地表塌陷过程的模拟。

进一步，离散元的微观耦合力学方程：

PFC中流固耦合计算的流动方程、压力方程和求解条件如下：

流动方程流体管道相当于一个平行通道，长度为L’、孔径为a，在垂直平面方向上为单位厚度，管道内的流速(单位时间内的体积)为：

K^‘—传导系数(cm·s^-1)；L’—管道的长度(mm)；P₂-P₁—相邻域的压力差；

压力方程

周围管道流入每个域的流量和为∑q，在单位时间步长Δt下，流体压力增量Δp(流入为正)为：

式中：K_f——流体的体积模数(kPa)；V_d——域的表观体积(mm³)；

应用显示求解方法，将流量方程应用于所有的管道，并将压力方程应用于所有的域之间交替求解；假设某个域内存在扰动压力ΔP_p，由于扰动流域里的流量可以从式(1)计算得：

由水流流入引起的响应压力变化

式中：N——域所黏结的管道数；R——域周围颗粒的平均半径(mm)；

当两者相等时可求出临界时间步长为：

本发明的另一目的在于一种实现所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法的计算机程序。

本发明的另一目的在于一种实现任意一项所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法的信息数据处理终端。

本发明的另一目的在于一种计算机可读存储介质，包括指令，当其在计算机上运行时，使得计算机执行所述的基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法。

本发明的另一目的在于一种实现所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法的基于离散单元法的覆盖型岩溶塌陷灾变演化模拟控制系统。

本发明的另一目的在于一种搭载所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟控制系统的计算机。

综上所述，本发明的优点及积极效果为：

本发明通过地质勘查对岩溶塌陷的特点、覆盖层性质、形成原因、致塌机理进行分析，归纳总结出覆盖型岩溶区典型塌陷的地质结构概化模型。在此基础上，利用基于颗粒流理论的离散单元(PFC)数值模拟方法对土洞的形成过程、覆盖型岩溶的塌陷过程(如图16覆盖型岩溶区土洞塌陷演化过程示意图，塌陷过程A、B、C、D)进行了模拟，并获取了初始土洞形成的临界地下水流速度和覆盖型岩溶塌陷的土洞临界平衡高度。

本发明的主要技术效果还在于：

现有的以―连续介质力学”为理论基础的数值模拟方法(例如：(FLAC)快速拉格朗日差分法)，仅仅能够模拟覆盖型岩溶区盖层土体(土洞)的整体大变形、挠曲和塑性流动等力学行为；

但土洞的盖层土体由破碎的固体颗粒组成，土的宏观变形主要不是由于颗粒本身变形，而是由于颗粒间位置即排列方式的变化造成的，传统的基于―连续介质力学”理论的数值模拟方法(例如：(FLAC)快速拉格朗日差分法)无法体现覆盖层土颗粒的动态迁移、潜蚀以及崩解过程与现象。

而本发明所述的基于―非连续介质力学理论”的离散单元法，克服了传统连续介质力学模型的宏观连续性假设，可以从微观角度对土的工程特性进行数值模拟，通过颗粒结构的微观参数研究来分析材料的宏观力学行为，体现塌陷过程中土颗粒潜蚀、迁移、崩落等基于“非连续介质假设”的力学行为,如图17基―非连续介质力学”理论的数值模拟则可以很好的体现土体颗粒迁移、崩落等塌陷过程中的真实现象。从而更加真实的模拟塌陷体由连续体到非连续体的渐进破坏过程。如图18基于―非连续介质力学”理论的离散单元数值模拟方法技术效果图。

附图说明

图1是本发明实施例提供的基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法流程图。

图2是本发明实施例提供的PFC2D计算过程图。

图3是本发明实施例提供的颗粒流接触关系示意图。

图4是本发明实施例提供的巴西劈裂实验示意图。

图5是本发明实施例提供的初始弹性模量图。

图6是本发明实施例提供的体积应变-轴向应变图。

图7是本发明实施例提供的应变-应力曲线图。

图8是本发明实施例提供的三种围压条件下偏应力与约束应力随轴向应变变化曲线图。

图9是本发明实施例提供的三种围压下的应力-应变曲线图。

图10是本发明实施例提供的岩溶塌陷地质Y-Z剖面图。

图11是本发明实施例提供的溶蚀沟槽两侧颗粒1位移图。

图12是本发明实施例提供的溶蚀沟槽两侧颗粒2位移图。

图13是本发明实施例提供的临界地下水流速度-裂纹数目变化曲线图。

图14是本发明实施例提供的临界地下水流速度—模型系统不平衡力变化图.

图15是本发明实施例提供的流量-运行时间曲线图。

图16是本发明实施例提供的覆盖型岩溶区土洞塌陷演化过程示意图。

图17是本发明实施例提供的基于“非连续介质力学”理论的数值模拟则可以很好的体现土体颗粒迁移、崩落等塌陷过程中的真实现象图。

图18是本发明实施例提供的基于―非连续介质力学”理论的离散单元数值模拟方法技术效果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

本发明利用适用于分析粒状集合体的破裂和裂纹发展问题、以及颗粒的流动(大位移)问题的微观力学程序颗粒流PFC(Particle Flow Code)建立覆盖型岩溶塌陷数值模型，模拟土洞的扩张、颗粒的剥落、地表的塌陷等，并从微观角度上分析岩溶塌陷过程中土体颗粒的位移、裂纹的扩展等问题。利用PFC软件定量化地揭示覆盖型岩溶塌陷的影响因素(地质条件、外界营力)，求解初始土洞形成、地表塌陷的临界阀值，并与理论值进行对比分析。其结论将提高岩溶塌陷防治水平，为治理岩溶塌陷灾害提供了理论依据。

图1，本发明实施例提供的基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法包括：

S101:运用PFC模拟渗流模型，结合土体颗粒位移曲线、裂纹增长曲线、不平衡力变化曲线判断是否发生渗流破坏；

S102:利用PFC数值软件模拟土洞初始发育阶段，得到土洞形成的临界地下水流速度，并与理论值作对比；

S103:运用PFC软件对覆盖型岩溶塌陷发展过程进行模拟：

对覆盖型岩溶塌陷过程中的应力变化情况进行模拟；

对覆盖型岩溶塌陷进行数值模拟，将塌陷过程分为四个阶段：第一阶段，土洞向上扩展，形成第一级土洞；第二阶段，土洞向两侧扩展，形成第二级土洞；第三阶段，土洞发展较为缓慢，覆盖层中多次发生小型坍塌，第三级土洞形成；第四阶段，土洞规模不在扩大，裂纹数目急速上升且形成垂连通面，地表发生塌陷；

覆盖型岩溶塌陷过程中的微观现象模拟；

利用PFC模拟软件计算土洞最大临界高度，与理论值对比。

一、下面结合覆盖型岩溶(土洞)塌陷的形成条件对本发明作进一步描述。

覆盖型岩溶地区，塌陷一般是指由于土层中土洞扩展，导致地面陷落而产生地表变形和破坏。根据已有的覆盖型岩溶区塌陷案例来看，岩溶塌陷的形成原因实质上是土洞(或岩洞)的抗塌力小于致塌力的结果。覆盖型岩溶塌陷过程可划分为三个阶段：第一阶段：碳酸盐岩地层在地下水的作用下形成形成溶蚀沟槽；第二阶段：溶蚀沟槽处的土体颗粒剥落、土体内部发生坍塌并发育成土洞；第三阶段：土洞规模逐渐扩大，导致地表失稳形成岩溶塌陷。

土洞的形成过程很缓慢，需要具备一定条件下方可形成，土洞的形成所必须具有的条件：①基岩为可溶蚀性岩体，具有溶蚀沟槽，且溶蚀沟槽开口向上——空间条件；②覆盖层具有一定厚度——物质条件；③地下水的变化——动力条件。三个条件缺一不可。

二、下面从颗粒流理论介绍及土体微观参数标定对本发明作进一步描述。

1、颗粒流理论概述

二维颗粒流程序(Particle Follow Code PFC2D)其理论基础是Cundall提出的离散单元法，一种非连续力学分析方法，基于此类方法来模拟颗粒之间的运动关系及其受力情况，此方法在岩石力学、构造地质、机械工程等中有广泛应用。

1.1离散元的理论概述

在岩土体的数值模拟中，土体颗粒的本构关系获取十分困难，但PFC不需要给材料参数定义宏观本构关系和对应的参数，可直接定义颗粒和黏结的几何和力学参数。其中，颗粒可以代表材料中的单个颗粒，例如土体颗粒，也可以代表黏结而成的固体材料，例如土层、岩块等。当黏结以渐进的方式破坏时，固体材料可以破裂。颗粒间的接触可选择以下模型：(1)线性弹簧或Hertz-Mindlin法则；(2)库仑滑块；(3)点接触、用平行的弹簧黏结。接触模型决定了宏观物质的基本力学特性、破坏形式等。

相对于有限元数值模拟来说PFC不但可以表征宏观物质的物理特性，而且可以反映其他方法无法实现的微观特性。对于离散软件而言，其优点：第一，它有潜在的高效率；第二，可以模拟颗粒大位移问题，对位移无限制；第三，模型是由颗粒黏结而成，模型可形成断裂，模拟材料破坏等，此特点是其他离散软件所不具有的。

1.2颗粒流分析塌陷问题合理性

覆盖型岩溶塌陷区其特点为：覆盖层岩性多为松散的粘土、砂砾土，且覆盖层厚度较小，基岩发育溶蚀沟槽、且地下水水位在覆盖层与基岩溶蚀沟槽交界处大幅度变化。分析区覆盖型岩溶塌陷其形成原因多为：降雨入渗、地下水水位下降引起的。对于覆盖型岩溶塌陷问题而言，塌陷过程中存在土体受力、变形、颗粒剥落等问题，覆盖层土体在塌陷中发生如下变化：土体颗粒在外力的作用下发生大变形——土体颗粒位移增大——覆盖层产生裂纹——裂纹发展形成连通面——颗粒沿着连通面剥落——覆盖层塌陷，在此过程中存在着大变形和大位移问题。现如今多用有限差分法FLAC3D对覆盖型岩溶塌陷问题进行分析，往往是把分析体当成一个整体来分析，对于灾害体为非连续介质时，有限差分法存在较大局限性。

颗粒流数值模拟方法(PFC)对覆盖型岩溶塌陷问题的分析中，PFC软件自带流体计算功能，能很好的模拟由降雨垂直入渗和地下水水位变化引起的潜蚀作用，并能进行流固耦合计算。颗粒流数值模拟方法 (PFC)适用于模拟固体力学大变形问题，可以得到土体变形——颗粒产生大位移——颗粒剥落——土体坍塌的全过程，且能较好的刻画地质构造、基岩溶蚀沟槽等。因PFC可对模型颗粒直接进行参数赋值，并根据不同位置、不同情况赋予不同数值，在覆盖型岩溶塌陷问题中可根据覆盖层和基岩的力学性质对模型参数给与不同数值。且PFC能显示裂纹发展过程和揭露土体塌陷过程，还能全面地给出从微观变化到宏观变化的各种参数响应信息，能运用微观变化解释宏观问题，从本质上分析塌陷问题。综上所述利用颗粒流方法能更有效的解决覆盖型岩溶塌陷问题。

2、颗粒流程序计算模型

2.1离散计算过程

离散颗粒运动遵循于力-位移定律和牛顿第二定律，通过力-位移定律，改变接触部分的接触力；PFC 的计算方法为时步迭代法，通过不断迭代求出正解，每一步计算更新颗粒单元、墙的位置和接触关系。在计算过程中运用牛顿第二定律，改变颗粒与墙体的位置，重新得出颗粒之间接触关系。如图2所示。

2.2力-位移方程

在离散元模拟过程中，颗粒和颗粒接触会产生力，运用力-位移定律将广义内力与接触处的相对运动联系起来。内力和力矩作用在两个颗粒的接触位置。其接触分为颗粒—颗粒接触，颗粒—墙的接触，接触类型如图3所示，其中R为颗粒半径、x为型心坐标、U表示颗粒之间接触所产生的重叠量d为颗粒流型心之间的距离，n为接触单位法向量。

颗粒和颗粒、颗粒和墙体接触单位法向量：

球-球接触：

球-墙接触：

颗粒和颗粒、颗粒和墙体接触产生的重叠位移：

球-球接触：

球-墙接触：Uⁿ＝R^[b]–d (4.4)

颗粒和颗粒、颗粒和墙体接触点的中心坐标：

球-球接触：

球-墙接触：

接触力F_i可分为法向分量F_i ⁿ和切向力分量F_i ^s

在颗粒之间的接触形成初期，接触力切向分量为0，随着模型运行，颗粒间的切向相对位移产生的接触力切向分量，接触力的切向相分量进行累加，且在模型运行过程中，接触点和位移时刻变化，接触力切向分量可表示为：

切向接触力增量ΔF_i ^s表示为

式中：K^s为接触点的切向刚度，为切向位移量。

2.3运动定律

在模拟过程中，颗粒单元的运动是由其质心所受的不平衡力和不平衡力矩决定的，其运动可分为平移运动和旋转运动。对于平移运动可用颗粒速度和颗粒加速度来描述，其方程式为：

其中：F_i为不平衡力，m为颗粒质量，g为重力加速度。

对于不平衡力矩与旋转运动的关系，可用角速度w和角加速度w’来描述，颗粒旋转运动方程为：

其中：M_i为不平衡力矩和为单元角动量。

M＝I w^’＝(βMR²)w’ (4.12)

式中：I为颗粒的主惯性矩，w’为角加速度，β为因子，当数值模型为二维模型时，β为1/2。对平移和旋转运动方程加进行中心有限差分法可得：

将式子(4.14)、(4.15)带入式(4.11)、(4.13)，求解得出t±Δt/2的速度为：

求解得出颗粒的中心位置为：

3、颗粒流的接触模型

在PFC中，材料的本构特性是由颗粒的实体接触本构模型决定的，常用的本构模型为黏结模型，且黏结模型可分为接触黏结模型和平行黏结模型。黏结模型多用于模拟土体、岩体。

3.1接触黏结模型

接触黏结为颗粒与颗粒之间只有力学关系。一旦颗粒之间形成接触黏结，颗粒之间的黏结强度小于最大黏结强度时，颗粒之间不会发生相对滑动，不产生相对位移。相互接触的颗粒之间的重心距离大于两颗粒的半径和时，表示颗粒之间脱离接触，颗粒之间产生的力为张力。

在PFC中，由法向黏结强度和切向黏结强度共同决定颗粒的接触强度，其力学单位为Pa。法向接触强度决定颗粒的抗拉强度，法向接触强度决定颗粒的抗剪强度。若颗粒之间的法向力大于最大法向接触强度时，颗粒之间的接触黏结破坏，法向、切向接触强度为0，颗粒间的接触力为0；当颗粒之间的切向力大于最大切向接触强度时，接触同样失效，但没有超过摩擦极限时，接触力保持不变。

3.2平行黏结模型

平行黏结模型用来描述颗粒集合存在胶结物使材料具有粘聚性，这种黏结建立在颗粒接触点，材料的截面可以为圆形或者是椭圆形，例如岩石材料。

平行黏结模型可在颗粒黏结处产生力和力矩，力和弯矩可分解为法向分量和切向分量。颗粒之间的任何方向上接触力大于最大黏结强度时，平行黏结破坏。

平行黏结开始形成时，总接触力和总力矩为0，在模型运行中，随着接触处的位移和旋转量的增加导致总接触力和力矩的增加，力和力矩的增加量将在下一步中与当前值进行叠加，此叠加过程在计算过程中不断重复。

通过对学者多次分析可知：在分析土工问题中，所有的土体由不连续的固体颗粒组成，土体宏观变形为颗粒之间的位移变化。运用颗粒流数值方法时，模拟对象为离散型土体颗粒时，其本构模型为滑动接触模型，对于强度较大具有水泥加固特性的材料，选用平行黏结模型，对于具有一定初始黏结强度的土体，选接触黏结模型。本发明中运用PFC模拟覆盖型岩溶塌陷时，覆盖层土体具有一定的初始强度，构模型选用接触黏结模型。

3.3颗粒间黏结破坏过程

PFC能够模拟试样、岩土体中裂纹的产生、发展直至试样的破坏，因为颗粒之间的本构关系为接触黏结模型；颗粒之间的接触可以产生微小的裂纹，模拟实际情况中岩体的微型结构面的产生，颗粒之间的黏结逐渐断裂形成大的规模时，可模拟岩体大裂缝的形成直至岩体破坏失稳等情况。

由上述介绍可知，颗粒间的黏结强度分为法向黏结强度和切向黏结强度，如果颗粒间的拉力和切向力大于法向黏结强度或法向黏结强度，颗粒黏结破会，产生裂纹。且裂纹的数量和位置受黏结的微观参数影响，颗粒位置和颗粒大小也将影响裂纹的发育、数目等。

4、PFC2D水流耦合模型

4.1离散元的微观耦合原理

PFC中用圆形颗粒代替土体颗粒，引入“域”和“管道”的概念来表征和模拟真实的流体。“域”是一系列封闭的颗粒链，链上的每个链接都是一个黏结接触。“域”代替水压进行模拟，将压力等效成体力施加在颗粒上，在模型计算过程中，颗粒所受等效体力不断更新。假想固体中流体的通道在颗粒接触处，并与颗粒接触处相切，称该通道为“管道”。利用“管道”和“域”的黏结模拟颗粒与流体的耦合作用。

4.2离散元的微观耦合力学方程

PFC中流固耦合计算的流动方程、压力方程和求解条件如下：

(1)流动方程

流体管道相当于一个平行通道，长度为L’、孔径为a，在垂直平面方向上为单位厚度，管道内的流速(单位时间内的体积)为：

K‘—传导系数(cm·s^-1)；L’—管道的长度(mm)；P₂-P₁—相邻域的压力差。

(2)压力方程

周围管道流入每个域的流量和为∑q，在单位时间步长Δt下，流体压力增量Δp(流入为正)为：

式中：K_f——流体的体积模数(kPa)；V_d——域的表观体积(mm³)。

(3)求解方法

应用显示求解方法，将流量方程应用于所有的管道，并将压力方程应用于所有的域之间交替求解。假设某个域内存在扰动压力ΔP_p，由于扰动流域里的流量可以从式(1)计算得：

由水流流入引起的响应压力变化

式中：N——域所黏结的管道数；R——域周围颗粒的平均半径(mm)。

(4)求解条件

保证模型运行稳定的条件是水流入引起的压力变化必须小于扰动压力，当两者相等时可求出临界时间步长为：

5、PFC关于土体的微观参数

5.1土体常规力学参数的影响因素

颗粒流数值模拟是通过调整微观参数来表示不同宏观物质的物理特性。因为PFC软件的宏观力学参数和模型的微观参数相对应起来十分困难，且宏观参数的影响因素较多且复杂，使得模型的力学参数与实际情况相差甚大。在相同模型大小、颗粒级配、组装模式、黏结方式等前体下，运用PFC进行双轴实验、单轴实验、剪切实验等，通过调节数值模拟的微观参数使数值模拟实验和物理实验具有相同的破坏形式、裂纹分布、应力-应变曲线，来建立宏观参数和微观参数之间的关系，这个过程称为参数标定。为了保证岩土体的宏观物理特性与微观力学参数相互因，故在进行岩溶塌陷数值模型之前将对岩土体的宏观参数与微观参数进行标定。

对于岩土的分析，涉及的常规力学参数有：粘聚力、内摩擦角、初始弹性模量、抗拉强度、泊松比等。这些宏观参数与颗粒间的强度、接触刚度、颗粒间的摩擦系数和黏结强度等有关，但微观系数对不同的宏观参数的影响不同，通过PFC数值软件模拟双轴实验，拉伸试验，达西实验等，对起主要影响作用的微观参数进行标定，得到与宏观相对应的初始弹性变形模量、接触刚度、黏结强度等。

土体常规力学参数和微观力学参数之间的关系并非是线性的，几个微观参数共同影响一个宏观力学参数，通过前人的大量分析，对于宏观力学参数与微观参数的标定中，以下几个参数的决定性较大：

(1)颗粒间的接触模量法向刚度K_n、切向刚度K_s；

(2)法向刚度和切向刚度的比率e_mod；

(3)颗粒间的接触法向强度σ_c(Pa)；

(4)颗粒间的接触切向强度τ_c(Pa)；

(5)颗粒间未黏结的摩擦系数μ。

通过前人的大量分析得出：岩土体的泊松比的与法向刚度和切向刚度的比率e_mod、试样的几何形态有关；岩土体的初始弹性模量的与颗粒的刚度有关；岩体的抗拉强度由颗粒间的法向黏结强度决定；岩土体的峰值强度由切向黏结强度决定，岩土体的内摩擦角由颗粒间的摩擦系数决定。

渗透系数的影响因素渗透系数K，是重要的水文地质参数，渗透系数说明岩石的渗透性能。通过达西定律可知：

取长度L范围内所有截面的流量总和，在Δt时间内：

则渗透系数K的表达式为：

且由式子(4.18)-(4.22)可知在PFC软件计算中，流体流量与传导系数K’，管道直径a，域的表观体积 V_d、时间t、管道数量N有关。通过比较宏观渗透系数计算公式(4.25)与(4.18)-(4.22)可知，宏观渗透系数与传导系数K’，管道直径a，域的表观体积V_d有直接关联，与时间t、管道数量N有间接关联。白若虚等^[44]通过PFC计算讨论可得渗透系数与传导系数K’呈线性相关，与管道直径a的立方成正比，与域的表观体积V_d、与时间t无关。本发明中设a为一定常数——0.8mm，主要考虑K’对渗透系数的影响。

三、下面结合具体实验对本发明作进一步描述。

6.颗粒流数值模拟中的微观参数标定

本发明运用PFC数值方法模拟的土力学实验包括：双轴压缩实验、巴西劈裂试验、达西渗透实验。通过PFC模拟得出的粘聚力、内摩擦角、初始弹性模量、抗拉强度、泊松比、渗透系数等参数与实际值相比较，对PFC软件中相对应的微观参数进行标定，需要标定的微观参数有：法向刚度K_n、法向刚度和切向刚度的比率e_mod、接触法向强度σ_c、接触切向强度τ_c、摩擦系数μ。

6.1分析区土体常规参数

对典型岩溶塌陷地区覆盖层进行取样分析，覆盖层岩性为粘土，其常规力学性质如下表所示：

表1岩土宏观力学参数

6.2双轴压缩实验和巴西实验

(1)双轴实验

双轴模型为四面墙体单元把颗粒约束在墙体内，把颗粒压缩成长方形集合。模型所受的应力通过墙体所受合力和移动位移来控制。墙体和颗粒之间存在力学关系，而墙体之间不存在力学关系，墙体位移随颗粒模型的形变而变化。模型的破坏标准为颗粒的接触强度降低到某一数值时，模型运行结束。

在墙体对颗粒的加载过程中，通过history机制对模型的应力、位移等变量进行监测、记录。本发明中，记录的变量有应力、应变、裂纹数量，应力包括轴向、径向应力，应变包括轴向、径向、体积应变。双轴实验主要标定的是宏观参数为摩擦角和粘聚力，其对应的微观系数为弹性模量、泊松比、粘聚力。

(2)巴西实验

岩石的抗拉强度是岩土力学的重要参数，本发明采用巴西劈裂实验获得岩土的抗拉强度。图4为巴西劈裂实验示意图，墙体水平，法向加载，墙体和颗粒为点接触，对上下墙体给定一个恒速对颗粒进行加载法向应力，加载过程分多步实现，通过history记录加载过程中颗粒的受力情况，在巴西试验中，颗粒强度降低到某一数值(一般是指颗粒所受最大峰值强度的0.8倍)时，模型运行结束，此时判断模型已破坏，模型所受最大峰值强度为岩土的抗拉强度。

6.3双轴实验标定变形模量和泊松比

双轴实验模型大小选定为50mm×100mm，颗粒半径为0.1—1.8mm，密度为1960kg/m³，因为模型颗粒的力学加载是通过墙体实现的，墙体刚度的分析不可忽视，在PFC软件自带的接触模型中，默认为墙体的刚度和颗粒刚度相同。对于颗粒间的法向刚度K_n、切向刚度K_s和接触黏结模量的K_n’、K_s’取值相同。双轴试验中墙体设置为光滑平面。实验过程中，先给模型赋予固定围压，然后在法向上墙体给与速度，给模型施加压力。在数值模拟过程中，监测的变量如下：

(1)mv_wsd为记录的轴向偏应力σ_d；

(2)mv_wea为记录的体积应变ε_v＝ε_x+ε_y；

(3)mv_mea为记录的轴向应变ε_a。

在标定过程中，开始多次调节参数K_n、K_s的数值来确定泊松比ν和初始弹性模量E，K_n、K_s调试过程中，为减少其他微观参数的影响，其他参数赋予数值较大。在围压恒定的条件下，泊松比ν和初始弹性模量E由式(4.26)、式(4.27)确定。

由图5所示，x轴为轴向应变，y轴为轴向偏应力，由图可知轴向偏应力与轴向应变成正比，应力随着应变的增加而增加(摩擦系数、黏结强度此时赋予的较大)，其斜率为初始弹性模量E，通过式(4.26)可知：

由图6所示，x轴为轴向应变，y轴为体积应变，由图可知轴体积应变与轴向应变成正比，斜率为泊松比ν，通过式(4.27)可知：

此时模型中其他参数取值较大，通过多次调试K_n、K_s得到的初始弹性模量和泊松比与宏观参数较好的匹配，此时的颗粒的法向刚度K_n＝1×10⁷N/m，切向刚度Ks＝5×10⁶N/m，e_mod＝2。

6.4巴西劈裂标定抗拉强度

标定法向刚度K_n、切向刚度K_s后，通过巴西劈裂实验确定接触法向黏结强度的数值。巴西劈裂模型为直径为50mm的圆盘，墙体以一定速度向模型施加荷载，接触断裂的裂纹发展到一定地步时，模型破坏，墙体停止加载；

如图7应变-应力曲线图所示，模型破坏时，通过多次调试模型可得，PFC模型中最大峰值强度为4.8kPa，此时接触法向黏结强度为8×10⁴Pa。

6.5双轴实验标定粘聚力和内摩擦角：

对粘聚力和内摩擦角进行标定时，需做大量的数值模拟实验，标定变形模量、泊松比、抗拉强度，确定K_n、 K_s、σ_c后，在此基础上继续做双轴压缩试验，围压σ_c分别为30kPa、80kPa、130kPa，在不同围压

对式样进行法向加压直至加载强度到达最大破坏强度，在此过程中做出三种情况下的应变-应力变化曲线图(如图8所示)。

根据实验结果，绘制三种情况下的Mohr圆。其中，PFC软件中的围压σ_c＝σ₃(最小主应力)，σ_a为双轴压缩实验获得最大峰值强度，最大主应力σ₁＝σ_a+σ_c，绘制出应力圆后，做出三个应力圆的包络线，根据摩尔-库伦理论公式：τ＝σtanθ+C，求解出θ和C，三种围压下的应力-应变曲线如下图9所示，根据图9可知，式样的抗剪强度粘聚力为30Kpa，内摩擦角为26°，此时的微观参数摩擦系数为0.4，切向黏结强度为2×10⁴ KPa.

6.6达西实验(标定渗透系数)：

达西渗流模型尺寸为20mm×20mm。故本发明中颗粒半径为0.1-1.8mm，调节模型使颗粒均匀分布。颗粒之间采用黏结接触，对上、下边界的颗粒进行固定，标记为不排水边界。通过PFC生成“域”和“管道”所组成的网络所示，圆点为“域”，线段为“管道”，“域”中可储存水压，管道传送流体水压。模型中不排水的上下边界没有“管道”。

模型左侧水压固定为2×10³Pa，右侧水压固定为0，模型中水压从左向右逐渐蔓延。对所有模型的左、右侧分别预留1.5mm进行加压，则实际的渗流途径为17mm.根据达西定律可得，本发明中用水压进行计算，计算可得：J＝1146.47。

因无法测量某一过水断面的流量，故选取模型L＝10mm的范围测量总流量，流体流经总面积为200mm。根据PFC软件计算的：Q＝1.33×10⁵mm³，根据达西公式计算可得：K＝5.7×10^-7cm/s。管道直径a为常数 0.8mm，此时的水力传导系数K’＝5×10^-11cm/s。

参数标定结果

表2岩土宏观力学参数及标定值

表3微观参数取值表

本发明具有以下几点：

(1)颗粒流的理论背景、基本特点、用途、基本假设、颗粒的计算模型、接触模型和流固耦合模型，得出运用PFC模拟粘土时应选用接触黏结模型；

(2)土体常规力学参数的微观影响因素；

(3)说明用PFC软件建立覆盖型岩溶塌陷数值模型之前，对土体进行微观参数标定的必要性；

(4)利用PFC数值软件模拟双轴实验、劈裂实验、渗流实验，对模型宏观力学参数进行标定，得到的微观参数为岩溶塌陷数值模型提供参数基础。

四、下面结合岩溶塌陷颗粒流数值模拟对本发明作进一步描述。

本发明对典型岩溶塌陷区进行分析，在此基础上对塌陷区的地质条件进行概化。并利用PFC软件建立覆盖型岩溶(土洞)塌陷数值模型，从微观角度上剖析及揭示土洞的形成及演化过程，得到覆盖层土体颗粒的位移特征及土洞形成过程中裂纹扩展情况，求出覆盖型岩溶区土洞形成的临界地下水流速和临界土洞高度。其模拟结果丰富了覆盖型岩溶塌陷理论分析体系并为预防覆盖型岩溶塌陷提供理论依据。

7.覆盖型岩溶塌陷数值模型的建立

在覆盖型岩溶区，伴随着地下水的持续潜蚀作用，土洞开始形成并逐步向上扩展，致使洞顶上覆盖层土体厚度逐步减薄，拱顶薄到不能支持上部土的重量时，发生突然塌落，形成地表塌陷。

7.1地质条件概化

典型覆盖型岩溶地质概化模型表层为第四系覆盖层，盖层以粘性/粉质粘土土为主，第四系与灰岩交界面起伏变化大，常形成溶沟和石芽，在个别溶沟处存在土洞。特别是在石炭二叠系灰岩中溶洞成层性明显，往往水平延伸较长，尤其在地下水排泄带附近，钻孔多揭露串珠状溶洞。受碳酸盐岩中节理裂隙、断层等构造结构面影响，在灰岩的表层风化带，垂直岩溶发育，在地下水的长期溶蚀作用下，从而形成溶蚀沟槽和石芽等岩溶形态。若灰岩纯度相对较高，且呈中厚层状，则沿岩层理层面、软弱夹层、碳酸盐岩与非碳酸盐岩交界面等原生结构面易形成层状溶蚀，钻孔揭露为串珠状溶洞。而塌陷区覆盖层厚度多小于10m，厚度小于5m的地段塌陷点密集，且岩性多为粘土。根据以上分析对分析区典型岩溶塌陷模型刻画如下： (1)塌陷模型的大小为10m×6m，覆盖层厚度为5m，覆盖层以下为基岩，本发明中基岩厚度取为1m； (2)覆盖层岩性为粘土，基岩为灰岩；

(3)分析区岩溶裂隙及沟槽较为发育，土洞以下为的溶蚀裂隙/沟槽开口宽度刻画为1.2m；

(4)初始土洞刻画为高为H，弦为D的圆弧(H，D在模拟中随着模拟工况不同而变化)。岩溶塌陷地质模型剖面如图10所示.

7.2数值模型假定

模拟对象为典型覆盖型岩溶，为分析覆盖型岩溶塌陷过程中的临界状态，运用PFC模拟覆盖型岩溶塌陷，并得到塌陷过程、临界阈值。模型建模中做如下假设：

(1)对典型塌陷区的覆盖层土体进行取样分析，覆盖层岩性为粘土、粉质粘或含角砾粉质粘土，故本发明中所有数值模型的覆盖层岩性为粘土。

(2)岩溶的基岩(灰岩)为均质各向同性的地下空间半无限体，基岩为不可溶蚀、不可变形、强度较大的岩体(不参与岩溶塌陷的计算)。

(3)数值计算过程中，地下水渗流作用促使覆盖层土体内水力裂隙发育，为岩溶塌陷的主要动力来源且在PFC模拟塌陷过程中不考虑溶洞内部的真空吸蚀作用。

(4)在模型模拟、计算过程中，参与计算的力主要为：覆盖层土体的自重应力、孔隙水压力和渗透力。

(5)运用PFC数值模拟软件对岩溶裂隙及沟槽进行刻画，刻画为一定深度、一定大小、一定形状的沟槽，且模拟过程中岩溶沟槽不发生任何改变。

7.3岩溶塌陷数值模型建立

本发明岩溶塌陷模型以野外岩溶塌陷灾害点为实例，取Y—Z方向上的剖面为分析对象，运用PTC数值软件建立数值模型，对覆盖型岩溶塌陷形成原因进行分析。地质模型概化为10m×6m，因模型过大，计算量庞大，故对岩溶塌陷数值模型进行等比例缩尺，数值模型尺寸设定为100mm×60mm。利用PFC生成覆盖型岩溶的过程如下：

(1)颗粒的生成：颗粒生成在墙体之内，墙对颗粒进行束缚保持一定的形状，并对与墙体相接触的颗粒进行固定，作为模型边界。

(2)初始应力消除：为使模型的颗粒性质与标定后的土体性质一样，在土体颗粒的参数进行赋值前对模型已存在的应力进行消除，运用PFC中自带的“伺服机制”对模型应力进行调整，调整后模型初始应力小于某一设定值，且分布均匀。

(3)悬浮颗粒的消除：与周围颗粒的接触数目为小于2的颗粒称为悬浮颗粒，对模型颗粒进行遍历，对悬浮颗粒进行标记并调整，直至所有颗粒的接触数目大于1。

(4)参数赋值：参考土体参数标定的结果，用已标定后的参数进对模型进行赋值。

(5)模型边界刻画：对赋值后的模型进行边界刻画，删除部分颗粒用来模拟初始土洞、溶蚀沟槽等。

(6)覆盖层塌陷计算：加载重力加速度，生成流体域，对覆盖层土体进行力学计算。

7.4数值模型边界及力学参数

PFC软件生成覆盖型岩溶数值模型，模型中心处的坐标为：(0，0)，模型范围为：X(-50，50)，Y(-30，30)；覆盖层的范围为：Y(-20，30)；基岩范围为：X(-30，-20)；平面X＝30为地面，设置为自由移动面，对模型边界进行设定：

(1)在X(-50，-49)范围内，约束为X＝-49的直线决定，位移、速度约束为0，标记为不排水边界；

(2)在X(49，50)范围内，约束为X＝49的直线决定，位移、速度约束为0标记为不排水边界；

(3)在Y(-30，-29.7))范围内，约束为X＝-49的直线决定，位移、速度约束为0，标记为不排水边界。

7.5力学参数及覆盖层参数折减

本发明中覆盖层和基岩的本构模型均为接触黏结模型，在本发明对覆盖层进行了参数标定，覆盖层的力学参数取值为：颗粒法向刚度为10⁷N/m，颗粒切向刚度1e7N/m，摩擦系数0.3，法向黏结强度8×10⁴Pa，切向黏结强2×10⁴Pa，传导系数为5×10^-11cm/s。在岩溶塌陷过程中，基岩不参与水力学计算，为防止基岩力学强度过小导致基岩破坏，基岩力学强度设置较大：法向黏结强度8×10⁴⁸Pa，切向黏结强度2×10⁴⁸Pa。模型中有左右加载水压的情况，为防止加压导致模型边界破坏，对左右边界的黏结强度赋予较大的数值，法向黏结强度8×10⁴⁸Pa，切向黏结强2×10⁴⁸Pa。

土洞发育过程中，土体强度在水流的软化作用发生变化，本发明中数值模型通过保持覆盖层土体的强度参数不变，加大系统重力加速度g，增大土体的重力来表示土体强度在水流的软化作用下所受到的影响。土体重力增加会致使土体正应力的增长速度大于剪应力，当正应力的增长速度大于剪应力时，覆盖层土体可能不发生塌陷，为了防止这种情况的发生，多次调节重力加速度的增长率，使重力加速度增长率适合整个模型。

7.6模型破坏的标准判断

岩溶塌陷过程是一个缓慢破坏过程，与土体在水流软化作用下强度变小有关。由于土体内部发生应力集中，土体颗粒所受应力大于最大法向黏结强度、最大切向黏结强度时，土体强度发生破坏，土体内部产生裂纹，裂纹贯通逐渐形成连通面，土体颗粒沿着连通面剥落、坍塌，直至整个区域塌陷。

覆盖型岩溶塌陷模拟过程中，地面是否稳定是根据覆盖层土样内部是否发生破坏来判断的，故土体的破坏判别十分重要，此次模拟过程中，土体是否破坏通过下面三种情况进行判别：一、裂纹大量出现且形成渗流通道；二、计算过程中系统不平衡力的不收敛；三、土体位移大幅度增加。岩溶模型同时满足

上面三种情况时，才可视为覆盖层土体内部发生了渗流破坏。

离散单元法数值模拟过程中，土样内部颗粒黏结发生断裂产生裂纹，或颗粒位移大幅度突变，并不能说明土样发生了渗透破坏^[57]。由于离散元颗粒流数值模拟方法自身的特殊性，用裂纹多少、位移法判断并不全面。在土洞的形成过程中，只有土样内部的裂纹形成一个完整的渗流通道，土体颗粒大规模流失，才能认为土体内部发生破坏；土洞阶段性扩展过程中，需观察颗粒位移的变化、系统不平衡力；在模拟过程中，土洞完全破坏时，裂纹数目应大幅度增长、土体颗粒位移应发生较大变化，系统不平衡力应不收敛。在整个覆盖型岩溶溶塌陷模拟过程中，应以裂纹的贯穿程度、系统不平衡及颗粒位移的大小作为模型破坏的判别标准。因此模拟过程中，在覆盖层土体内设置的监测变量为：裂纹数目、系统平均不平衡力、颗粒位移。

8、土洞形成临界地下水流速模拟

岩溶塌陷灾害与地下水水位变化有关，地下水位在基岩面与覆盖层交界处发生强烈变化，致使基岩面附近土体强度变小、在地下水的浮力作用下，土体之间的有效应力变小。若地下水大幅度变化，此时在土岩交界面附近溶蚀沟槽开口处易产生土洞。为分析地下水位下降时，在土岩分界面附近土洞发育起始阶段的临界地下水流速，用PFC进行模拟，在此过程中分析覆盖层中不同位置的颗粒位移、裂纹扩展情况、系统不平衡力的变化情况，计算得出土洞在初始形成阶段开始扩展的临界地下水流速。

8.1模型工况设定

为模拟抽水条件下，地下水位突降时，覆盖层土体颗粒在基岩裂隙溶蚀沟槽顶部开口处发生剥落并形成土洞的临界条件，运用PFC建立数值模型，此模型说明如下：

(1)模型尺寸为100mm×60mm，颗粒半径大小为0.1-1.8mm，且均匀分布，颗粒之间采用接触黏结。

(2)土洞上覆盖层厚度为50mm(PFC主要力学计算区域)，岩性为粘土，基岩地层为10mm(不参加力学计算)。此模型中，基岩裂隙统设为宽12mm，高为13mm，覆盖层土体中发育高为0.3mm的裂隙。

(3)为模拟由抽水引起的水位突变导致盖层土体在基岩裂隙溶蚀沟槽顶部开口处发育土洞，此数值模型的边界条件设置为定水头边界，模型运行中边界水位恒定，由于模型较小忽略位置水头，以水压大小代替地下水水位大小。

(4)对模型左右两侧的颗粒给定恒定水压，且水流主要作用在覆盖层中，土体颗粒施压范围中，范围为： X的范围(-50，-35)和(35，50)，Y的范围(-20，12)；土体在Y(12，30)的范围内无水流作用，只受重力作用。

(5)模型中水体最终排泄点为溶蚀沟槽，此处水压大小设置为0，排泄点范围为：X(-12，12)，Y(-25，-18)。

8.2数值模型运行过程及结果

假定在裂隙处存在一口抽水井进行抽水，为模拟得到土洞形成的临界地下水流速度，对两侧水压逐级调整，直至在土岩交界面附近溶蚀沟槽开口处的土体开始发育土洞。对水压进行差值递增，在调试过程中，当水压大小为2.25×10⁵Pa时，盖层土体颗粒之间的接触发生断裂、颗粒剥落，此时的水压作为加载水压。

设定水压大小后，当模型运行到2.4×10⁵时步时，在溶蚀沟槽开口处的土体颗粒位移增加、土体内部产生裂纹，系统不平衡力不收敛，此时，覆盖层土体内部发生了渗流破坏。其破坏过程如下：土体颗粒间的接触黏结被破坏形成裂纹，裂纹数目迅速增长形成连通面，土体颗粒沿着连通面运动并在重力作用下掉落，此时，覆盖层土体发生部分坍塌，在溶蚀沟槽开口处形成土洞。临界地下水流速度-数值模拟破坏图中，裂纹出现在溶蚀沟槽开口处附近的土体内部，土体颗粒沿着裂纹连通面剥落、覆盖层内部发生坍塌形成初始土洞。

8.3裂纹扩张及模型破坏判定

通过上述分析可知，以裂纹数目、颗粒位移、系统不平衡力其中的一个参数作为土体是否发生破坏的判别依据是不准确的，本发明将讨论裂纹的贯穿程度、系统不平衡及颗粒位移三者间的联系，并得出模型发生渗流破坏的时间。

(1)模型不同位置颗粒位移

对在土岩交界面附近溶蚀沟槽开口处的土体颗粒的位移进行监测，监测颗粒分别为溶蚀沟槽中间位置的颗粒1和溶蚀沟槽两侧的颗粒2。本模型中设定颗粒位移大于0.5时颗粒自动删除，代表颗粒已经剥落。颗粒位移变化图如图11、图12所示，在240000时步之前，颗粒1、2的位移随着时间的增长而增长，且颗粒2的增长速率较大；模型运行到240000时步时，颗粒2的位移急剧增加，位移值大于临界值，此颗粒已剥落。由于土体坍塌产生的水压扰动颗粒1位移小幅度下降；240000时步后，溶蚀沟槽开口两侧处的土体颗粒1的位移随着时间的增长而增长，而颗粒2已剥落则无位移曲线。

由上述分析可知：抽水作用下，在土岩交界面附近溶蚀沟槽开口处的土体颗粒位移逐渐增大，中间位置的颗粒开始剥落，土洞开始发育并有向上扩展的趋势。因覆盖层土体颗粒剥落而产生了扰动域，导致土岩交界面附近溶蚀沟槽开口处两侧颗粒位移发生突变，随后两侧颗粒位移继续增加直至剥落。

(2)裂纹扩展与贯通

模型塌陷过程中土体颗粒间的接触力链图所示，线条代表接触力链且力链越粗说明接触力越大、土体颗粒所受外力越大。短线为颗粒接触黏结发生破坏所产生的裂纹。

如图13临界地下水流速度—裂纹数目变化曲线图所示：裂纹产生于覆盖层与溶蚀沟槽交界面处的土体内，在2400000时步之前，模型产生的裂纹个数为1；模型运行到240000时步时，裂纹数目急剧增长且快速发展为连通面，此时土体颗粒发生剥落、坍塌，颗粒间的接触力消失。在裂纹形成之后，裂纹附近的接触力链逐渐变粗说明此时发生了应力集中、并形成明显的应力拱。

(3)模型系统不平衡力

对特殊位置的颗粒位移进行了分析，得到模型破坏时颗粒位移的变化趋势，但用颗粒位移的变化趋势来判别模型是否破坏的做法具有一定局限性，本发明中也对系统不平衡力进行了监测。

如图14中a所示，在模型建立初期，颗粒刚生成时，颗粒之间的重叠量较大，导致系统不平衡力较大，运用“伺服机制”对模型进行调整，调整后模型的不平衡力大小在60N左右小幅度来回波动，说明在此条件下模型达到了自然平衡。如上图b所示，在较长一段时间内，系统不平衡力波动幅度较小，并保持在某一范围值内；220000时步时，系统不平衡力小幅度增长；240000时步时，系统不平衡力急剧增长，不平衡力的数值大小达到了10⁶N，随后在10⁵N附近波动，最后恢复到平衡状态。

综合对颗粒位移、系统不平衡力、裂纹数目进行分析，可得：在240000时步时，颗粒位移、系统不平衡力、裂纹数目皆有较大的变化，且三者之间相呼应，可认为240000时步为模型发生渗流破坏的时间点，此时，裂纹急剧增长发育成连通面、土体颗粒顺着连通面剥落、不平衡力急剧增长，覆盖层土体部分坍塌，土洞开始发育。

8.4临界地下水流速度模拟结果：

此模型模拟土洞的发育过程，并求取土洞初始形成的地下水临界地下水流速度。本发明对模型出口处的流量进行监测，监测出口处的横断面积A＝48mm，流体记录的时间步长为0.1s，出口流量和时间的关系如上图所示，综合对颗粒位移、裂纹数目、系统不平衡力的分析可知240000时步时，土体颗粒剥落、覆盖层发生部分坍塌，此时覆盖层土体中开始发育土洞，如流量—时间图15所示，240000时步流量先急剧减少后急剧增加，此时模型监测的流量大小为：5×10⁵mm³，由达西定律可得：

带入数值、单位换算可得模型破坏时的临界地下水流速度为v＝0.104cm/s。

8.5结果分析及对比

由粘土的临界地下水流速度公式：

计算理论临界地下水流速度，其中n₀＝0.2，e₀＝0.25，d₀＝0.25mm(在总模型中，颗粒粒径在d<0.25mm 的范围内的颗粒数目较多)，μ＝1.01×10^-6kPa/s，γ_s＝19.6kN/m³，γ_w＝9.86kN/m³，G＝2.5，A＝10^-31，α＝0.5，计算可得：d₁＝1.33×10^-5，ν_cr＝0.113cm/s，

误差：

数值模拟得到的临界地下水流速度为：v＝0.104cm/s；理论值为：ν_cr＝0.113cm/s，两者之间较为接近，理论值与模拟流速的误差小，为8％，通过对比模拟值和理论值可得：PFC计算得到的土洞扩展临界地下水流速度是合理的、可靠的，从而用PFC分析岩溶塌陷问题是可行的，为其他条件下的岩溶塌陷问题提供理论依据。

9、岩溶塌陷土洞临界高度模拟

覆盖型岩溶塌陷的主要动力来源为地下水水位变化。分析区内覆盖层厚度不大，土体力学强度小，地下水水位波动频繁，导致覆盖层土体易发生塌陷。在地下水的作用下，土体颗粒被水流带走，土洞逐渐扩展，土洞规模扩展到某一范围时，覆盖层顶板发生瞬时塌陷，造成地面塌陷。为得到在大气降雨或地表水垂直入渗条件下，覆盖层失稳坍塌的过程，用PFC软件进行模拟分析。用PFC数值模拟软件得到塌陷过程中覆盖层颗粒位移、裂纹数目，系统不平衡力的变化，并得出土洞扩展过程中的最大临界高度。

9.1模型工况设定

为模拟在降雨条件下，覆盖层中存在初始土洞的情况下，土洞的扩张、塌陷过程，并求出在地面塌陷之前的最大土洞高度，运用PFC数值软件进行模拟，模型介绍如下：

(1)模型尺寸为100mm×60mm，颗粒半径大小为0.1-1.8mm，且均匀分布，颗粒之间采用接触黏结。

(2)覆盖层厚度为50mm(PFC主要力学计算区域)，基岩地层为1mm(不参加力学计算)。此模型中，溶蚀沟槽统设为长12mm，宽为5mm。

(3)初始土洞的圆心为(0，-40)，半径为8mm，跨度为49mm，高度为13mm。

(4)对处于地面附近的颗粒进行加压处理来模拟地表水的入渗作用，考虑到边界对模型的影响，本模型中假设水流已进入土体，且在某一范围内产生连续的、恒定的水压，模型加压范围为：X(-45，45)， Y(30，20)。

(5)模型的水体最终排泄点为溶蚀沟槽，此处水压设置为0，其范围为：X(-12，12)，Y(-25，-23)。

9.2模型过程及分析

模型中土体颗粒的颜色代表着此颗粒的位移大小，红色颗粒显示此颗粒的位移较大，蓝色颗粒显示此颗粒的位移较小，颗粒位移越大说明土体受到扰动且扰动力越大。

可得出，土体塌陷初始阶段在土洞顶端最先形成连通面，颗粒沿着连通面掉落，土层受到扰动导致模型的初始平衡状态发生变化，土洞顶端的颗粒在扰动的作用下逐渐剥落，土洞向上扩张。土洞扩展至某一阶段时，模型呈现出相对稳定状态，土洞高度保持不变且达到最大临界高度。对模型持续加压，一段时间后，覆盖层发生瞬时塌陷，塌陷从覆盖层内部迅速扩展至地表致使地表发生塌陷。地表塌陷后形成了一个完整的塌陷坑，坑壁为直筒型，且塌陷坑两侧呈对称形状。上述过程可描述为：雨水入渗——土洞扩张——地表塌陷。

本模型中设定颗粒位移大于0.5时自动删去，代表颗粒已经剥落，通过对覆盖型岩溶塌陷过程图进行以下分析：

(1)模型运行初期，地表水垂直入渗地层，地表附近的土体颗粒受外力的影响较大，因此地表附近土体颗粒的位移较大。

(2)模型运行到15000时步时，水压逐渐到达溶蚀沟槽处，在垂向渗透压力的影响下，土洞顶部的颗粒位移逐渐增大，颗粒间的接触黏结发生破坏，游离土体颗粒在重力作用下剥落、土体内部发生部分坍塌，土洞向上扩张；如此同时，左右两侧的土颗粒受到扰动、土体变形较大，扰动区向两侧逐渐扩大，且破坏土体厚度以外的土体发展为第一级土层顶部；

(3)模型运行到16500时步时，覆盖层在垂向渗透压力和重力作用下，土洞两侧变形较大的土体颗粒发生剥落，土洞向两侧扩张，发展为第二级土洞；

(4)模型运行到330000时步时，在垂向渗透压力和重力作用下，扰动变形区的土体颗粒已全部剥落、坍塌，土洞继续向上扩张，此时土洞发展到覆盖层中间位置，为三级土洞；此阶段发育过程用时较长，土洞达到了一个相对稳定的状态；

(5)模型运行到335000时步时，覆盖层内部发生瞬时塌陷，在较短时间内(模拟运行时间在335000— 360000时步内)，土洞扩张至地表表面导致地表塌陷。地表塌陷过程用时较短，塌陷速度快，且地表塌陷并不是覆盖层整体坍塌而是从覆盖层内部逐步发展到地表的。

以上分析表明：岩溶塌陷过程中，土洞扩张是从土洞顶部逐级扩展的，其发育过程为：向上扩张→两侧扩张→继续向上扩张→地表塌陷。并得到当土洞发展到某一阶段时，土洞不在扩张，模型达到临界状态，此时覆盖层直接塌陷，且坍塌过程是从内部发展至地表的。

9.3裂纹扩张及模型破坏判定

(1)模型不同位置颗粒位移

对岩溶塌陷模型中特殊位置的颗粒位移进行监测，颗粒位置分别为：土洞洞顶附近的颗粒3，土洞右侧附近的颗粒4，覆盖层中间位置的颗粒5，地表附近的颗粒6，对塌陷模型过程进行分析。因本模型中设定颗粒位移大于0.5时自动删除此颗粒，故测量颗粒位移时，位移到达0.5时，测量自动停止。

分析可得：14000时步时，土洞向上发展、覆盖层发生局部崩塌；16500时步时，土体发生崩塌，土洞向两侧扩展；330000时步时，土洞继续向上扩展；335000时步时，土洞不在扩展，覆盖层发生瞬时塌陷。颗粒位移图证实了塌陷过程，并且很好的说明了覆盖型岩溶塌陷过程中土洞的发展具过程具有阶段性和方向性，也验证了土洞发展到某一阶段不在扩展而是直接塌陷。

覆盖型岩溶塌陷过程中颗粒之间的接触力链、裂纹发展阶段图可知：(1)模型运行到14000时步，地表水入渗土层，土洞顶部的土层发生应力集中，在应力集中处，颗粒所受合力大于颗粒的接触黏结强度导致土体结构破坏，土洞顶部最先出现裂纹；(2)模型运行到16500时步时，土洞两侧的土体中发生了应力集中并形成了明显的应力拱，土体中出现了大量的裂纹，裂纹逐渐贯通形成了连通面，土体颗粒沿着连通面剥落，土洞向上扩展；(3)模型运行到330000时步时，土洞两侧的应力拱逐渐消失，裂纹增长缓慢并形成新的连通面，导致土洞向两侧扩张，土洞拱趾处发生应力集中；(4)在335000—350000时步内，模型中产生大量裂纹且裂纹迅速扩展、贯通形成新的连通面；(5) 模型运行360000时，裂纹急剧上升并形成了垂直连通面，土体沿着垂直连通面发生滑动，整个覆盖层迅速发生塌陷；(6)覆盖层完全塌陷后，模型达到稳定，裂纹不在继续增长，且应力拱消失。

以上现象表明：(1)裂纹形成的原因为：土体颗粒接触处发生了应力集中导致颗粒所受应力大于颗粒间的接触黏结，黏结发生断裂形成裂纹；(2)土洞的发育过程可分为两个阶段：裂纹形成阶段和裂纹形成连通面阶段。

覆盖型岩溶塌陷过程为：覆盖层中发生应力集中导致土体结构发生破坏并形成裂纹，裂纹逐渐贯通成连通面，土体顺着连通面剥落，土洞发生扩张，此过程循环直至覆盖层塌陷，模型塌陷后应力集中消失达到平衡状态。此过程应力与裂纹的耦合可近似认为：“应力平衡→应力集中→裂纹形成→裂纹形成连通面→颗粒剥落→应力平衡”，伴随着覆盖层塌陷，模型应力达到稳定

(2)模型不平衡力

模型不平衡力变化图所示，塌陷模型未发生破坏时，不平衡力较小，且浮动的范围较小，此时模型相对稳定；模型运行到14000时步时，不平衡力发生较大的浮动，其最大值达到了75000N。不平衡力的急剧增长说明了塌陷模型稳定性发生了变化，土层内部发生了渗透破坏，土体内部初次发生渗透破坏的时间段为14000—20000时步；土层内部初次发生渗透破坏后，系统不平衡力波动幅度较小，塌陷模型恢复到相对稳定状态，模型稳定期的时间段为20000—60000时步；在60000—320000时步段内，系统不平衡力多次发生小幅度波动，说明模型在此时间段内发生了小型渗透破坏；如图d所示，时间范围为 320000—350000时步，系统不平衡力波动频率高、幅度大，在350000时步附近不平衡力最大值达到了800000N，覆盖型岩溶模型发生坍塌。

综上对颗粒位移、裂纹发展规律、系统不平衡力的分析对比可知：(1)不平衡力波动幅度较大，渗透破坏程度越大。(2)覆盖型岩溶的塌陷过程中可分为四个阶段：第一阶段(14000时步时)，土洞向上扩展，形成第一级土洞，此时裂纹形成小型连通面；第二阶段(16000时步时)，第二级土洞形成，裂纹增长较快且贯通成大型连通面，覆盖层土体颗粒发生较大规模的崩塌；第三阶段(330000时步时)，第三级土洞形成，此过程中耗时较长，土洞发展较为缓慢，覆盖层中多次发生小型坍塌，此阶段可得到土洞最大临界高度；第四阶段(330000—360000时步之间)，土洞规模不在扩大，裂纹数目急速上升且形成垂连通面，地表发生塌陷。

9.4土洞临界高度模拟结果

由上述模拟可得，在覆盖型岩溶模型中存在初始土洞的情况下，土洞发育过程具有阶段性。当土洞扩展到某一程度时，覆盖层土体不在发生渗透破坏、土洞高度达到最大临界值，在雨水入渗和土体重力作用下覆盖层将发生瞬时塌陷。由前人分析的可得地表发生瞬时塌陷前的土洞高度为临界土洞高度，由上述模拟结果可知，模型运行到330000时，此时土洞高度为最大临界值；

可知，土体发生破坏前的临界土洞形状等效为圆弧形，此时压力拱(圆弧)的跨度D＝48mm，测量模型中基岩面到土洞顶端的距离为22.9mm。覆盖层土体厚度为H，土洞的极限平衡高度h，由土洞临界平衡高度的应用前提为H>2h，此模型中覆盖层土体厚度H＝55mm，土洞极限平衡高度2h＝45.8mm，符合前提条件H>2h。

9.5结果分析及对比

(1)临界高度模拟值与理论值

由王滨提出了临界土洞的极限平衡高度的计算方法，得出极限临界土洞高度和土洞跨度之间的关系式，如式所示：

其中：h_max为土洞的极限平衡高度(m)；D为土洞的跨度(m)；f_k为盖层土体的坚固性系数，其中中砂、细砂：

f_k＝0.5；砂质黏土：f_k＝0.6；粉质粘土、粘土：f_k＝0.8，

此模型中覆盖层土体为粘土f_k取为0.8。

其中覆盖层土体为粘土，f_k＝0.8，带入D＝48mm到公式(5.2)中，可得：

H＝55>2h_max＝49.6

由计算可得，理论土洞的极限高度与实际土洞极限高度的误差为7.6％，误差小于15％。

(2)模拟塌陷现象与实际现象

PFC数值模型塌陷图为近似直筒式塌陷，因其边界条件影响，塌陷形状不太规则，以实际塌陷图对比， PFC模拟的岩溶塌陷在宏观现象上符合实际情况。由上文分析可知，覆盖层塌陷前土洞的跨度D约为13m、临界土洞高度约H为6m，故土洞跨度D和最大临界高度H的比值为I：

I_理论＝0.52

通过比较土洞跨度D和最大临界高度H比值的实际值、理论值、模拟值，可得三者之间较为接近，说明PFC计算结果是符合理论与实际的。通过验算说明PFC计算得到的岩溶塌陷过程中的土洞临界高度的是符合理论的、可行的，并说明用颗粒流数值模拟方法得出覆盖型岩溶塌陷的裂纹发展、颗粒剥落、地表塌陷等现象是符合实际的。

本发明运用PFC数值软件，模拟降雨或地表水垂直入渗导致覆盖层土体在土岩交界面附近溶蚀沟槽开口处发育土洞，模拟地下水渗流引起的地表塌陷，在分析分析过程中得到如下结论：

运用PFC模拟渗流模型时，对于模型是否发生渗流破坏，应结合土体颗粒位移曲线、裂纹增长曲线、不平衡力变化曲线对模型进行判断，只有满足这三种情况时才可确定模型发生了渗流破坏。

利用PFC数值软件模拟土洞初始发育阶段时，得到土洞形成的临界地下水流速度，并与理论值作对比，两者之间相差较小，说明模拟结果与理论结果一致。说明PFC软件模拟岩溶塌陷问题是可行的，且所呈现出来的微观过程是合理可信的，从而可实现了其他方法(物理模型、解析法等)无法实现的微观破坏过程的目的。

运用PFC软件对覆盖型岩溶塌陷发展过程进行模拟，其结果表明：

①覆盖型岩溶塌陷过程中的应力变化情况为：应力平衡→应力集中→裂纹形成→裂纹形成连通面→颗粒剥落→应力平衡。土体中发生应力集中，颗粒所受应力大于土体强度时，接触破坏、形成裂纹，颗粒沿着连通面剥落后，模型达到短暂的稳定状态。

②对覆盖型岩溶塌陷进行数值时，可把塌陷过程分为四个阶段：第一阶段，土洞向上扩展，形成第一级土洞；第二阶段，土洞向两侧扩展，形成第二级土洞；第三阶段，土洞发展较为缓慢，覆盖层中多次发生小型坍塌，第三级土洞形成；第四阶段，土洞规模不在扩大，裂纹数目急速上升且形成垂连通面，地表发生塌陷。

③得到覆盖型岩溶塌陷过程中的微观现象：如裂纹扩张、土体颗粒位移、土洞发育、地表塌陷等过程。

④利用PFC模拟软件计算得到土洞最大临界高度，与理论值对比，误差较小，符合理论情况。

⑤在微观过程中，土洞高度达到最大临界值时，地表迅速塌陷，用时较短，且塌陷从内部延伸至地表，并不是覆盖层直接塌陷。

现有的以―连续介质力学”为理论基础的数值模拟方法(例如：(FLAC)快速拉格朗日差分法)，仅仅能够模拟覆盖型岩溶区盖层土体(土洞)的大变形、挠曲和塑性流动等力学行为；

下面结合效果对本发明作进一步描述。

土洞的盖层土体由破碎的固体颗粒组成，土的宏观变形主要不是由于颗粒本身变形，而是由于颗粒间位置即排列方式的变化造成的，传统的基于“连续介质力学”理论的数值模拟方法(例如：(FLAC)快速拉格朗日差分法)无法体现土颗粒的动态迁移、潜蚀以及崩解过程。

而本发明所述的基于“非连续介质力学理论”的离散单元法，克服了传统连续介质力学模型的宏观连续性假设，可以从微观角度对土的工程特性进行数值模拟，通过颗粒结构的微观参数研究来分析材料的宏观力学行为，体现塌陷过程中土颗粒潜蚀、迁移、崩落等基于非连续介质假设的力学行为,如图17基于“非连续介质力学”理论的数值模拟则可以很好的体现土体颗粒迁移、崩落等塌陷过程中的真实现象。从而更加真实的模拟塌陷体由连续体到非连续体的渐进破坏过程。如图18基于“非连续介质力学”理论的离散单元数值模拟方法技术效果图。

在上述实施例中，可以全部或部分地通过软件、硬件、固件或者其任意组合来实现。当使用全部或部分地以计算机程序产品的形式实现，所述计算机程序产品包括一个或多个计算机指令。在计算机上加载或执行所述计算机程序指令时，全部或部分地产生按照本发明实施例所述的流程或功能。所述计算机可以是通用计算机、专用计算机、计算机网络、或者其他可编程装置。所述计算机指令可以存储在计算机可读存储介质中，或者从一个计算机可读存储介质向另一个计算机可读存储介质传输，例如，所述计算机指令可以从一个网站站点、计算机、服务器或数据中心通过有线(例如同轴电缆、光纤、数字用户线(DSL)或无线 (例如红外、无线、微波等)方式向另一个网站站点、计算机、服务器或数据中心进行传输)。所述计算机可读取存储介质可以是计算机能够存取的任何可用介质或者是包含一个或多个可用介质集成的服务器、数据中心等数据存储设备。所述可用介质可以是磁性介质，(例如，软盘、硬盘、磁带)、光介质(例如，DVD)、或者半导体介质(例如固态硬盘SolidState Disk(SSD))等。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法，其特征在于，所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法包括：

在覆盖型岩溶区典型塌陷的地质结构概化模型基础上，利用颗粒流理论数值模拟方法对土洞的形成过程、覆盖型岩溶的塌陷过程进行模拟，获取初始土洞形成的临界地下水流速度和覆盖型岩溶塌陷的土洞临界平衡高度。

2.如权利要求1所述的基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法，其特征在于，所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法具体包括：

第一步，运用PFC模拟渗流模型，结合土体颗粒位移曲线、裂纹增长曲线、不平衡力变化曲线判断是否发生渗流破坏；

第二步，利用PFC数值软件模拟土洞初始发育阶段，得到土洞形成的临界地下水流速度，并与理论值作对比；

第三步，运用PFC软件对覆盖型岩溶塌陷发展过程进行模拟：

对覆盖型岩溶塌陷过程中的应力变化情况进行模拟；

对覆盖型岩溶塌陷进行数值模拟，将塌陷过程分为四个阶段：第一阶段，土洞向上扩展，形成第一级土洞；第二阶段，土洞向两侧扩展，形成第二级土洞；第三阶段，土洞继续发展，覆盖层中多次发生小型坍塌，第三级土洞形成；第四阶段，土洞规模不在扩大，裂纹数目急速上升且形成垂连通面，地表发生塌陷；

覆盖型岩溶塌陷过程中的微观现象模拟；

利用PFC模拟软件计算土洞最大临界高度，与理论值对比。

3.如权利要求1所述的基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法，其特征在于，第一步的模拟过程包括：应力平衡→应力集中→裂纹形成→裂纹形成连通面→颗粒剥落→应力平衡。

4.如权利要求1所述的基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法，其特征在于，覆盖型岩溶塌陷过程中的微观现象模拟中，包括裂纹扩张、土体颗粒位移、土洞发育、地表塌陷过程的模拟。

5.如权利要求1所述的基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法，其特征在于，离散元的微观耦合力学方程：

PFC中流固耦合计算的流动方程、压力方程和求解条件如下：

流动方程流体管道相当于一个平行通道，长度为L’、孔径为a，在垂直平面方向上为单位厚度，管道内的流速为：

K‘—传导系数(cm·s^-1)；L’—管道的长度(mm)；P₂-P₁—相邻域的压力差；压力方程周围管道流入每个域的流量和为∑q，在单位时间步长Δt下，流体压力增量Δp(流入为正)为：

式中：K_f——流体的体积模数(kPa)；V_d——域的表观体积(mm³)；

应用显示求解方法，将流量方程应用于所有的管道，并将压力方程应用于所有的域之间交替求解；假设某个域内存在扰动压力ΔP_p，由于扰动流域里的流量计算得：

由水流流入引起的响应压力变化

式中：N——域所黏结的管道数；R——域周围颗粒的平均半径(mm)；

当两者相等时求出临界时间步长为：

6.一种实现权利要求1～5任意一项所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法的计算机程序。

7.一种实现权利要求1～5任意一项所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法的信息数据处理终端。

8.一种计算机可读存储介质，包括指令，当其在计算机上运行时，使得计算机执行如权利要求1-5任意一项所述的基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法。

9.一种实现权利要求1～5任意一项所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟方法的基于离散单元法的覆盖型岩溶塌陷灾变演化模拟控制系统。

10.一种搭载权利要求9所述基于离散单元法的覆盖型岩溶塌陷灾变演化模拟控制系统的计算机。