煤层开采后覆岩稳定性物理模型及数值模拟

2019-04-09张谦张延军于子望张通张雨罗银飞

世界地质 2019年1期

张谦，张延军，于子望，张通，张雨，罗银飞

1.吉林大学建设工程学院, 长春130026；2.青海省环境地质重点实验室, 西宁810007；3.青海省环境地质勘查局, 西宁810007；4.青海省地质环境保护与灾害防治工程技术研究中心, 西宁810007

0 引言

煤炭是中国重要能源之一，与国民经济的发展有着密切联系，随着中国经济的不断发展，煤炭作为主要能源其消耗量也在与日俱增[1,2]。煤层开采过程中引发的上覆岩层移动及地表沉陷等问题严重威胁着人民的生活和生产。因此对煤层开采后上覆岩层变化规律的研究具有重大意义。现阶段，对上覆岩层变化规律的研究方法主要为相似材料模拟和数值模拟。由于岩土介质本构的复杂性及结构与岩土介质之间的相互作用，相似材料模拟的模拟结果相对贴近实际，一直是研究覆岩稳定性的重要手段[3--6]，但相似材料模拟存在着实验费用高，场地约束大等问题[7--9]。采用数值模拟具有可以灵活地改变模型的形状及参数，受客观条件影响小，节约经费且能处理复杂的材料本构关系等优点，技术上与相似材料模拟优势互补，结果上相互印证。

随着FLAC软件的不断成熟，近年来其已经成为众多学者研究地表沉降规律的一种重要手段。李培现等[10]用FLAC3D对老采空区的地基稳定性进行了研究。郭俊廷等[11]用 FLAC3D对非充分采动下地表移动变形值进行分析，并获取了地表的移动参数。李想等[12]用 FLAC3D软件求解地表移动带范围，优化了圈定移动带的方法。周振亮等[13]用FLAC求解出条带开采对地表变形值的影响。何标庆[14]采用FLAC软件对大型采空区群的稳定性进行模拟。一些学者[15，16]将FLAC模拟地表移动的结果与理论公式计算结果及INSAR监测数据进行对比，进一步证明了FLAC模拟的准确性。笔者通过建立物理模型和采用FLAC3D软件进行数值模拟两种方法，将鄂尔多斯某煤矿的物理模型开挖结果与数值模拟结果对比分析，得出该煤矿下部煤层开采对上覆岩层稳定性的影响范围及规律。

1 地质概况及模型建立

研究区位于鄂尔多斯高原东北部，地貌总体上波状起伏，河流不甚发育，属于典型的构造剥蚀高原地貌。地势开阔呈波状起伏，枝状沟谷发育其中；南部地势平坦，多被风积砂覆盖，并有洼地、湖盆散列其中。

1.1 物理模型建立

在几何相似的多个系统中，进行相同性质的物理过程，当几何上的对应点在相应时间点上的表现具有一定的一致性时，称这样的物理过程为相似现象。具有相似现象的物理系统在几何形态上具有相似性，遵循着同样的物理定律。

相似系数是原型与模型之间相同物理量之比，通常用C表示。在煤层相似材料模拟中常用的相似系数有：

几何相似系数：

(1)

式中：δp为原型位移；δm为模型位移；为Lp原型长度；Lm为模型长度。

应力相似系数。

(2)

式中：σp为原型应力；σm为模型应力。

容重相似系数：

(3)

式中：λp为原型容重；λm为模型容重。

为准确建立研究区模型，收集了研究区内的工程地质、环境地质、水文地质资料，根据研究区内的钻孔资料，以及参考前人的资料[17]，最终确定物理模型的相似常数为：CL=200、Cλ=1.5、Cσ=300。原地层平均容重为λ1=26.32 KN/m3，因此模型的相似配比材料容重设计为λ2=17 KN/m3。由于原型横断面高290 m，最终模型整体尺寸设计为长×高×宽=1.2 m×1.45 m×0.12 m。

为简化模型铺设，将部分岩层简化为均布荷载铺设于模型的上部，由于上覆岩层的岩性变化引起下部岩层荷载分布的动态与非均布的特点，模型的上边界加载遵循的原则：岩性均一的覆岩可简化为均布荷载加载于模型的上边界，覆岩中有突变岩层时，只有突变岩层的上部岩层可以简化为均布荷载。按照此原则最终将原型简化为十层模型。模型由上至下分布如表1。

根据相似材料性质、模型特点、实际地质情况及前人对采空区模型实验的经验[18]，物理模型采用石英砂、细砂做骨料，石灰、石膏作为胶结材料，由下至上逐层建立模型，通过控制胶结材料的使用量控制岩层强度，以每次2～3 cm的厚度进行模型装填，由于石膏具有遇水极易凝结硬化的特点，加入浓度为1.0%的硼砂溶液作为缓凝剂，并在每层之间加入云母粉模拟实际岩层之间的原生裂隙，模型制作完成后进行一周的晾置，得到最终模型。

表1 相似模型地层分布Table 1 Stratigraphic distribution of similar model

模型建立完成后，采用在顶层表面布置千分表的方式对实验数据进行监测。在模型最上方非等距布设千分表，靠近模型中轴处，千分表的布设间距较小，随着远离中轴，千分表布设距离增大。开采方法采用水平方向壁式采煤法，开采实际宽度200 m(模拟开采宽度100 cm)的煤层。

图1 物理模型图Fig.1 Physical model diagram

1.2 数值模型建立

通过研究区地质数据统计资料，确定地层标高及厚度；参考钻孔资料确定煤层埋深，最终建立地质概念模型。通过工作区内钻孔物理、力学测量结果，结合实验室力学实验最终确定数值模型参数。为方便模型建立与计算，将物理性质相似的非关键层简化为同一层，采用Midas描述并建立研究区数值模型(图2)，参数设置如表2。

表2 数值模型各层参数Table 2 Layer parameters of numerical model

图2 模型建立及网格剖分Fig.2 Modeling and meshing

2 模拟结果及分析

2.1 物理模型模拟结果及分析

实际生产中采用壁式开挖法对3--1煤层进行开采，开挖煤层实际厚度为4 m(模型厚度2 cm)，实际开采宽度200 m(模型宽度100 cm)，实际开采时间为5 a。实验中，由于模型开挖层厚度较小(2 cm)，因此采用小刀进行开挖，将距模型右边界20 cm处设为右端点，在右端点处以10 cm/h的速度由右向左对3--1煤层进行开挖，一次开挖时长为2 h。每次开挖完成后，每隔1 h测量并记录冒落带和裂隙带高度，记录千分表监测数据，验证模型是否达到稳定，若模型达到稳定，则进行下一次开挖。

图3 物理模型开挖效果图Fig.3 Physical model excavation effect diagram

实验过程中，第三次开挖，即开挖长度为60 cm时，可以看到模型出现裂隙，在第三次开挖过程中伴随部分坍塌与脱落现象的出现(图3)。从采动岩体移动和离层裂隙发育情况可以看出，随着工作面的推进，岩体裂隙自上而下逐步发展。在采动作用下，距离开挖层较近的层面首先开裂，随着开采范围的增大，距离开挖层较远的岩层也开始产生裂隙。随着工作面的进一步推进，在前一工作面推进距离条件下的采动岩体裂隙网格上又叠加了后一工作面推进距离条件下的采动岩体裂隙网络，从而使采动岩体裂隙分布更加复杂化，不仅有新的裂隙产生，而且原有裂隙也发生了一系列的扩展、闭合以及张开。

实验发现，随着开挖的进行，模型的裂隙发育率逐渐增大，开挖0.5 m时，模型的裂隙发育率为23.11%，当开挖至1.0 m时，模型的裂隙发育率为30.41%。经测量，裂隙带高度为32.5 cm，约为采高的16.25倍。煤层采出后，冒落带仅出现在开采区域，冒落块体以外的岩体不受影响。当工作面开采至1.0 m时(实际200 m)，覆岩冒落高度为7.5 cm(实际15 m)，是采出煤层厚度的3.75倍。冒落带的左侧垮落角36°，右侧垮落角57°。实验开采完成时，裂隙并未延伸到含水砂层，但是考虑地层原有裂隙的存在，3--1煤层开挖还是会对上方的含水砂层产生较大的影响。

图4 冒落带、裂隙带、底板位移图Fig.4 Displacement diagram of caving zone, fissure belt and floor slab

根据地表千分表制作底板位移图，采用网格测量法，制作冒落带、裂隙带图。由图4可以看出，底板的最大位移近似等于采高，最大位移为2 cm，裂隙带位移远小于底板位移，最大位移为0.9 cm，冒落带最大位移为0.51 cm。以上3个层均具中间位移大，两端位移小的特征，近似呈以采空区中心线为对称轴的类似“V”字型。

2.2 数值模型模拟结果及分析

2.2.1 摩尔库伦本构模型下数值模拟分析

数值模拟的主要目的是计算煤层开挖对上覆岩层造成的影响，首先应消除重力作用对上覆岩层的影响。因此，在开挖前，需要对初始应力进行求解，将初始应力解出的位移和应力值清零后再进行开挖，以此确定开挖造成的地表形变。

达到初始应力过程中的最大不平衡力的变化趋势如图5所示，从图5中可以看出，最大不平衡力先增加后逐渐变小，最后趋于稳定，最终最大不平衡力数值<1.0×10-5，说明模型已经达到初始平衡。

图5 最大不平衡力变化图Fig.5 Variation diagram of maximum unbalanced force

根据研究区勘察报告中的采矿区开采情况，本次数值模拟采用走向长壁式开挖法对3--1煤层进行开挖，设计年采煤量288万t，平均采高2.02 m，工作面长度250 m，年进尺400 m，开挖走向由南向北开挖，计算五年内采空区对上覆岩层稳定性的影响。

采用FLAC3D对采空区开挖模拟的竖向位移如图6所示。由图6可以看出，煤层开采后，煤层上方岩层逐渐悬空，地下形成采空区，此时顶板岩层可能会出现冒落、坍塌的现象，远离采空区，这种现象也逐渐变弱。岩层的弯曲、破环现象一般沿着岩层的层理法线方向发生，伴随着这些现象发生的还有部分离层现象与断裂的产生。随着开采的进行，地层的表面也开始出现这些现象，此时地表形成下沉盆地。由图6可以看出，开挖对地表的影响范围要大于采空区实际尺寸，地表下沉形状近似椭圆状。并且这种现象随着时间的推移，变得愈发明显，在第五年时，地表沉降及位移达到了最大值。地表最大沉降量也随时间推移逐渐增大，第一年地表最大沉降量为0.093 m，第二年为0.301 m,第三年为0.500 m，第四年为0.676 m，第五年为0.833 m。

由于自重应力的存在，在采煤过程中，岩层的初始平衡状态遭到破坏，应力重新调整、分布，最终达到新的平衡状态。FLAC3D应力分布云图通过颜色深浅显示出工作区的各个区域位移值的大小关系，越趋近于红色的区域，应力集中现象越明显。由剪应力分布图与垂向应力分布图中可以看出，在开采的第一年，上覆岩层中出现应力集中现象，随着煤层开采的推进，在采空区及覆岩所在的区域均有剪应力形成，且伴随着开采面积的增大，剪应力的范围与应力集中现象也逐渐地增大，在第五年，最大剪应力已经达到3.63×106Pa。随着上覆岩层应力集中现象愈发明显，上覆岩层逐渐破坏，当应力增加到一定程度，上覆岩层出现裂隙，导致岩层失去隔水性，最终发生破坏。

对比五年塑性区变化可以看出，当开采面积增大，应力随之增大，塑性区的范围也逐年增加。从垂向应力可以看出，由于开挖，地下应力重分布，采空区及上覆岩层周围出现应力集中现象，在采空区两侧顶点处这种现象尤为明显。由应力状态可以推出，在第五年时，由于开挖作用，上覆岩层底板在重力及覆岩的共同作用下，在法线方向产生弯曲与移动，变形达到一定程度后，顶板内部应力剧增，应力增大到超过岩石的极限抗拉强度时，岩石发生破坏，塑性区出现，上覆岩层也随之出现断裂、破碎。

2.2.2 修正剑桥模型下3--1煤层开挖模拟

采用修正剑桥本构模型，对研究区开展数值模拟。建立模型的尺寸、地层分层及开挖方式均与摩尔库伦本构模型相同，同样选择开挖3--1煤层，工作面长度为250 m，年进尺400 m，开挖走向为由南向北，其模拟结果如图7、8所示。

从图9可以看出，采用修正剑桥模型和库伦模型模拟得到的位移分布形态基本相似，但沉降量模拟结果上，二者差异较大，修正剑桥模型模拟至第五年时，地表最大沉降量为0.377 m。出现这种结果的原因可能是修正剑桥模型主要适用于土层的模拟而摩尔库伦模型主要适用于岩层的模拟。

图6 3--1煤层开采位移分布云图Fig.6 Displacement distribution cloud map of 3--1 coal seam

图7 第五年I-I剖面塑性区分布剖面图Fig.7 Distribution profile of plastic zone in I-I profile in the fifth year

图8 第五年I-I剖面剪应力以及垂向应力分布图Fig.8 Shear stress and vertical stress distribution map of I-I section in the fifth year

图9 修正剑桥模型下开挖3--1煤层第五年位移分布云图Fig.9 Correction of the fifth year displacement distribution of 3--1 coal seam under Cambridge model

3 理论公式验证

为进一步验证数值模拟的准确性，将模拟的围岩压力结果与理论公式中的普式压力拱理论计算结果进行横向对比。普式理论认为，松散岩体中作用于硐顶的围岩压力由自重应力提供，工程中为了方便计算，往往将硐顶的最大围岩压力简化为均布荷载而不计硐轴线的变化引起的围岩应力变化。因此，硐顶的最大围岩压力公式：

(4)

侧向压力公式：

(5)

式中：f为坚固性系数；b、b1为拱的矢高；a为侧壁稳定时平衡拱的跨度；a1为自然平衡拱的最大跨度。

图10 煤层开挖中轴线上覆岩层垂向最大应力对比图Fig.10 Vertical maximum stress diagram of overburden strata on axis of coal seam excavation

对比普式压力拱理论计算结果与摩尔库伦模型数值模拟计算结果发现，两者计算结果中最大竖向应力的变化数量级为105Pa。两者的模拟曲线整体拟合效果较好，对于形态上的差距，这可能是由于普式压力拱理论在计算等速开采时，计算结果只受开采宽度的影响，因此其计算的硐顶最大竖向压力表现为线性增加。而数值模拟计算的结果由于考虑了材料的弹塑性变形以及围岩摩擦力的影响，其计算的最大竖向应力曲线则表现出一定的缓变性。两者计算结果的对比进一步证明了采用摩尔库伦模型进行数值模拟的结果更加准确。