固体充填开采地表沉陷数值模拟研究

2015-12-04康育鹏孟海梅

河南城建学院学报 2015年1期

康育鹏，孟海梅

(河南理工大学测绘与国土信息工程学院，河南焦作454003)

煤炭是我国的主要能源，在我国的一次能源消耗中占65%左右，我国的煤炭资源虽然丰富，但经过长期大量开采正在逐年减少，尤其是在我国的东部地区，煤炭资源正在逐步枯竭，剩余储量的一半以上属于“三下”(建筑物下、水体下、铁路公路下)压煤。据不完全统计，全国生产矿井“三下”压煤量达137.9亿t，其中建筑物下压煤87.6亿t，村庄下压煤占建筑物下压煤的60%以上［1］。煤矿“三下”压煤问题是当前和今后一段时期内煤炭行业发展面临的主要问题之一，关系到煤炭资源的合理开发与矿井的可持续发展。

固体充填开采是近些年发展起来的一种绿色开采技术，它是利用矸石等固体充填物充入采空区，占据上覆岩层的下沉空间，以阻止顶板的下沉［2－3］。相比于传统的垮落法开采，它可以有效减小采动损害，控制上覆岩层与地表的移动变形，大大减少矸石等固体废弃物的排放，进而解决煤矿的“三下”压煤问题。

由于使用FLAC3D等数值模拟软件进行开采沉陷的模拟与概率积分法得到的预计结果比较相近，所以FLAC3D等数值模拟软件可以被广泛应用于开采沉陷的预计与可视化方面［4］。固体充填开采虽然可以有效减小采动损害，但仍不可避免地造成一定程度的地表沉陷，因此需要对固体充填开采所引起的地表沉陷进行预测。本文尝试使用FLAC3D软件对固体充填开采进行数值模拟研究，并应用其内嵌的FISH程序语言，编制相关程序以提取出固体充填开采地表的沉陷数据，结合MATLAB强大的数据可视化功能完成对固体充填开采地表沉陷的可视化操作。

1 数值模拟与模型建立

1.1 数值模拟及FLAC3D简介

数值模拟也叫计算机模拟，它是通过数值计算和图像显示的方法，借助计算机来获得满足工程要求的数值解，对工程问题和自然界各类问题进行研究。数值模拟现在已经成为矿山开采沉陷研究的重要手段之一。

FLAC3D是三维快速拉格朗日分析(Fast Lagrangian Analysis of Continua)的简称，FLAC3D程序自美国ITASCA咨询集团推出后，已经成为目前岩土力学计算中的重要数值方法之一。该程序是FLAC二维计算程序在三维空间的扩展，用于模拟三维土体、岩体或其他地质材料达到强度极限或屈服极限时发生的破坏或塑性流动的力学特征，特别适用于分析渐进破坏失稳以及模拟大变形［5］。FLAC3D是采用显式算法来获得模型全部运动方程的时间步长解，从而可以追踪材料的渐进破坏和垮落，这对研究开采的时间效应和空间效应是十分有效的。此外，该程序提供了多种模拟材料的本构模型，可以输入和改变整体或局部的模拟材料类型和参数，增强了程序使用的灵活性，用来提供采动区域的垮落过程和开采中的充填过程［6］。FLAC3D拥有自己的内嵌程序语言FISH，使用户能够自定义变量和函数，扩大了FLAC3D的应用及用户自有的特色。FLAC3D还具备比较完善的后处理系统，用户可以提取和查看所需要的计算信息，并且可以生成和输出图形、切片、动画等多种形式的后处理结果。

1.2 地质采矿条件

以某矿区地质采矿条件为模型，模拟开采区域位于某村庄下，所开采的煤层为5煤层，煤层开采的平均厚度为5.0 m，平均采深为695 m，煤层倾角为0°～2°。由于煤层的倾角很小，可以将该煤层看作水平煤层。煤层的上覆岩层以砂岩和砂质泥岩为主，其中第四系松散层厚度约为186 m。因为煤层开采的区域位于村庄下，而该村庄建筑物较为密集并且抗变形能力不强，容易受到采动影响而发生破坏，所以采用固体充填开采的方法进行开采，并且在开采之前需要进行开采沉陷预计。各岩层的物理力学参数如表1所示。

1.3 数值模型建立

假设试采工作面长度为400 m，工作面推进距离为1 000 m。工作面采用走向长壁固体充填开采的方法，充填为即采即充，充填率为100%。结合矿区的地质采矿条件和试采工作面的尺寸等相关信息，运用FLAC3D数值模拟软件建立大小为3 000 m×2 400 m×715 m的三维模型，该模型共有24 800个矩形网格单元和261 792个单元节点。三维数值模型如图1所示。

图1 三维数值模型

2 等效采厚与模拟开挖

2.1 固体充填开采的等效采厚

固体充填开采之所以能够减小岩层及地表的移动变形，就是因为矸石等固体充填物占据了上覆岩体的下沉空间，阻止了顶板的垮落式下沉，使得顶板岩层最终下沉较小，这就相当于减小了煤层的开采厚度［7］。固体充填开采的等效开采厚度就是充填工作面煤层顶板从开始移动到最终稳定时所移动的距离。

根据固体充填开采的工艺过程以及顶板岩层的移动下沉机理可知，等效采厚与充填前顶板的移近量、充填的未接顶量、充填率以及固体充填体的压缩率密切相关。固体充填开采的等效采厚的计算公式为:

式中:Me为煤层的等效开采厚度;S1为充填前顶板的移近量;S2为充填的未接顶量;S3为充填体最终的压缩量。

其中S1、S2和S3的计算公式为:

式中:M为煤层的实际开采厚度;η为充填前顶板的下沉率(%);ρ为充填率(%);μ为充填体的压缩率(%)。

所以固体充填开采的等效采厚的计算公式可写为:

由于固体充填开采所引起的岩层与地表移动变形是充填工作面的等效采厚引起的，而固体充填开采的等效采厚一般都明显小于煤层的实际采厚，所以固体充填开采可以有效控制岩层与地表的移动变形。在对固体充填开采所引起的地表沉陷进行预计时，可以把它看作是开采厚度为等效采厚Me的煤层在垮落法开采条件下的情况。沉陷预计的方法可以采用概率积分法，数值模拟分析法等。

2.2 模型的模拟开挖

煤层的实际开采厚度M=5 m，充填工作面的尺寸为1 000 m×400 m。假定充填前顶板的下沉率η=5%，充填率ρ=100%，固体充填体的压缩率μ=10%。结合固体充填开采的等效采厚计算公式，得到的等效采厚Me=0.725 m。

在模拟工作面开挖时，只需要将所开挖的区域定义为空模型(Null Model)。而在模拟工作面充填开采时，需要给之前开挖的部分区域赋予新的材料本构模型，如摩尔－库仑模型，双屈服模型等，并设定对应模型的材料参数和初始条件，然后进行解算即可。模拟固体充填开采的FLAC3D程序为:

在以上的FLAC3D模拟开挖程序中，可以在solve这一行命令之前加上如“hist id=1 unbal”和“hist id=2 gp zdisp 1500，1200，715 ”这样的命令，用于采样最大不平衡力和节点(1500，1200，715)在 Z方向上的位移随时间步的迭代变化，并可使用FLAC3D的plot绘图命令显示出来。如果需要采样的记录较多时，用户可以输入“print history”命令，显示出什么样的记录号对应什么变量或参数。在解算完成后，用户可以通过输入诸如“print gp disp id 257702”或“print gp disp range x 1500 y 1200 z 714 715”这些命令得到id=257702或range命令包含的相应区域内的节点位移。

在模型的模拟开挖解算完成后，可以运用apply命令在模型上边界的部分区域施加一定的竖向应力，用来模拟建筑物荷载作用下的情况。比如在开挖后模型的上边界中央区域内施加100 kPa的一般建筑工程荷载后，就可以得到模型的上边界增大了约1.5 cm的竖向位移。考虑到一般村庄建筑物的荷载较小(＜100 kPa)以及岩土体的固结变形等因素，可以认为村庄建筑物自身荷载对于煤矿开采沉陷区域内地表的移动变形影响很小，某些情况下甚至可以忽略不计。

由于FLAC3D具有比较强大的后处理功能，在模型的模拟开挖和施加建筑物荷载等解算结束之后，用户就可以使用一系列的plot绘图命令将所需要的图形显示出来。分析云图如图2～图6所示。

图2 模型的沉降位移云图

图3 模型的最大位移云图

图4 走向观测线沉降位移云图

图5 倾斜观测线沉降位移云图

3 地表沉陷数据提取与可视化

3.1 观测站布设及地表沉陷数据的提取

为了全面、准确地反映固体充填开采对地表移动变形的影响，需要在三维数值模型上模拟建立等间距的剖面线状观测站或网状观测站。由于煤层的平均采深为695 m，根据表2所示的开采深度与测点间距的关系，设定测点间距为30 m［8］。因为FLAC3D软件和其内嵌的FISH程序语言只能查询或输出网格节点和单元体的空间位置等相关信息，所以在先前建立三维数值模型时，就将各个层组的网格划分为长宽均为30 m的单元格，以便对测点间距为30 m的各地表移动观测站进行沉陷数据的提取。

当模型上建立的移动观测站是剖面线状观测站时，就仅需要提取走向观测线和倾斜观测线上各观测站上的沉陷数据。假定走向观测线的Y坐标为1 200，倾斜观测线的X坐标为1 500，利用FISH语言编制以下程序:

图6 施加100 kPa建筑物荷载后的模型沉降位移云图

表2 地表移动观测站测点密度表

在模型开挖的解算完成后，运行以上程序，就可以在xyw.log文件中得到走向观测线和倾斜观测线上各观测站的下沉值。

当模型上建立的移动观测站是网状观测站时，则需要提取地表上开采影响区域内所有间隔为30 m的各移动观测站的沉陷数据。由于FISH语言中的table函数只能产生含有X、Y这两个字段的表，因此在提取网状观测站的沉陷数据时，就不能利用剖面线状观测站提取沉陷数据时采用的table、xtable、ytable等这样的函数。针对这一问题，可以通过提取模型开挖后地表上间距为30 m的各节点的三维坐标和位移来解决，对应的FISH程序为:

运行以上程序，就可以得到模型开挖后地表上各节点(移动观测站)的三维坐标和位移。通过选择xyz.log文本中的X、Y坐标和对应节点Z方向上的位移(Z－Dis)，就可以整理得到地表上各移动观测站的沉陷数据。

通过提取各个地表移动观测站的沉陷数据，得到移动观测站的最大沉降值为56 mm，各观测站之间的最大倾斜值为0.102 mm/m。即使在施加了100kPa的一般建筑工程荷载后，移动观测站的最大沉降值仅为72 mm，各观测站之间的最大倾斜值也仅为0.232 mm/m。结合我国一般砖木结构建筑的临界变形值［8］和《建筑地基基础设计规范》［9］(GB50007－2011)中所规定的建筑物地基的最大沉降量(见表3)可知，各观测站的最大沉降量和观测站之间的最大倾斜值都远远小于规定的最大沉降量和临界变形值。因此采用固体充填开采的方法可以有效减轻或消除煤炭开采对于地表建筑物的影响，适用于村庄等建筑物下采煤的情况。

表3 建筑物地基的最大沉降量和地表的临界变形值表

3.2 利用MATLAB实现后处理与可视化操作

MATLAB是美国MathWorks公司研发的目前世界上使用最广泛的科学计算软件之一，它具有矩阵计算、数值分析、数据可视化以及交互式程序设计等强大功能［10］。在这里，本文利用MATLAB强大的数据可视化功能，结合提取的地表移动观测站的沉陷数据，运用MATLAB的cftool曲线拟合工具箱以及plot3、mesh、contour、clabel等MATLAB绘图函数，实现固体充填开采地表沉陷的可视化，绘制并输出相关的图形(见图7～图10)。