具“上黏下砂”二元结构覆盖层的隐伏岩溶区在我国大中城市广泛分布，如武汉市白沙洲地区、山东泰安市城区等。在此类岩溶区进行工程建设，引发了大量覆盖型岩溶塌陷案例。从物质结构角度而言，都属于砂漏型岩溶塌陷^{[1 − 2]}，即黏性土和岩溶地层之间发育有足够厚度的松散砂性土层时，溶蚀洞隙为砂性土颗粒提供漏失通道、储存以及运移空间，砂颗粒逐渐漏失，形成初步的砂漏空洞。随着砂颗粒漏失量的增加，在自重、荷载等作用下，黏性土层失稳垮塌，在地表形成规模较大的塌陷坑。以武汉市为例，引发砂漏型岩溶塌陷有降雨或江水位波动等自然因素，更多的则由钻探、桩基施工、基坑抽降水等工程活动诱发^{[3 − 6]}。主要诱发机制之一是岩溶水位降低，形成自上而下的渗透力作用于覆盖层砂土，当渗透力或水力梯度超过临界值时漏砂启动^{[7 − 8]}。

针对渗流作用下砂漏型岩溶塌陷机理的研究，较为常用的手段是物理模型试验。例如，较早的有雷明堂等^[9]通过物理模型试验研究武昌区包含二元结构覆盖层在内的多种地层结构类型的临界水力坡降，认为无论是何种地质结构土层，岩溶水位下降是导致覆盖层破坏的主要原因。近年来，吴庆华等^[10]采用恒压取样与剖面压力监测技术，定量研究武汉市砂土覆盖层的岩溶塌陷全过程以验证“渗流液化论”塌陷机理，并研究黏土层对岩溶塌陷的阻隔作用；刘才华等^[11]等利用物理模型试验探讨不同水力梯度和黏砂层厚比条件下二元结构覆盖层塌陷模式。除物理模型试验手段，极限平衡模型和数值模拟也是较为常用的研究手段。例如，熊启华等^[12]针对武汉市长江Ⅰ级阶地“上土下砂”覆盖型岩溶发育区特点，分析提出武汉市长江Ⅰ级阶地岩溶塌陷的成因过程，并结合“潜蚀致塌理论”建立相关极限平衡力学模型，给出土洞形成的临界地下水降幅以及岩溶塌陷发生的数学表达式；屈若枫等^[13]基于PFC（particle flow code）离散元对武汉市4种岩溶地质结构覆盖层塌陷过程进行塌陷形态数值模拟研究，认为上黏下砂地层塌陷机制为下部砂土漏失与上部黏性土弯折拉裂复合型，塌陷坑上部为井状、下部为漏斗状，下部存在呈疏松漏斗状的扰动区。综上可知，目前关于岩溶水位降低导致砂漏型岩溶塌陷机理的研究，多集中于探索临界水力梯度指标获取或塌陷过程宏观演化形态评价，而漏砂启动的宏细观力学机理或力学评价模型还有待进一步研究。

依据陆家街岩溶塌陷实例，开展岩溶水位降低诱发“上黏下砂”二元结构覆盖型岩溶塌陷的物理模型试验与PFC（particle flow code）-CFD（computational fluid dynamics）数值模拟试验，探索塌陷破坏的宏细观力学机理，并构建渗流作用下漏砂启动土拱极限平衡力学评价模型，为二元结构覆盖层岩溶区塌陷防灾监测提供参考依据。

1.   陆家街塌陷实例概况

陆家街岩溶塌陷位于武昌白沙洲，地面平坦，地势低洼，地面高程22～24 m。塌陷区上覆第四系全新统松散堆积物（Qh^al），下伏下三叠统大冶组（T₁d）灰岩，为覆盖型碳酸盐岩分布区。第四系全新统松散堆积物具河流相二元结构特征，总厚度27.46～29.86 m。上部为人工填土和粉质黏土，厚5.65～13.61 m；下部以粉砂及细砂为主，层厚12.84～23.82 m。下三叠统大冶组泥质灰岩及角砾状灰岩，岩层产状倾向北，倾角60°，充填方解石脉并有溶蚀现象，具较多小溶孔及0.1～0.8 m的小溶洞。塌陷区孔隙水位埋深3 m，岩溶水位一般低于孔隙水位，如图1所示。

图  1  陆家街地面塌陷地质剖面图

1—粉质黏土；2—钻孔编号/孔深；3—粉细砂；4—灰岩； 5—溶洞；6—第四系冲积物；7—下三叠统大冶组；8—地下水位

Figure  1.  Geological profile of Lujia Street collapse

下载: 全尺寸图片 幻灯片

根据前期调查资料，陆家街岩溶塌陷前长江水位变化较大，加剧了孔隙水与岩溶水的水位差，促使砂土颗粒向溶隙、溶洞中流失，使土洞迅速扩大，从而发生塌陷。因此，该成因机制主要为孔隙水与岩溶水之间水位差导致砂漏型地面塌陷。

2.   物理模型试验

2.1   物理模型试验装置

图2为模型试验装置，由模型主箱、溶腔模拟箱（副箱）、水位控制系统、监测系统共4部分组成。

图  2  物理模型试验装置及其水位控制系统

Figure  2.  Physical model test device and its water level control system

下载: 全尺寸图片 幻灯片

（1）模型主箱：长1.2 m、宽0.6 m、高0.6 m，分为3个区域，中间填土以模拟覆盖层，左右两侧为水槽，用以隔绝砂土、平稳水流、观测水位。底部开有一圆状开口，并可用球阀控制开口大小（如图3所示阀2），并与下部副箱相连，模拟覆盖层土体从开口发生漏失至副箱。

图  3  模型试验装置原理示意图

Figure  3.  Schematic diagram of the model test device

下载: 全尺寸图片 幻灯片

（2）溶腔模拟箱（副箱）：副箱为整体透明有机玻璃箱。箱顶与主箱岩溶通道相连，用以收集漏失泥沙，判断试验进程。

（3）水位控制系统：由溢流箱与升降机组合而成。其中1#溢流箱与主箱水槽相连用于控制孔隙水位（H₁），2#溢流箱与副箱相连用于控制岩溶水位（H₂），原理示意图如图3所示。

（4）监测系统：用可拆卸的导管将测压管引伸至土层内，用数字孔隙水压力计监测孔隙水压力（ky1）。同时在覆盖层表面和黏砂层交界处布置位移计（wy1—wy4），用以监测沉降变形（图4）。

图  4  模型试验监测布置（单位：cm）

Figure  4.  Monitoring arrangement of model test (unit: cm)

下载: 全尺寸图片 幻灯片

2.2   试验方案与结果分析

根据陆家街塌陷地质条件，可概化为二元结构覆盖层，几何参数为：上部黏性土层6 m，下部粉细砂层厚24 m，孔隙水位埋深3 m。采用相同地层土体作为相似材料，根据模型箱大小，按1∶60几何相似比填筑模型，砂层填筑40 cm，黏土层填筑10 cm，调整1#溢流箱使孔隙水位埋深5 cm，初始时调整2#溢流箱使岩溶水位与孔隙水位一致，从而对应初始水位差为0。保持1#溢流箱不动，调整2#溢流箱使每间隔2 cm逐步降低岩溶水位，每降低一次停留1～2 min观察溶腔中是否漏砂，直到出现漏砂现象或地表塌陷为止。

由图5可知，在水位差达到12 cm之前，岩溶开口中心处的孔隙水压力随岩溶水位降低而有所降低（图5中红色框对应的局部放大图），副箱中有轻微漏砂现象。而当降至12 cm临界水位差时，孔隙水压力下降明显，后续呈波动式下降，且观察到副箱中剧烈漏砂。在漏砂持续约47 min后地面环状裂缝贯通，地面塌陷形成。

图  5  不同水位差条件下物理模型试验孔隙水压力曲线

Figure  5.  Pore water pressure curve of physical model test under seepage conditions with different water level difference

下载: 全尺寸图片 幻灯片

当达到临界水位差后，漏砂启动并逐步向上发展，如图6所示。首先黏砂交界处靠近岩溶开口正上方的沉降曲线wy3缓慢增加，约25 min时沉降量突然增大，对应图5中孔压波动加大，说明此时漏砂量突增。随后wy2开始增加，而wy1几乎保持不变，说明漏砂扩展引起黏土层底板脱空由中心向四周逐步减弱。考虑到黏土层与砂土层变形协调能力的差异，因此wy4开始增大时间晚于wy2和wy3，但由于其监测的是正中心地表沉降，后期沉降量也大于wy2和wy3，最后塌陷形成时wy4也出现陡增现象。

图  6  不同水位差条件下物理模型试验位移监测曲线

Figure  6.  Displacement monitoring curve of physical model test under seepage conditions with different water level difference

下载: 全尺寸图片 幻灯片

当地表出现环状裂缝时说明地面塌陷形成，并停止模型试验。裂缝宽度5～10 mm，环状裂缝呈圆形，直径约39 cm，按1∶60几何相似比进行换算后约为23 m，与陆家街实际塌陷坑形态与影响范围接近，如图7（a）所示。把模型主箱中的水排干后，再开挖揭露的土洞形态为“漏斗状”，是为典型的砂漏型岩溶塌陷变形破坏现象，如图7（b）所示。

图  7  模型试验宏观变形破坏现象

Figure  7.  Model test macroscopic deformation and failure phenomenon

下载: 全尺寸图片 幻灯片

3.   考虑土拱效应的塌陷机理

3.1   PFC-CFD耦合数值模型构建

（1）细观参数标定

细观参数标定可分为力学细观参数标定和流体细观参数标定。

力学细观参数标定采用双轴压缩试验确定，其中砂土模型采用线性接触模型，黏土模型采用平行黏结接触模型。为了模拟土颗粒的不均匀性并同时满足后文以原型尺寸建立的岩溶塌陷数值模型计算效率和精度要求，颗粒单元半径采用均匀分布，颗粒半径范围通过调整后取值见表1。双轴压缩实验分为两步：第一步，对土样施加围压进行预压，使颗粒密实且受力均匀；第二步，轴向加压，进行双轴压缩直至土样数值模型破坏。每一组土样至少在3种不同围压下进行双轴压缩试验，得到应力应变曲线。再根据峰值强度作应力摩尔圆，分析黏聚力与内摩擦角，通过试算调整细观参数至模拟计算的宏观参数与室内试验所得宏观参数一致。

表  1  数值模拟细观参数

Table  1.  Numerical simulation mesoscopic parameters

材料类型	密度/（kg·m−3）	颗粒半径/m	孔隙率	法向刚度/（106 N·m−1）	切向刚度/（106 N·m−1）	法向黏结强度/kPa	切向黏结强度/kPa	摩擦系数	黏滞系数
砂土	2650	0.06~0.30	0.21	20	10	—	—	0.7	—
黏土	2700	0.06~0.30	0.24	20	6.7	40	40	0.5	—
水	1000	—	—	—	—	—	—	—	0.001
注：二维数值模拟孔隙率经验公式：$n_{2\mathrm{D}}=0.42n_{实际}^2+0.25n_{实际}$。

下载: 导出CSV 

| 显示表格

通过监测模型的流体流速与孔隙率，并结合Kozeny-Carman渗透系数计算公式，调整相关流体参数，当测得的渗透系数与室内试验测得的宏观渗透系数相同时，即完成对流体细观参数的标定。

最终得到岩溶塌陷PFC-CFD数值模拟力学和流体细观参数如表1所示。

（2）几何模型

按照案例原型建立长120 m高30 m的模型，按表1中所示颗粒半径范围生成26632颗粒并采用“落雨法”生成模型，如图8所示。上部为黏土层，采用平行黏结接触模型（绿色部分），下部为砂土，采用线性接触模型（蓝色部分），分层设置红色标记，以便观察后续覆盖层变形情况。为平衡效率与精度，选取2.5 m×2.7 m大小的流体单元尺寸，形成12行50列共600个流域单元（包含流域边界单元120个和流域计算单元480个）。力学边界条件为：将岩溶开口设置在模型底部中心位置，为自由边界，模型左右侧及底部（除岩溶开口位置）设置刚性墙体。

图  8  颗粒模型及流体单元划分

Figure  8.  Particle model and fluid unit division

下载: 全尺寸图片 幻灯片

（3）调控流体边界实现渗流模拟

根据陆家街塌陷实例，塌陷处孔隙水位埋深3 m，即孔隙水位27 m。据此，流域网格顶部边界为零孔压边界，左右侧设置静水压力边界，即从顶部到底部压力值从0 kPa线性增加到270 kPa，如图8所示粗网格表示流体边界，包含的颗粒不在流域计算范围内。而在底部岩溶开口处流域网格ID号为（24，0）与（25，0）流体网格设置岩溶水压力边界，底部岩溶开口处流体边界设置分为3个阶段：

①初始无水位差阶段：第一阶段底部流体边界水压力与岩溶开口处流体网格孔隙水压力值统一设置为270 kPa，在自重作用下运行至稳定，获得静水条件下孔压分布见图9（a），不同深度处孔压数值与理论值相符，误差仅0.8%。

图  9  静水与水位差环境下覆盖层渗流场

Figure  9.  Seepage field of overlying layer under hydrostatic and water level difference condition

下载: 全尺寸图片 幻灯片

②岩溶水位逐步降低阶段：减小岩溶开口部位流体网格的孔压值以模拟岩溶水位降低，孔压值根据岩溶水位换算得到，每计算1.2万步岩溶水位降低1 m，直至模型发生漏砂。图9（b）岩溶水位为20 m，即水位差为7 m时对应孔隙水压力分布云图，对应流速分布见图9（c）。

③保持临界水位差致塌阶段：保持岩溶开口部位流体网格水压力值为刚好出现颗粒漏失时对应的岩溶开口处水压力值，直至发生地面塌陷。

3.2   数值模拟结果分析

通过岩溶开口附近单元的孔压监测曲线（图10绿色线）和流速监测曲线（图10蓝色线）可知，在水位差小于8 m情况下，随着水位差增加，孔压呈阶梯形下降，而流速呈阶梯式增加。水位差达到8 m后，砂颗粒开始漏失（图10红色线所示累计颗粒掉落数开始增加），孔压和流速值出现剧烈波动，原因在于颗粒漏失导致监测单元内颗粒迁移碰撞，流体与颗粒力学不平衡所致。另外，如图10所示地表中心处沉降可知，当达到临界水位差时沉降加速。在达到临界水位差之前地表也存在微弱沉降，说明水位差增加导致渗透力增大，颗粒之间变得更加紧密，引起地层沉降。

图  10  孔压、颗粒漏失量、流速、位移数值模拟监测图

Figure  10.  Numerical simulation monitoring results of pore water pressure, particle leakage, flow rate, and displacement

下载: 全尺寸图片 幻灯片

结合物理模型试验结果的临界水位差为12 cm，按1∶60几何相似比换算后为7.2 m，与数值模拟揭示的临界水位差8 m接近，误差约10%，说明该数值模型准确性较高。

静水环境下模型颗粒在重力作用下岩溶开口位置有少许颗粒掉落，但很快达到稳定，形成一个结构稳定的土拱，如图11（a）所示。越靠近岩溶底板，力链越粗，如图11（b）所示；当水位差达8 m时，土洞部位力链相较其他区域更为粗且集中，如图11（d）所示。土洞顶部力链发生断裂，土拱结构开始发生破坏，开始漏砂，如图11（c）所示。由于漏砂导致黏土层底部脱空，黏土层呈简支梁结构，其内力表现为上部为压应力而下部为拉应力（黏土层中红色力链代表拉应力）。保持水位差8 m不变运算至400000步时，如图11（e）所示，颗粒漏失量达到一定程度而使得漏砂扩展范围较大，对应黏土层简支梁结构跨度较大，当黏土层拉应力超过极限抗拉强度时，则发生地面塌陷。此时，黏土层无明显压应力区，漏砂扩展范围内的砂层颗粒力链稀疏，土拱效应消失，说明砂层扰动破坏已发育完成，如图11（f）所示。

图  11  不同演化阶段的覆盖层变形破坏与力链图

Figure  11.  Deformation, failure, and force chain of overburden at different evolutionary stages

下载: 全尺寸图片 幻灯片

3.3   基于土拱极限平衡的漏砂启动分析

砂土中土拱效应广泛存在，最早由太沙基通过活动门试验提出，后续大量学者进一步研究发现其形成有如下3个必要条件：①土颗粒间具有摩擦力；②存在支撑拱脚；③土体内部产生相对位移且应力重分布。粉细砂土颗粒之间存在摩擦力和岩溶开口两侧基岩作为支撑拱脚显然满足前两个条件。根据图5所示试验现象表明在达到临界水位差之前确实存在少量漏砂，说明岩溶开口附近砂层颗粒产生了相对位移，部分颗粒下沉位移较大而漏失，其余砂土颗粒位移较小并相互之间“楔紧”形成土拱，发挥了拱结构承载能力，即满足产生土拱效应的第3个条件。

在岩溶通道上方砂层内形成土拱后，若水位差继续增大，向下的渗透力突破了土拱承载能力后就会引起大量砂颗粒漏失，进而可能引起覆盖层地面塌陷。因此，探究岩溶开口处砂颗粒土拱力学平衡状态对于评价孔隙水与岩溶水的水位差引起渗透力作用下的漏砂启动具有重要意义。

不同于文献关于土洞稳定性极限平衡分析是选择土洞顶板整体进行受力平衡分析，本文建立如图12所示的土拱，取土拱单位体积土体作受力分析。受竖直向下的力有：自重（G）、渗流力（$F_{\mathrm{J}}$）、上覆地层对该单元的压力（$F_{\mathrm{v}}$）。受竖直向上的力有：岩溶水对该单元水压力（$F_{\mathrm{w}}$）和两侧土体对其摩擦力（T）。受水平方向力主要为侧向土压力（$F_{\mathrm{h}}$）。

图  12  土拱单位土体受力示意图

Figure  12.  Schematic diagram of unit soil force of soil arch

下载: 全尺寸图片 幻灯片

土体受力单元采用极限平衡理论定义土拱稳定性系数（$F_{\mathrm{s}}$）如下：

当$F_{\mathrm{s}}$>1时，土拱处于稳定状态；$F_{\mathrm{s}}$=1时，土拱处于临界状态；$F_{\mathrm{s}}$<1时，土拱处于不稳定状态，砂层漏失启动。

其中各分力计算公式如下：

将式（2）—（7）代入式（1）整理后得：

上式中除侧压力系数（$\xi$）和岩溶开口处砂层水力梯度（J）需要讨论确定外，其余参数根据勘察和室内试验资料确定：水的重度（$\gamma_{\mathrm{w}}$）取10 kN/m³、饱和砂层有效内摩擦角（$\varphi$）为32°、黏土层天然重度（${\gamma _1}$）为19.12 kN/m³、黏土层浮重度（${\gamma '_1}$）为7.24 kN/m³、砂层浮重度（${\gamma '_2}$）为8.51 kN/m³、黏土层上部非饱和土体厚度（${h_0}$）为3 m、黏土层饱和土体厚度（${h_1}$）为3 m、砂层厚度（${h_2}$）为24 m，如图4所示。

针对侧压力系数而言，从PFC-CFD数值模型试验中提取岩溶开口处测量圆的竖向应力和水平应力曲线（图13），漏砂启动前（即水位差8 m以前），土拱水平应力、竖向应力、二者之比计算得侧压力系数均随水位差增大而增大，说明随着水位差增大，渗透力随之增大，进而强化土拱应力水平。根据常用的Jaky公式^{[14 − 19]}求解静止侧压力系数（$\xi_静$），即${\xi }_{静}=1-\mathrm{sin}\varphi = 0.47$，与图13中所示砂层土拱中侧压力系数相比都较小，其原因在于土拱效应表现为土体颗粒受左右两侧土体的相互挤压更加明显，水平应力会比正常静止状态下侧向压力偏大。

图  13  应力和侧压力系数随模拟时步变化曲线

Figure  13.  Curves of stress and lateral pressure coefficient changing with simulation time step

下载: 全尺寸图片 幻灯片

根据达西定律，用前文数值模拟获得的流速值（图10蓝色曲线所示）除以砂层饱和渗透系数值$4.3 \times {10^{ - 6}}\;{{\text{m}} / {\text{s}}}$，即可计算出不同水位差条件下岩溶开口处砂层水力梯度值。将侧压力系数和水力梯度随水位差变化情况绘制于图14并线性拟合，将拟合关系代入式（8）便建立了土拱稳定性系数与水位差之间的函数关系如式（9）所示：

图  14  水力梯度、侧压力系数、稳定系数与水位差关系

Figure  14.  Relationship between hydraulic gradient, lateral pressure coefficient, stability coefficient, and water level difference

下载: 全尺寸图片 幻灯片

如式（8）可知，侧压力系数越大代表土拱侧向压力越大，即向上摩擦力越大，抗塌力越大；水力梯度越大代表作用在土拱的渗透力越大，致塌力越大。如图14所示，当把勘察和室内试验数据代入式（9）计算土拱稳定性系数随水位差变化情况，可以看出随着水位差增加，侧压力系数和水力梯度都增大，但水力梯度增大速率明显大于侧压力系数，说明随着水位差的增加，虽然土拱竖向致塌力和抗塌力都在增加，但致塌力增大更明显，因此当水位差增大到8 m时稳定性系数减小到1，与数值模拟结果一致，并与物理模型试验按几何相似比换算后临界水位差7.2 m也接近。

4.   结论

（1）物理模型试验和数值模型试验得到的陆家街塌陷临界水位差分别为7.2 m和8 m，误差约10%。在达到临界水位差后出现剧烈漏砂，砂层底部孔隙水压力骤降，地表沉降增加。

（2）在水位差逐渐增大至临界水位差之前，岩溶开口处砂层流速增大而孔隙水压力减小，土拱水平应力、竖向应力以及侧压力系数增大。

（3）基于土拱极限平衡的漏砂启动分析发现随水位差增加，致塌力和抗塌力都增加，但致塌力增加速率明显大于抗塌力增加速率，当致塌力大于抗塌力时土拱失稳而漏砂启动。

（4）文中所得临界水位差只针对陆家街岩溶塌陷案例，后续可开展不同总层厚、黏砂层厚比、初始水位条件下模拟试验研究。