采空区地表沉降影响因素的研究

王彦海 时文峰 李清泉 陈波

(三峡大学 宜昌443000)

引言

在过去的几十年里，在大量的煤炭资源被开采的同时还造成了大规模、大面积的采空区。 这些采空区极易造成上覆岩体的冒落、弯曲乃至断裂，并且容易使围岩的力学强度降低，故而导致采空区上方的建筑物地基承载力降低，严重的则会造成地表塌陷、大面积的农田损毁、森林植被的破坏、水土流失加剧和土地沙漠化程度的增高等[1，2]，使得矿区居民的财产乃至生命安全受到严重威胁。 针对采空区地表塌陷与变形的研究问题，国内外学者采用了包括概率积分法、岩石力学理论和数值模拟分析等技术和方法进行了一系列的研究[3-5]。 但由于前两者方法的理论公式较为复杂，且计算过程较为冗长，在实际的施工中难以精准测得各项参数，以及其他众多的多变因素是此类方法无法考虑到的。 近年来，随着力学机理研究的深入和计算机技术的发展，数值模拟方法在采空区沉降变形研究中得到了广泛应用。如陈赞成、柴静静、赵付万、张德鹏等学者[6-9]就基于ANSYS 有限元分析软件，对实际工程中的采空区进行仿真计算并取得相应成果，对运用ANSYS 计算采空区的地表沉降做出贡献和指导。但这些研究内容仅为针对某一实际工程进行展开，所得结论不具备一定的普遍性，对于其他的工程而言也不具备一定的参考性。 又如孙超、容宇等学者[10-12]就基于ANSYS 对采空区地表沉降影响因素展开研究，得出了各因素与采空区地表沉降的影响关系。 但上述这些研究都有着考虑因素较少、且没考虑到各因素变化下的相互影响作用、设计工况较少和对土层及采空区的建模较为简略单一等缺点。

对此，本研究选取了数值模拟分析法，通过运用ANSYS 有限元软件对采空区进行数值模拟分析，利用了多种土壤参数模拟出较为复杂的土壤环境，并模拟设计出共计29 种计算工况，对采空区地表变形进行计算分析，得出其地表变形与采深、采宽、采厚以及上覆岩性之间的规律。分别研究了单因素变化下与多因素变化下的各自规律并得出结论，且深入地对采空区地表变形影响因素进行分析，得出了各因素对地表沉降塌陷速度的影响，并且还得出各因素的影响主次关系。 因此对有限元数值模拟方法在采空区地表沉降研究中的应用以及采空区地表的处理和矿区安全生产有一定指导意义。

1 采空区数值模拟分析方案设计

根据在采空区地表变形问题的研究中发现，影响采空区稳定的因素包括矿层的埋深条件、覆岩的力学性质、采空区的几何尺寸等。 其中对采空区地表变形影响尤为显著的有：采掘方式、顶板管理方法、回采率的高低等采空区工程活动[13]。 在工程的实践中，运用数值方法研究一个已经探明且确定的采空区的地表变形，其数值计算模型的建立将会变得非常简单。 但遗憾的是，在工程中对于采空区的几何尺寸不但往往无法探明清楚，而且对于其回采率以及诸多的顶板管理方法也无法精确掌握。 故而想要详细地了解采空区的地下挖开情况是比较困难的[14]。 为了能够便捷地探究采空区地表变形的影响因素，本文研究简化了采空区的实际模型。 将影响采空区地表变形的因素简化为4 个，并均给出4 个不同的水平。 对于多因素多等级的实验工况，其实验工作量巨大，且各个因素之间还可能会有交互作用，相互产生影响。 在研究多因素多等级的实验设计方法中，正交实验设计是一种高效、快速且经济的多因素实验设计方法。 故本研究选取了正交实践设计，运用正交实验表设计出四因素四水平的正交工况共计16 个。 分别计算出各工况下采空区地表沉降数值，并得出各因素对采空区地表沉降的影响趋势函数曲线。 为验证正交实验设计运用在本研究的准确性，又设计出在其余三因素等级不变的情况下，某一单因素等级变化的单因素实验工况共计13 个。 并以单因素实验工况的计算结果为基准验证正交实验工况计算结果的准确性。 以此完成各个因素对采空区地表变形的研究。

1.1 采空区的模型建立

为寻求各因素对采空区地表沉降的影响结果，本研究利用了模拟工况法，选取了采深、采厚、采宽和上覆岩性作为采空区地表沉降的主要的影响因素并将其划分为4 个等级，将其进行汇总见表1。

表1 釆空区地表变形影响因素和等级Tab.1 Factors and grades of surface deformation

其中上覆岩性1 级对应材料折减系数取1.0，上覆岩性2 级对应材料折减系数取0.95，上覆岩性3 级对应材料折减系数取0.9，上覆岩性4 级对应材料折减系数取0.85。

1.采空区地表变形影响分析正交工况组合

为探究各影响因素对于采空区地表沉降的影响，本研究运用正交实验设计研究，将4 个主要影响因素的4 个等级按照正交模型表组合，得出正交工况组合见表2。

表2 多因素正交工况L16(44)Tab.2 Multi-factor orthogonal condition table L16(44)

2.采空区地表变形影响单因素工况组合

为论证多因素正交实验设计的准确性，根据表1 中采空区地表变形因素和等级，将影响因素包括采深、采宽、采厚和上覆岩性4 个因素按4个等级组合起来，在其余因素不变的情况下，研究单因素对采空区地表变形影响规律。 单因素工况组合见表3。

表3 单因素工况组合Tab.3 Single factor working condition combination

为研究采深对采空区地表变形的影响规律，工况 1、2、3、4 保持了采厚，采宽以及上覆岩性不变，将采深逐步增加等级。 为研究采厚对采空区地表变形的影响规律，工况 1、5、6、7 保持采深，采宽以及上覆岩性不变，将采厚逐步增加等级。 为研究采宽对采空区地表变形的影响规律，工况 1、8、9、10 保持采深，采厚以及上覆岩性不变，将采宽逐步增加等级。 为研究上覆岩性对采空区地表变形的影响规律，工况1、11、12、13 保持采深，采厚以及采宽不变，将上覆岩性逐步增加等级。

2 ANSYS模型建立

由于采空区地表的实际沉降问题较为复杂，模型建立难度高，故本文研究将其简化为平面应变问题，进行弹塑性数值模拟分析[15]。 首先用到了solid45 单元建立土壤模型，为了模拟采空区周围岩层的真实应力状态需用到ANSYS 单元生死功能，每次的模拟计算过程分为两步进行：第一步先建完土壤模型后，将土壤模型和采空区模型分别赋予不同的单元号并均保留参数用于模拟开挖前岩层仅有自重作用时的应力状态，以获取模拟计算区域的初始应力状态以及仅由重力而引起的初始位移； 第二步通过ANSYS 单元生死功能，将采空区部分的单元“杀死”，模拟矿体被采出的情况，并计算出此状态下的应力和变形。 最后通过ANSYS 后处理中的荷载步相减求出采空区地表仅由开采矿体引起的真实变形。

ANSYS 计算模型在水平方向取800m，垂直方向取200m，一共分为五层：地表覆盖土层为黄土层，厚度为20m； 第二层为钙质泥岩层，厚度为35m； 第三层为砂质泥岩，厚度为20m； 第四层为中粒砂岩，厚度为20m； 最后一层为粗砂岩，厚度为105m。 采空区的大小和位置由上述工况确定。 由于模型较大，为使计算结果更加精确，对采空区所在的单元采用了ANSYS 网格局部细化功能。 模型底部采用固定边界施加全部约束，模型两侧采用滚轴边界仅施加垂直于该面的约束，取其开采深度为100m、采宽为100m、采厚为2.5m 的工况，并将各材料赋予不同属性并用不同颜色标出的模型示意如图1 所示。 各层的岩土力学参数见表4。

图1 ANSYS 建模示意Fig.1 Schematic diagram of ANSYS modeling

表4 各岩土层及填充材料物理力学参数Tab.4 Physical and mechanical parameters of each rock layer and filling material

3 计算结果与分析

3.1 正交工况结果汇总分析

根据表2 采空区计算模型正交模型表L16(44)运用ANSYS 计算出正交模型下16 个工况的结果并将其汇总见表5。

表5 正交工况L16(44)计算结果Tab.5 Calculation results of orthogonal working conditionL16(44)

将正交工况4 个影响因素中每个因素的4 个等级的4 个工况采空区地表沉降的极差计算值进行汇总，地表沉降正交极差分析见表6。

表6 地表沉降正交极差分析Tab.6 Surface settlement orthogonal range analysis

通过其中极差一栏的数据大小比较，可得出如下结论：

(1)采深因素影响分析：采深的4 个等级其ANSYS 模拟计算得出的采空区地表最大沉降由0.99mm 降低至 0.42mm。 减小比例依次为：11.1%、44.3%、14.3%。 其沉降大小随着采深的增加而明显减小，其原因可能由于采空区顶板及其上部的部分岩层与整体分离，破碎成小块岩块而不规则地填充了采空区，且采空区上方裂缝带上本身可以自成平衡压力拱，采空区越深，其压力拱效果则越显著。

(2)采厚因素影响分析：采厚的4 个等级其ANSYS 模拟计算得出的采空区地表最大沉降由0.34mm 增加至 1.09mm。 增大比例为：50%、64.7%、29.7%。 其原因可能由于随着采厚的增加会引起较大的塌落高度，随之裂缝带、弯曲带的影响范围则会上升，对地表沉降的影响将必然增加。

(3)采宽因素影响分析：采宽的4 个等级其ANSYS 模拟计算得出的采空区地表最大沉降由0.25mm 增加至 1.42mm。 增大比例为：16%、179%、75%。 其原因可能由于随着开采宽度的增加，采空区工作面宽度超过其极限值时，控制层将会因拉应力超过其抗拉强度而断裂，会随着下覆岩层的冒落而下沉，地表沉降会明显增加。

(4)上覆岩性因素影响分析：上覆岩性的4个等级其ANSYS 模拟计算得出的采空区地表最大沉降由0.64mm 增加至0.91mm。 增大比例为：-15%、29.6%、30%。 其原因可能由于随着上覆岩性的折减越大时，对应的岩性越为软弱，在其他条件相同的情况下，坚硬的上覆岩层地表沉降值将小于软弱的上覆岩层地表沉降值。

(5)根据表 6 可知采深、采厚、采宽、上覆岩性四者的极差数据有：1.18 ＞0.74 ＞0.57 ＞0.37。 则在本模拟环境下，影响地表位移的4 个因素中，按照其对地表沉降影响的主次排序依次为：采宽＞采厚＞采深＞上覆岩性。

(6)给出多因素变化下正交工况采空区地表沉降随采深、采厚、采宽以及上覆岩性变化的趋势如图2 所示。

图2 正交实验下各因素各等级沉降值Fig.2 Settlement values of various factors under orthogonal experiment

综上所述，多因素变化下正交实验工况的计算结果完全符合单因素工况的各因素对采空区地表沉降的影响规律。 说明正交实验工况计算数据较为准确。

3.2 单因素工况计算结果汇总分析

通过ANSYS 软件计算出共13 个单因素工况各工况的最大垂直移动、水平移动，见表7。 4个单因素影响下各自变化对地表沉降及水平移动的影响如图3 所示。

表7 单因素工况计算结果汇总Tab.7 Summary of calculation results of single factor working conditions

图3 各单因素变化的影响Fig.3 Effect on various factors

通过上述单因素变化下工况计算结果可得，其结果与正交实验所得结果趋势一致，表明正交实验设计工况的可行性与准确性。 且根据控制变量法下的单因素工况结论，还能得出以下结论：

(1)随着采深的增加采空区地表沉降减少的速率越慢。

(2)随着采厚的增加采空区地表沉降增加的速率越快。

(3)随着采宽的增加采空区地表沉降增加的速率越快。

(4)由于折减系数差值较小，导致上覆岩性对于采空区地表沉降的影响很小，这里难以体现出其对地表沉降变化速率的影响。

4 结论

本文研究运用了有限元分析软件ANSYS 对模拟设计出的正交实验工况进行数值模拟分析，分别构建了各工况下的采空区三维有限元计算模型，且单因素变化下工况组合和多因素变化下的正交工况组合所计算出的结果均表现一致。 对各工况采空区地表沉降以及水平移动的结果进行了比较分析，从而得出以下结论：

1.对于探究采空区的地表稳定性而言，运用数值模拟的方法来计算其地表变形数值是可行的，只要能给出所需的计算参数，即可模拟出采空区地表的各项形变，以预先了解采空区地表变形大小，从而可以为后续工作提供必要的指导。

2.在本模拟环境下，各影响因素对采空区地表变形的影响趋势为：采深对采空区地表沉降的数值呈负相关，且随着采深的增加采空区地表沉降减少的速率越慢。 而采厚、采宽和上覆岩性均对采空区地表沉降的数值呈正相关，且随着采厚与采宽的增加，采空区地表沉降增加的速率越快。

3.在本模拟环境下，通过对正交工况组合的计算结果，建立了极差分析表得出各影响因素对地表沉降影响的主次排序依次为：采宽＞采厚＞采深＞上覆岩性。