基于LBM-DEM 细观数值模拟的水力诱导覆盖型岩溶地面塌陷发育过程分析

陶小虎，叶 明，龚建师，王赫生，胡晓雨

（1.中国地质调查局南京地质调查中心，江苏 南京 210016；2.自然资源部流域生态地质过程重点实验室，江苏 南京 210016；3.佛罗里达州立大学地球、海洋与大气科学学院，佛罗里达州塔拉哈西 32306；4.江苏省水资源服务中心，江苏 南京 210029）

0 引言

覆盖型岩溶地面塌陷的形成是自然环境变迁的一个组成部分，而人类活动对岩溶发育地区自然平衡状态的干扰增多，大大加快了这一地质灾害发生的进程，由此引发的地面塌陷问题日趋频繁、后果严重，已成为岩溶地区可持续发展的一大障碍[1-2]。在众多覆盖型岩溶塌陷形成过程中，由地下水位变化引起的岩溶塌陷居多，地下水流提供了最主要的动力。地下水过量开采、水库蓄水等引起的地下水位变化为主要诱发因素[3-6]。

长期以来，国内外学者对地下水位变化引起的覆盖型岩溶地面塌陷问题一直十分关注并开展了卓有成效的研究。在数学解析方面，基于渗透力学、极限平衡理论、牛顿摩擦定律等方法建立简单力学模型，获得宏观规律，确定岩溶塌陷的机制及其临界条件的研究[7-9]。在物理实验方面，覆盖型岩溶塌陷模型试验研究从单一致塌因素研究到多因素组合、从简单的二维定性研究到三维定量分析以及从真实土到透明土可视化分析都体现了试验研究的巨大进步[10-12]。在数值模拟方面，基于连续介质理论，在岩溶塌陷临界条件、主控因素、塌陷模式等内容的研究上取得了较好的成果[13-14]。然而由于受到当前科学技术的局限、土洞发育以及塌陷过程的非连续、非均匀、流固耦合等特性，现有的物理实验和数值模拟方法很难对水力诱导覆盖型岩溶塌陷进行系统的岩溶塌陷发育和塌陷过程的细观研究，对各影响因子之间的定量关系研究也受到了限制[15]。

离散元方法（discrete element method,DEM）是根据离散物质本身的离散特性建立数值模型，与连续介质理论方法相比，DEM 可以从细观角度完成大变形分析，更好地描述结构破坏过程[16-17]。格子Boltzmann 方法（lattice Boltzmann method，LBM）是一种新兴的计算流体学方法，能够从介观尺度上描述流体的运动和模拟多相多尺度问题。与传统流体力学计算方法相比，LBM在程序实现上更为简单，能够实现与大量移动颗粒的完全耦合，很好地模拟多相多尺度耦合问题，尤其是在模拟复杂边界条件的多孔介质水流中，更能体现出其相比于传统流体力学计算方法的优势[18]。因此，为了能够从细观上对覆盖性岩溶塌陷模拟研究，本文将采用DEM模拟土体颗粒运动，选择库仑定律来判别土体颗粒之间的黏性接触是否受到破坏，采用LBM 中D2Q9 模型[19]模拟土体孔隙间水流运动，应用Bouzidi 等[20]反弹边界格式和动量交换法实现LBM 与DEM 双向耦合，建立覆盖型岩溶塌陷二维流固耦合细观数值模型，并对承压水下降工况下的覆盖性岩溶塌陷过程进行探索研究。

1 数值方法

本节主要介绍LBM-DEM 双向流固耦合的方法，主要包括固体对流体的作用、流体对固体的作用、耦合模型时间步长同步。

1.1 固体对流体的作用

为了计算固体颗粒与流体之间的相互作用的准确性，特别是流体对固体颗粒的作用力的计算，本文采用具有二阶精度且更为简单的线性插值的Bouzidi 插值反弹边界格式（图1），假设每个时步末节点A处-α方向的分布函数是由虚拟流体节点D处的粒子向右迁移，经过点C反弹运动至节点A所得，其中一个时间步长所迁移长度为一个空间步长，即 |DC|+|AC|=|AB|。取q=|AC|/|AB|，fα(r,t)为在r位置离散速度在 α方向上粒子密度分布函数。当q<时，虚拟节点D处于流体节点E、A之间，可由流体节点E、A处的分布函数插值得到虚拟节点D的分布函数；当q≥时，虚拟节点D处于流体节点A与固体节点B之间，可以通过节点A处的碰撞后的分布函数进行线性插值得到点A 处下一时步的-α方向的分布函数。

图1 Bouzidi 插值反弹边界示意图Fig.1 Schematic view of Bouzidi interpolation bounce-back diagram

1.2 流体对固体的作用

动量交换法能够精确有效且便捷的处理流体流动对固体的力学作用的计算[21]。因此本文将采用动量交换法评估二维流体对固体的力和力矩。在评估流体对固体颗粒的作用力和力矩时，固体颗粒中的分布函数不参与计算。

式中：rf——颗粒边界流体节点位置；

Fσ——流体节点处流体对固体颗粒的作用力；

Tσ——流体节点处流体对固体颗粒的作用的力矩；

δx——LBM 中格子步长；

eα——速度方向；

δt——时间步长；

fα(r,t)——在r位置离散速度在 α方向上粒子密度的分布函数；

rσ——作用力点到颗粒圆心的矢量。

1.3 耦合模型时间步长同步

在耦合模型中，涉及到两套方法的时间步长，作如下考虑：

当LBM 时间步长 δt小于或等于DEM 预估时间步长 δt时，DEM 时间步长取LBM 时间步长，δt=δt；

当LBM 时间步长 δt大于或等于DEM 预估时间步长 ∆t时，令n=，将DEM 时间步长更新为 ∆t=。

2 模型优化

2.1 二维颗粒的水力半径

实际覆盖层土体中，流体能够穿过由多个颗粒形成的多孔介质，即使相互紧密接触的颗粒之间也存在渗透路径。然而在二维LBM-DEM 耦合模型中，三维颗粒会被简化成二维圆形，相互紧密接触的圆形之间不存在任何通道，阻止了流体的流动，这与实际情况相悖。此外，三维颗粒与二维圆形的受力也不相同。因此在二维LBM-DEM 耦合模型中需要对颗粒进行处理。

Cui[22]从流体等效拖曳力出发，将三维球体颗粒等效为圆柱体，从而得到二维圆形颗粒的水力半径rh≈0.8R，其中R为颗粒的实际半径，该式适用于雷诺数较大（Re>1 000）的情况。本文将推导较小雷诺数的情况。

当雷诺数较小时，由斯托克斯阻力公式[23]可知三维中流体作用于圆球的力为：

同样，二维中，流体作用于圆盘的力为：

式中：U——实体颗粒相对于流体的速度；

F——三维实体颗粒收到流体作用力；

U2D——二维圆形颗粒相对于流体的速度；

F2D——二维圆形颗粒收到流体的作用力。

因此当球体颗粒在静水中自由沉降时，阻力应与速度成正比增大，逐渐达到匀速运动，作用于球体的重力、浮力以及阻力保持平衡，将二维和三维计算中的平衡速度进行等效，则有：

本文建议LBM-DEM 耦合模型模拟多孔介质时水力半径因子取：

2.2 颗粒直径与格子步长比优化

由于流体是在孔隙介质中流动，并对细小的孔隙介质产生力的作用，为保证流体运动和流体作用力达到精度要求，需要选择细小颗粒影响下的格子网络最优步长[24]。定义颗粒直径与格子步长比 β=D/δx，其中D为颗粒直径。本文通过平板间流体对圆形阻碍物作用力的收敛性分析，以选择计算合理的 β。数值模拟取平板长度为0.15 m，宽度为0.01 m，流体为20 °C 下的水，其密度为1×103kg/m3，动力黏滞系数µ为1.010×10-6kPa·s。左端流入压强3.4 Pa，右端流出压强为3.3 Pa。颗粒直径为1 mm，位于x=0.070 m、y=0.005 m 处。对3 种不同直径D（1.0 mm、1.2 mm、1.5 mm）的颗粒，计算分析不同β 下流体对颗粒的作用力，计算结果如图2 所示。

图2 颗粒所受到水流的水平作用力随β 的变化Fig.2 The variation of horizontal force on solid particle introduced by water flow with different β

颗粒受力计算值随着β 的增大而逐渐收敛。当β≥8 时，流体对颗粒作用力与参考值（最小步长0.05 mm工况）比较，相对误差可控制在0.5%以内。因此，对覆盖层多孔介质中进行流固耦合计算时，为确保计算的准确性和收敛性，在进行流体计算时，模型中经过水力半径因子缩放后的最小颗粒直径与格子步长比需满足β≥8。

3 覆盖型岩溶塌陷的建模

如图3 所示，覆盖型岩溶塌陷的流固耦合仿真计算共包括4 个部分：DEM 初始设置、LBM 初始设置、塌陷模型初始设置、塌陷模型计算。

图3 覆盖型岩溶塌陷的建模流程Fig.3 Modeling flowchart of covered karst ground collapse formation

（1）DEM 初始设置：锁定LBM 程序；确定模型的范围和边界，形成模拟空间；模型的左右边界、下边界为Wall 边界，模型内颗粒不能穿越此类边界。岩溶开口在前期设置中假设为Wall 边界，在后期修改为开放边界。上边界为开放边界，即颗粒可以自由出入该边界。在空间内按照四边形规则排布或者随机生成颗粒，并赋予颗粒力学参数（杨氏模量和内摩擦角，为了颗粒间紧密，此时不设黏性键）；引入局部阻尼[25-26]，在重力作用下将颗粒压实，经过足够长时间的循环，消除土体颗粒内部的不平衡力，让其达到相对稳定状态。

（2）LBM 初始设置：在DEM 初始设置完成后，激活LBM 计算程序，锁定DEM 程序（固体颗粒位置固定）；按照模拟范围生成离散网格；施加水力边界条件，计算固体颗粒间的水的流动情况，直至达到稳定；打开DEM 程序，修改DEM 程序中的时间步长，保证DEM和LBM 之间能够实时的同步交换，在水力和重力共同作用下，运行足够步数循环，使得土体达到相对稳定状态。

（3）塌陷模型初始设置：对已经生成的模型中的颗粒之间存在的接触按照参数黏性力学参数配置黏性键，经过足够长时间的循环，消除土体颗粒内部的不平衡力，完成颗粒流模型的黏性设置（黏性键只配置一次，后续计算中新生成的接触或者颗粒黏性键破坏后的接触不再拥有黏性属性）；去除模型中岩溶开口处的墙体，改为开放边界。开口处的颗粒由于失去墙体的顶托力在重力作用和应力释放影响下产生细微变形，少部分颗粒向下脱落，并形成平衡拱，运行足够步数，使得颗粒达到相对稳定状态。至此，完成对覆盖型岩溶塌陷模型的初始状态设置。一般情况下，拆除岩溶开口处的墙体（Wall）后的计算属于动力计算，离散元系统中的阻尼应设置为0，但由于岩溶开口的形成一般是非常缓慢的扩展过程，而模型中为了方便，将岩溶开口处的墙体直接拆除，容易造成模型系统内颗粒的流失，与实际状态不符。因此，在此步骤下可以仍然考虑继续采用阻尼来避免这一问题。

（4）塌陷模型计算：模型完成初始状态设置后，确保去掉阻尼（设为0）；根据计算需求修改流体边界条件；进一步进行耦合计算至稳定或者塌陷。（图4）

图4 覆盖型岩溶塌陷LBM-DEM 模型Fig.4 Schematic view of LBM-DEM model of covered karst ground collapse

在实际中黏土颗粒非常小，模拟数量巨大，普通计算机能力有限而不能有效模拟。颗粒放大法是目前较为可行的处理方式，已广泛应用于工程研究中，该方法将颗粒放大，降低模型中的颗粒数目，使原物理模型问题能在合理有效的时间内解决[27-30]。因此本文假设模型中的颗粒为多个黏性土颗粒结合的土体，以增加颗粒的半径。如图5 所示，模型空间范围取15.9 cm × 8.0 cm，岩溶开口大小为0.6 cm，颗粒直径取0.15 cm，对应LBM 中颗粒缩放后的直径为0.12 cm（水力半径因子取0.8），LBM 空间步长取 δx=0.000 1 m，即β 为12，生成127.2 万个网格，1 274 391 个格点。固体颗粒与流体具体参数主要参照文献[17]、[22]、[25]（表1）。在LBM中水头的计算可以用等效水压代替，令隔水层上方水压为Pt，下方岩溶开口处水压力为Pb，其余边界为水量平衡边界或不透水边界。

表1 计算模型参数Table 1 Summary of simulation model parameters

图5 承压水位下降工况覆盖型岩溶塌陷过程Fig.5 The process of covered karst ground collapse formation with declining confined water level

4 覆盖型岩溶地面塌陷的讨论

为了方便模拟，本次模拟工况为初始时刻模型中各处水头相等，等效成水压为933.33 Pa，Pb以降速50 Pa/s下降500 Pa 的情况，模拟结果如图5 所示。

（1）覆盖型岩溶塌陷形成过程

在初始时刻，模型中各处的水头相等，水流处于相对静止状态，对土体仅有向上的浮力作用。图5（a）为模型完成初始状态设置后的颗粒分布情况。

随着承压水的水位下降，岩溶开口处的水头低于土层中的水头，土层中孔隙水向下流动，对土体颗粒造成向下的力。由于水位下降初期水流对岩溶开口上方土颗粒的作用力较小，对颗粒未造成较大影响，随着承压水位的继续下降，孔隙水对土颗粒的作用力逐渐增加，岩溶开口上方土体颗粒逐渐向洞口处移动，并脱离内侧相邻土体颗粒，见图5（b）。如图5（c）所示，随着水力的作用，岩溶上方失稳土体颗粒增多，土洞不断扩展，但是此时黏性土层表面并不会出现明显的下沉。土洞不断扩大，在自身重力以及水力作用下，表层土体受到破坏并突然塌陷进下方的洞穴，见图5（d），而形成覆盖型岩溶塌陷，此后坠入土洞的土体在水力作用下逐渐通过岩溶开口流失，见图5（e），并最终形成落水洞。

比较图5（b）（c）（d）（e），土洞扩展由初期土体颗粒运动呈现单颗粒向下剥落，逐渐演变成块状失稳剥落，块体内的土体颗粒之间仅发生极小的相对位移甚至不发生相对运动，发生失稳破坏的块状构造不断变大，土洞扩展加速，最终导致土洞上方土层失稳。隔水层中的土洞从开始扩展到土洞顶部达到极限厚度是一个逐渐加速破坏的过程，顶板失稳塌陷是土洞扩展的特例。

（2）塌陷过程中的水力特性

在承压水位下降的影响下，岩溶开口处的颗粒优先发生失稳，发生较大运动，从而影响附近的水压分布。如图6 所示为塌陷过程中岩溶开口附近高程水压随时间变化曲线，图6 中y表示为塌陷过程中岩溶开口中心线不同高度（底部y=0 mm）。发生颗粒剥落时，颗粒向下运移，原位置和上方附近水压有陡降，原因来自两方面：一方面是颗粒向下运动，原位置需要得到周围地下水补偿，从而引起周围水压释放；另一方面，颗粒失稳后会在重力作用下加速下降，这个速度一般会比水流速度大得多，会对颗粒附近的局部范围的水流造成扰动，加速水流，从而导致颗粒上方附近的水压降低。这使得上方颗粒附近形成较大的水压差，从而加大上方颗粒塌陷的可能。

图6 塌陷过程中岩溶开口中心线不同高程处水压随时间变化曲线Fig.6 The variation of water pressure at different elevations of karst opening centerline during covered karst ground collapse formation

（3）地下水对颗粒作用力分析

对岩溶开口中心线上7 个颗粒（0377、0588、0799、1854、3331、4808、6285，位置从下至上排序）受到地下水垂向作用力随时间变化作出曲线，如图7 所示，其中0377、1854、3331、4808、6285 五个颗粒在垂向上等距，0377、0588、0799 为垂向上相邻的颗粒。

图7 塌陷过程中岩溶开口中心线颗粒受地下水垂向作用力随时间变化曲线Fig.7 The variation of vertical force on the particles at karst opening centerline during covered karst ground collapse formation

从图7 可以看出，岩溶开口附近的土体颗粒受力失稳剥落后，导致内部颗粒受到地下水的作用力明显增大。例如，颗粒0377 在地下水垂向作用力为-0.073 2 N（负值表示方向向下）时失稳剥落，导致颗粒0588 受到的地下水垂向作用力陡增至-0.144 N 从而失稳，进一步导致颗粒0799 上地下水垂向作用力突变为-0.129 2 N。颗粒一旦开始下降后，受到向下的水流作用力逐步减小。模型表层的颗粒因水流实现贯通，造成所受到的水流作用力紊动。

由此不难推断出，当岩溶开口上方土洞周围存在强度较弱、甚至对土洞内部产生黏性拉力土体时，该部分土体可以在较小水力作用下失稳，但其失稳后会造成较强的瞬间水流作用力，容易使得土洞边沿土体失稳。土洞发展过程中，对土体颗粒的作用力会越来越大，从而导致剥落力呈增大趋势。

（4）位移变化分析

图8 所示为塌陷过程中岩溶开口中心线不同高程处颗粒的垂向位移随时间变化曲线。从图8 中可以看出，岩溶开口中心线上土体颗粒从下向上逐级剥落，高程越高的土体颗粒位移变化较为一致，这说明岩溶开口上方土层从下向上，以颗粒剥落逐步演变为块状剥落，直至盖层整体失稳，整个过程是加速过程。

图8 塌陷过程中岩溶开口中心线颗粒的垂向位移随时间变化曲线Fig.8 The variation of vertical displacement of the particles at karst opening centerline during covered karst ground collapse formation

5 结论

（1）承压水位下降情况下地下水主要对岩溶开口处土层的颗粒产生向下的作用力。在克服土颗粒之间的黏结作用后，岩溶开口处土洞附近的颗粒在重力和水力作用下剥落。土体颗粒的剥落容易造成土颗粒原位置和上方位置处水压的陡降，从而造成较强的水力坡降，使得地下水对内部颗粒作用力陡增，容易进一步引起上方颗粒在地下水作用力和重力作用下失稳，从而使得承压水位下降引起的塌陷从土体颗粒失稳至土层塌陷过程加速。

（2）本文所建立的流固耦合模型对进一步从细观上研究水力驱动的覆盖型岩溶地面塌陷的发育特征以及发育的临界条件、定量分析各影响因子之间的关系具有理论及实际意义。