王文轩,夏 源

(1.桂林理工大学 环境科学与工程学院;广西 桂林 541004;2.桂林理工大学岩溶地区水污染控制与用水安全保障协同创新中心,广西 桂林 541004)

0 引言

地下水活动导致岩溶塌陷非常普遍,前人研究后总结出几种致塌机制,如潜蚀致塌,真空吸蚀致塌等[1],真空吸蚀致塌是重要的致塌机理之一。随着人类活动的加剧,生产生活中大量抽排地下水[2],导致岩溶空腔中产生气真空值,此时若岩溶空腔上部贯通致岩基面,与覆盖层接触,且覆盖层气密性较好,导致部分处于平衡状态的岩溶土洞发生地面塌陷[3]。这种现象被称为“真空吸蚀”[4]。针对该原理,学者们进行的大量研究,姚春梅等建立岩溶塌陷预警系统模型对岩溶地质灾害进行预警[5];熊启华等[6]结合塌陷成因建立真空吸蚀相关力学模型,并对岩溶塌陷进行分析预测;张鑫等[7]通过理论分析和室内模型试验,定量化的研究水位下降对岩溶塌陷的影响。钟瑶[8]利用FLAC3D数值模拟软件模拟水位下降时岩溶塌陷演化过程,从应力应变角度分析了岩溶塌陷的发育过程;于林弘等[9]通过改变土洞开口大小及水位降幅,对三元覆盖型岩溶塌陷模式进行模拟;方先知等[10]利用有限元软件对岩溶体稳定性进行分析;张晓宸等[11]运用Modflow软件,分析了地下水位变化对岩溶塌陷的影响。虽然学者已经对真空吸蚀导致的岩溶塌陷进行了大量研究,但使用最多的是模拟连续性、均匀性和小变形的FLAC3D、GMS等软件,基于网格方法,可以模拟出塌陷的发生,而后期的大变形破坏则不能模拟,而离散元法可以模拟非连续性、不均匀性和大变形破坏,此类研究较少,本文使用离散元软件MatDEM对岩溶塌陷进行模拟研究。

1 真空吸蚀致岩溶塌陷分析

1.1 覆盖层土体极限平衡公式

在实际的工程中,不同的地下水位变化导致土洞中出现不同真空度,对覆盖层的稳定产生不同的影响。根据文献[12],气压差ΔP一般情况下取值范围在0～50 kPa之间,所以分别选取0、10、20、30、40、50 kPa进行分析,通过不同的气压差模拟不同水位下降深度对岩溶塌陷的影响,根据选取不同的气压差,本实验将真空吸蚀分为六种工况进行模拟分析。

目前,国内对于分析覆盖层土体在岩溶塌陷过程中的受力条件与稳定状态进行了研究,本文通过前人[13]的方法并结合二维离散元模型进行分析,基本假定如下:(1)在真空吸蚀作用下,塌陷体剪切破裂面从拱脚处垂直发展至地表;(2)考虑到覆盖层多为黏性土体孔隙率较小,地下水骤降条件下,外界气体未能及时补给,洞内气压暂时保持不变;(3)地下水位始终位于岩基面以下,土体容重始终不变;(4)在不受外力作用下自然塌陷所形成的土洞体积较小,进行力学分析时忽略此处脱落土体。

如图1所示,土体重度G=γSH;水位下降引起的负压F=ΔPS;覆盖层破坏体侧面摩擦力

图1 覆盖层土体受力分析图

(1)

其中K0为土体侧压力系数;γ为覆土天然容重;H为覆盖层厚度;C为土体黏聚力;φ为内摩擦角;D为覆盖层施加气压的面积;W为二维离散元模型发生侧面位移的周长;为岩溶裂隙宽度,取W=3 m。

由此得到基于二维离散元模型得到的覆盖层土体极限平衡状态时平衡公式:

(2)

式(2)右侧分子部分表示土体极限应力值,分母部分为实际受力情况。为了对岩溶塌陷进行预测评价,将塌陷点的基本数据代入数学表达式计算,当塌陷系数K>1,理论上不会产生岩溶地面塌陷;当塌陷系数K<1,理论上发生岩溶地面塌陷。具体计算结果如表1所示。

表1 六种工况下的塌陷系数K值表

由表1可知,当不考虑其他外界因素时,仅考虑土体的自重压力及气压差时。气压差为0、10、20 kPa覆盖层土体受到的应力达不到极限值,处于稳定状态,岩溶地面塌陷很难发生。当气压差为30、40、50 kPa时覆盖层土体受到的应力达到极限值,发生地面塌陷。

1.2 塌陷临界气压分析

为得到不同覆盖层厚度的塌陷临界气压差ΔP,令式(2)中的K=1,得到覆盖层厚度H与塌陷临界气压差ΔP的关系(图2)。如图3所示,塌陷临界气压差与覆盖层厚度呈正相关,且随着覆盖层厚度的增大临界气压差增大的速率增加;当覆盖层厚度大于9.16 m时,覆盖层土体始终保持稳定,覆盖层厚度小于6.71 m覆盖层土体始终处于失稳状态。

图2 覆盖层厚度与塌陷临界气压的关系

图3 岩溶塌陷离散元模型

2 建立离散元模型

2.1 模型尺寸

如图3所示,岩溶塌陷二维离散元模型长×宽为30 m×12 m,根据矿区水位地质条件,结合数值模拟经验,建立岩层厚度为4 m,覆盖层厚度为8 m,岩溶裂隙宽度为3 m。颗粒单元共37 470个,其中固定单元12 435个(岩土层),活动单元25 035(覆盖层),将岩层单元颗粒设置为固定单元减小计算时间。

2.2 材料参数

设置合适的力学参数在离散元数值建模中设置非常重要,MatDEM离散元随机堆积模型的宏观力学参数与微观力学参数存在解析解,即转换公式[14]。线弹性接触模型的五个微观力学可以通过转换公式,由材料的五个宏观力学性质计算得到

一般情况下,根据转换公式计算得到的紧密堆积模型力学性质小于理论值,为了获取模型颗粒合适的微观参数,通过MatDEM的MatTraining文件实现自动训练材料,经过多次数值测试后得到模型最终的材料参数。根据文献[15]可以得到岩溶覆盖层宏观力学参数值,经过材料训练得到时颗粒单元适合颗粒单元的宏观力学参数值如表2。

表2 覆盖层的力学参数

2.3 数值模拟过程

在MatDEM离散元软件的空箱子中生成一系列具有一定速度的颗粒单元,随后进行重力沉积压实,建立离散元模型;利用MatDEM的切割函数C-tool将模型切割为所需要的岩溶塌陷模型,将岩层设置为墙单元减少计算时间,随后赋材料于覆盖层平衡模型,重新计算模型颗粒的受力状态,使模型重新平衡。

在赋予材料后进行通过过滤矩阵筛选,得到岩溶裂隙上部地表的颗粒单元,根据不同的内外压力差赋予,筛选出来表层颗粒相应的体力;根据不同的真空值对岩溶裂隙正上方的覆盖层施加一个短时间竖直向下的压力来模型岩溶裂隙内外气压差导致的压力,以模拟空蚀对岩溶覆盖层的影响,随后进行自然迭代。如图4所示,在地表颗粒单元上施加竖向体力,模拟真空吸蚀原理对覆盖层土体的作用。

图4 真空吸蚀模型图

3 塌陷数值模拟预测分析

针对六种模拟工况,利用MatDEM进行数值模拟,图5为模拟结果,从中可以看出当气压差为0、10、20 kPa覆盖层土体未发生位移,处于稳定状态,岩溶地面塌陷未发生。当气压差为30、40、50 kPa时覆盖层土体失稳,岩溶地面塌陷发生。数值模拟结果与力学解析法相符。

图5 岩溶塌陷演化过程图(单位:m)

3.1 位移分析

如图5所示,在工况一、二、三情况下岩溶塌陷状况类似,在经过自然塌陷后施加0、10、20 kPa的竖向压力均不能造成覆盖层发生地面塌陷,且根据单元连接图可以看出土体颗粒未发生明显位移,此时覆盖层处于稳定状况,但是施加竖向应力后,岩溶土洞开始发育;在工况四中,覆盖层在气压差为30 kPa的竖向压力作用下,岩溶裂隙处覆盖层向上发育,部分土体脱落。但未形成较大规模,所以覆盖层土体未发生明显竖向位移;在工况四、五、六中均发生地面塌陷。

经过对图5 d-f的单元检测可以得到岩溶地面塌陷的范围(表3)。如表3所示,当发生岩溶地面塌陷时,随着ΔP增加,岩溶地面塌陷宽度增大。

表3 岩溶地面塌陷范围表

3.2 应力分析

由于模拟工况较多,本文选取未发生地面沉降的工况三,和发生地面塌陷的工况五进行应力分析,图6与图7分别是工况三和工况五情况下岩溶覆盖层应力状态。

图6 工况三覆土层应力图(单位:Pa)

图7 工况五覆土层应力图(单位:Pa)

由图6 a可知,覆盖层土体仅发生了内部塌陷,形成岩溶土洞,土洞周围颗粒在水位升降结束后重新排列,应力状态也随之重新分布,在拱效应的作用下土洞上部土体水平应力均大于相同深度土体的水平应力,土洞空腔中的塌落土体的水平应力接近于0,即该处土体未受到形成土洞所产生的水平应力;图6 b为竖直方向应力状态图,整体上竖直方向上应力状态随着土层深度的增加而增大,土洞周围的土体单元经过重新排列和应力重分布后,在洞趾处出现应力集中,土洞上层土体由于拱效应的影响,该处土体的竖向应力明显低于相同深度土体的竖向应力,

由图7 a可知,覆盖层土体失稳,发生地面塌陷,土体水平方向应力整体上随着土层深度的增大而增大,塌落土体经过土体单元重新排列后,水平方向应力重新分布呈拱状,土层表面水平方向上应力较小或者为拉应力状态;图7 b为竖直方向应力状态图,整体上竖直方向上应力状态随着土层深度的增加而增大,塌落土体竖向应力明显小于同一土层深度的其他土体。

4 结语

(1)基于二维离散元模型及前人的研究得到真空吸蚀致岩溶塌陷的解析公式,并利用MatDEM离散元软件进行模拟。离散元数值模拟在模拟真空吸蚀致岩溶塌陷中可以较为准确地预测岩溶塌陷的稳定性,

(2)在覆盖层厚度与岩溶裂隙宽度不变的情况下,发生岩溶地面塌陷的塌陷范围与气压差呈正相关。

(3)岩溶裂隙宽度一定时,塌陷临界气压差ΔP和覆盖层H厚度呈正相关,且随着覆盖层厚度的增大塌陷临界气压差增大的速率增加;当覆盖层厚度H大于或小于一定值后覆盖层稳定不再受气压差影响。