岩溶地面塌陷稳定性数值模拟分析及防治对策：以深圳市某区域为例

蔡建斯，张洪岩

(深圳市房地产评估和发展研究中心(深圳市地质环境监测中心)，广东 深圳 518034)

1 区域地质背景

深圳市地处广东省南部，东临大亚湾和大鹏湾，西濒珠江口和伶仃洋，南边隔深圳河与香港相联，北部与东莞、惠州两城市接壤。深圳市东部某区域属于丘陵谷地-盆地地貌带，下伏地层主要有石炭系下统大塘阶石磴子组(C1ds)，为浅海台地相碳酸盐岩地层[1]，岩性主要为深灰色白云质灰岩，局部受热液变质或动热变质作用影响，常形成深灰色、灰色、灰白色、白色结晶灰岩、大理岩化灰岩、大理岩等，分布于山间谷地，上部大部分为第四系覆盖，总厚度>700 m。

2 资料收集及分析

本研究收集研究区域钻孔45个，其中遇见溶洞的钻孔11个，遇洞率24.4%，钻孔岩溶率5%～10%，该区域微风化石灰岩岩层表等高线见图1。根据《岩溶工程地质》中岩溶发育程度等级及其指标判别，该区域属于岩溶强烈发育区，需要对该区域进行稳定性分析。

图1 深圳市某区域微风化石灰岩层表等高线图Fig.1 Contour map of the surface of slightly weathered limestone strata in a certain area of Shenzhen city

3 数值模拟分析条件设定

研究区岩溶较发育，其中岩洞埋深、大小、顶板厚度均不等，土洞埋深较深，在30～40 m，例如：ZK010钻孔，溶洞埋深在23.4～33.0 m，顶板厚4.4 m；ZK025钻孔，溶洞埋深在8.5～12.8 m，顶板厚1.8 m；ZK028钻孔，溶洞埋深在19.1～22.5 m，顶板厚3.9 m；ZK034钻孔，溶洞埋深在15.9～29.6 m，顶板厚0.3 m。为对研究区内溶洞稳定性进行基本的判断，根据上述溶洞发育情况，结合城市工程建设动荷载等级，采用数值模拟研究方法对区域内溶洞现状进行稳定性计算及分析[2]。设定模拟条件如下所述。

1) 研究溶洞上方不同级别地表荷载对溶洞稳定性的影响。设定4个荷载级别，分别为100 kPa/m2、300 kPa/m2、500 kPa/m2、700 kPa/m2，模拟得出各级别荷载下溶洞的稳定性状况。

2) 研究溶洞大小对溶洞稳定性的影响。设定4个级别的溶洞规模，分别为4 m×6 m、8 m×12 m、12 m×18 m、16 m×24 m，模拟得出不同大小溶洞的稳定性情况。

3) 研究溶洞顶板厚度对溶洞稳定性的影响。设定两种不同的溶洞顶板厚度，分别为1 m和3 m，模拟分析溶洞顶板厚度对溶洞稳定性的影响情况。

4) 研究溶洞埋深对溶洞稳定性的影响。设定溶洞埋深为14 m、17 m、19 m，模拟分析不同埋深下溶洞自身的稳定状况。

4 模型的建立

溶洞的真实发育情况非常复杂，其几何形态具有不确定性和复杂性[3]，仅凭勘探资料难以对溶洞的实际几何尺寸进行推测。在采用软件FLAC3D进行数值模拟时[4]，如果根据溶洞的实际情况进行建模，那么模型将会非常复杂且不易实现。在实际建模时，根据勘探资料，一方面要尽可能的贴合溶洞发育形态及发育规模[5]；另一方面对溶洞的形状进行简化，用椭圆形来拟合溶洞发育形状，方便模型的建立[6]。

在建模中所用的力学模型，对于土层和岩层都采用摩尔-库伦模型(Mohr-Coulomb)。静力模拟中，所需参数主要包括重力加速度、密度、变形模量、泊松比、内聚力、内摩擦角、抗张强度等。参数主要是根据实测资料和《工程地质手册》。计算体积模量和剪切模量，由变形模量及泊松比求得，其公式见式(1)和式(2)。

(2)

式中：K为体弹性模量；G为剪切模量；E为变形模量；μ为泊松比。

数值模拟的成果分析，主要是根据各种情况下不平衡力曲线、塑性区的分布、最大主应力图、位移等值线图、速度矢量图、位移矢量图等。其中塑性区的分布结果是比较重要的，直接反应溶洞的稳定性和破坏形式[7]；另外，位移矢量图显示溶洞和上部覆盖层土体的移动情况。

5 数值模拟条件选取

5.1 地层概况

根据区域内勘查资料的地层情况，主要为：0～3 m为人工填土；3～10 m为粉质黏土；10～16 m为含砾粉质黏土；16～19 m为微风化灰岩；19～31 m为溶洞；31～50 m为微风化灰岩。

5.2 参数选取

本次模拟参数的选取，主要根据土工试验数据、岩石力学试验以及结合一定的经验值获得，具体参数详见表1。溶洞处于微风化灰岩中，但根据钻孔资料和工程经验，可以确定溶洞周边的岩体裂隙较发育，且有溶蚀现象，因此，对溶洞周边岩体的黏聚力和弹性模量参数进行折减[8]，折减后分别为5.0e5 Pa、2.0e9 Pa，其余参数与微风化灰岩相同。

表1 模拟参数一览表Table 1 List of simulation parameters

5.3 边界条件

根据场地条件，研究区内溶洞有3 m厚的顶板，上部有16 m厚的覆盖层，不考虑地下水的作用。实际模拟面积为80 m×1 m×50 m，并将溶洞简化为12 m×1 m×18 m的椭圆，位于地表下-10 m处，溶洞四周封闭，位于模型中间。模型竖直方向为Z轴，坐标零点位于模型底部中点。数值模拟过程中，主要使用了位移约束边界[9]，约束底部和四周位移(位移为零)，地面为自由边界[10](图2)，边界条件设置如下：位移约束边界，分别约束X=-40 m和X=40 m所决定的面，使之X方向位移为零；Y轴向上为正，向下为负，约束Y=0 m和Y=1 m所决定的面，使之Y位移为零；约束Z=-50 m所决定的面，使Z轴底部位移为零；自由边界，范围Z=0 m的地面。

6 荷载大小对溶洞稳定性的影响

6.1 数值模型

为研究荷载大小对溶洞稳定性的影响，取图2中溶洞尺寸进行分析研究。

6.2 荷载加载工况

为了分析溶洞在不同荷载下的变形破坏情况，本次模拟分析在顶部(地面)上施加了100 kPa/m2、300 kPa/m2、500 kPa/m2、700 kPa/m2四个级别的荷载。

6.3 结果分析

6.3.1 自重条件下的模拟分析

图3为初始应力条件下最大不平衡力与计算时步的关系，在两万多步时达到平衡。模型由开始的最大不平衡力逐渐演化到0，中间有两处突变，这是由于建模过程中改变参数导致[11-12]。

图2 数值模型Fig.2 Numerical model

图3 最大不平衡力Fig.3 Maximum unbalanced force

6.3.2 不同荷载条件下的模拟结果

荷载越大，对下部的影响深度越深。在100 kPa/m2荷载下，溶洞顶板上部的土层发生了将近1 cm的竖向位移；而在700 kPa/m2荷载下，溶洞顶板处的位移发生了近20 cm的竖向位移，如图4所示。荷载施加的越大，在溶洞顶板处出现的拉应力集中区逐渐扩大，在两腰处的压应力集中区的压应力值越大。在荷载不断增大的过程中，塑性区的范围越来越大，在溶洞顶板处的塑性区也不断地增大，但在该种情况下，溶洞顶板处的塑性区并未贯通，溶洞都处于稳定状态。

图4 施加700 kPa/m2荷载条件下塑性破坏区Fig.4 Plastic failure zone under 700 kPa/m2 load

7 溶洞大小对溶洞稳定性的影响

7.1 数值模型

为研究溶洞大小对溶洞稳定性的影响，取轴长4 m×6 m、8 m×12 m、12 m×18 m和16 m×24 m四种溶洞尺寸进行分析研究。

7.2 结果分析

1) 在相同荷载下，各尺寸溶洞模型的竖向位移比较近似，溶洞尺寸对位移值得影响较小。

2) 溶洞尺寸越大，最大主应力范围越大，其中包括最大压应力和最大拉应力。

3) 在300 kPa/m2荷载下，尺寸1的塑性区分布最小，而在尺寸2、尺寸3、尺寸4三种情况下，塑性区的分布范围都较大，但都未贯通。

4) 在500 kPa/m2荷载下，竖向位移和最大主应力的分布与300 kPa/m2荷载下的情况类似。

5) 在500 kPa/m2荷载下，各尺寸下塑性区都较大，有贯通的趋势，但都未贯通，见图5。

图5 施加500 kPa/m2荷载下各尺寸的塑性区分布Fig.5 Distribution of plastic zone of different sizes under 500 kPa/m2 load

8 溶洞顶板厚度对溶洞稳定性的影响

8.1 数值模型

为研究溶洞顶板厚度对溶洞稳定性的影响，取溶洞顶板厚1 m和3 m两种情况进行研究。地层概况、土层参数的取值和边界条件均与第6章节的数值模拟条件选取相同。

8.2 结果分析

1) 顶板越薄，荷载对下部的影响深度越深。

2) 两种情况下，溶洞两侧的压应力集中区的范围大小相差不大，顶板厚1 m的压应力集中区的值要大于3 m厚的情况；而在溶洞顶板处，拉应力集中区主要处于顶板处，顶板越薄，拉应力集中区越小，但拉应力的值却越大。

3) 在顶板厚1 m的情况下，塑性区贯通，而在3 m的情况下塑性区并未贯通，见图6。

9 溶洞埋深对溶洞稳定性的影响

9.1 数值模型

研究不同深度下溶洞的稳定性情况，取14 m、17 m、19 m三种埋深的溶洞进行分析。地层概况、土层参数、边界条件与第6章节数值模拟条件选取相同。

9.2 结果分析

1) 上部荷载对下部的影响深度受溶洞埋深影响，溶洞埋深越深影响深度越深，导致这样的结果主要是受下部岩层厚度的影响。

2) 溶洞埋深越深，顶板处的拉应力集中区值越大，这是因为溶洞埋深越深受上部覆盖层的竖向压力越大。

3) 溶洞埋深越浅受上部荷载影响越大，在埋深14 m的情况下塑性区已经贯通，见图7。

图6 不同顶板厚度在500 kPa/m2荷载下塑性区分布Fig.6 Distribution of plastic zone of different roof thickness under 500 kPa/m2 load

图7 不同埋深在500 kPa/m2荷载下塑性区分布Fig.7 Distribution of plastic zone under 500 kPa/m2 load with different depth of burial

10 结果分析

1) 上部荷载对溶洞的稳定性起着至关重要的作用。荷载越大，影响深度越深，对溶洞的稳定性影响越大。

2) 溶洞尺寸的大小对溶洞稳定性也有影响的。溶洞越小越稳定，溶洞越大越不稳定，但尺寸大小对溶洞稳定性的影响较小。

3) 溶洞顶板对溶洞的稳定性起着决定性的作用。溶洞顶板厚度越薄溶洞稳定性越差，溶洞顶板厚度越厚溶洞稳定性越好。

4) 发生岩溶地面塌陷时，首先是顶板出现塑性破坏，拉张破坏，并沿着与溶洞直径近似大小的范围内向上延伸，直到与上部土体发生塑性破坏，整个塑性区贯通，发生塌陷，这与岩溶地面塌陷多以圆形或椭圆形的形状相符。

5) 溶洞的埋深决定着溶洞的稳定性，溶洞埋深越深受上部荷载影响越小。

6) 本次模拟分析只考虑最安全模式，即在地表施加荷载的模式下进行，未考虑工程建设过程中的基坑开挖情况。如果进行基坑开挖，将会直接减小上部覆盖层厚度和增加侧向动荷载，岩溶地面塌陷灾害更易发生。

11 防治对策

11.1 保持岩溶水系统中输入和输出的平衡状态

严格控制地下水的开采，汛期有针对性的控制地下水的补给输入，避免岩溶水系统中局部地段的输入和输出量的急剧变化，导致改变水力坡度和流速，加剧水流的冲刷侵蚀作用和对溶洞内填充物的冲蚀及迁移作用，有效减缓或防治塌陷的出现。

11.2 避免各种附加荷载破坏土体的极限平衡状态

在岩溶的高发育区，要避免各种地面动荷载破坏其极限平衡状态[13]。在工程建设开工前，采取一定的工程措施提高该区域的稳定性，如回灌填充或加固等方式后，再进行工程建设。

11.3 防止地表水的集中渗入

定期检查、维护地表水系统的防渗漏措施，防止地表水集中、大量渗漏所产生的动能破坏土层的稳定状态，以免对塌陷区造成进一步的影响[14]。

11.4 建立长期监测点

在岩溶高发育区内，选取人口密集区，如医院、学校、集市等区域设立监测点，条件允许时设立自动监测系统，长期监测地下水的变化幅度及地面沉降的异常情况。一旦发生异常情况，及时进行预警预报。