坝下深部煤层开采坝体移动变形规律的数值模拟

2022-04-22秦洪岩题正义

煤矿安全 2022年4期

张峰，秦洪岩，题正义

（1.山西工程技术学院采矿工程系，山西阳泉 045000；2.辽宁工程技术大学矿业学院，辽宁阜新 123000；3.华北科技学院，北京 101601）

库坝下深部煤层开采，既要保证煤层开采的安全性，也要保证坝体的安全运营，掌握煤层开采过程中坝体的移动变形机理，才能最大限度地对库坝下煤炭资源进行合理开发利用[1]。但煤层埋藏深度大、地表库坝土体结构小，煤层埋藏深度与坝体高度比值较大，常规数值模拟、相似材料模拟和理论计算等方法难以推算坝体的移动变形规律。采用数值模拟方法，按照边长5 m 至少要划分300 多万个单元网格，根据求解时间与网格数大致呈N4/3的正比关系[2-3]，完成一次模型的开挖计算大约需要6.8 d，完成整个模拟需要的时间之长是非常惊人的。利用相似材料按照1∶200 的相似比进行模拟计算，建立至坝体的全尺寸模型过大，现有装置满足不了[4-6]。依据煤层上覆岩层的矿压显现规律推断坝体的移动变形，由于煤层埋藏深度过大，顶板岩层结构复杂，推算结果也难以保证其准确度[7-8]。为此，提出将地表和坝体结构的移动变形分开研究，定量地分析地表移动变形下坝体结构的移动变形，获取坝体结构随工作面不同推进长度下的移动变形规律，建立坝下深部煤层开采坝体结构局部放大的时空序列反演机制，能够为库坝下煤层开采坝体移动变形机理的分析和坝体的安全性评价提供理论指导依据。

1 分析方法的构建

1.1 方法的提出

坝下煤层开采后，采空区覆岩的移动变形自下而上，首先波及到坝体周围的地表，然后带动坝体一起移动变形。将此传递过程分解为2 个部分[9-12]，一是工作面开采引起的覆岩及地表的移动变形，得到随工作面推进形成不同尺寸的下沉盆地；二是逐渐扩大的下沉盆地带动坝体的移动变形，得到随下沉盆地不断扩大后坝体不同部位的移动变形。如果能够舍去工作面开采和下沉盆地形成这一过程，按照下沉盆地形成的过程直接开挖盆地，从而观测分析坝体的移动变形和破坏，将会大大缩小模型尺寸、减少运算时间和提高小目标的分析精度，坝下开采覆岩及地表移动变形过程分解如图1。

1.2 模型构建及实施步骤

由图1 可以看出，如果直接在地表进行开挖，坝体的移动变形会立刻发生相应的移动变形，甚至破坏。这与深部工作面开采引起的缓慢、连续变形不相符的，为了使模拟能够再现煤层开采引起的坝体变形破坏，在地表以下保留一定厚度的缓冲岩土层，根据采空区覆岩变性破坏的“三带”结构，留设的岩土层厚度至少要大于垮落带的高度，才能保证地表不切冒。

图1 坝下开采覆岩及地表移动变形过程分解示意图Fig.1 Decomposition diagrams of mining overburden and surface movement and deformation process under the dam

用于挖掘的盆地选择由数值模拟或地表移动变形预计得出的下沉盆地。开挖步距可分成2 部分：一是在坝体下沉之前设置较大的步距；二是在坝体下沉时设置较小步距。简化模型的开挖顺序如图2。

图2 简化模型的开挖顺序Fig.2 Excavation sequence of simplified model

简化模拟分析的实施步骤如下：

1）依据概率积分法和数值模拟计算出煤层开采后地表的最终下沉盆地，选择最大盆地作为开挖盆地的大小。模型的长宽高依据下沉盆地的影响范围进行确定。

2）按照计算结果中不同推进长度时地表下沉曲线影响范围确定开挖范围，对模型中的下部岩土进行开挖，随时观察地表的移动变形。

3）针对模拟结果中不同推进距离下的地表下沉曲线，逐次进行开挖、观测、记录坝体移动变形规律，直至坝体移动变形稳定。

4）根据各阶段坝体的移动变形状态，总结分析随着工作面推进、地表下沉变化时坝体的移动变形规律。

2 数值模拟

根据局部放大分析方法的内涵，结合地表移动变形预计结果和数值模拟得到的工作面推进过程中地表的下沉盆地，通过浅部开挖下沉盆地，模拟坝体受采动影响下的变形破坏规律，检验分析方法的可行性。

2.1 模拟方案

依据地表移动变形预计影响坝体走向的最大范围为570 m，确定模型的尺寸为600 m×2 300 m×87 m。坝体下覆留设30 m 厚的岩土层，中部留设开挖体，开挖体下覆留设40 m 左右岩层用于消除底板对开挖体的影响。

模型开挖尺寸依据工作面与坝体中心位置之间的距离划分2 种：一种是当间距大于200 m 时，此时开采还没有对坝体产生影响，为节约运算时间，选择工作面每推进100 m 时地表移动变形范围作为开挖步长；另一种是间距在200 m 以内时，为深入观测坝体移动变形，选择工作面每推进5 m 时地表移动变形范围作为开挖步长。

利用概率积分法预计坝下工作面开采后地表任意点的下沉值，将任意点的下沉值对应的三维空间坐标导入到Midas 建模软件中[13-14]，绘制出开采坝下工作面地表移动变形的下沉盆地，按照确定尺寸将下沉盆地从模型中分割出来，在开挖体中分割出工作面方向和推进方向均为曲面的开挖步距，然后按照工作面的推进方向对开挖体依次进行分割，划分网格。岩石力学基础参数见表1。

表1 岩石力学基础参数Table 1 Basic parameters of rock mechanics

在大坝走向方向的坝体中心线上布置114 个观测点，测点间距为5 m，测线长570 m，观测坝体的移动变形情况。

2.2 坝体变形破坏规律及特征

观测不同开挖步距下坝体的位移变化、应力变化和塑性区变化[15-16]，分析坝体的移动变形规律，检验分析方法的可行性。

2.2.1 位移分析

根据工作面不同推进长度下的计算结果，提取坝体测线中点随着工作面开采的监测数据，绘制的坝体随开采变化的下沉曲线如图3。

图3 坝体中心位置的下沉曲线Fig.3 Settlement curves at the center of dam body

由图3 可以看出，工作面推进470 m 时，坝体开始受采动影响，此时坝体中心位置的最大下沉值为0.092 m；推进至1 800 m 时，坝体测线中心位置达到该条件下的最大下沉值7.17 m。

2.2.2 应力分析

由于土体发生变形破坏后其塑性破坏区域无法恢复至原始状态，为了更好地说明坝体的变形破坏过程，选择坝体第1 次达到全部塑性破坏时作为重点分析对象，详细说明坝体的移动变形过程，根据位移变化分析结果，选择工作面推进至500、600、700、800、850、900、910、920、930、2 000 m 的应力变化云图进行分析，坝体中心位置的应力变化云图如图4（图中左侧为迎水坡，右侧为背水坡）。

由图4 可以看出：

图4 坝体中心位置的应力变化云图Fig.4 Diagrams of stress change at the center of dam body

1）在工作面推进500 m 时，坝体迎水坡侧干砌石上端出现约1 m 长的压应力集中区（负为压应力，正为拉应力），最大值为-0.38 MPa，其余区域的应力值基本保持不变。

2）当工作面推进至600 m 时，坝体背水坡侧坝脚处出现应力集中，最大应力值为0.054 MPa。迎水坡侧坝肩位置（干砌石段）应力集中区域沿斜坡向下转移，最大应力值为-0.32 MPa，其它区域的应力值均呈现增大的变化，且为压应力。

3）当工作面推进至700 m 时，坝体背水坡侧坝脚处的应力集中区域依然存在，且应力值由0.054 MPa 增大至0.099 MPa。迎水坡侧（干砌石段）应力集中区域基本保持不变，其上覆土体的应力值为正值。其它区域的应力值继续增大。

4）当工作面推进至800 m 时，背水坡侧坝脚处应力集中区域增大，应力值继续增大（0.378 MPa）。迎水坡（干砌石段）应力集中区域的应力值由-0.316 MPa 增大至-0.693 MPa，应力集中区域也呈现增大的变化趋势，但迎水坡侧坝脚位置（坝基土体）出现了小范围的应力集中现象，最大应力值为0.01 MPa，其它区域的应力值继续增大。

5）当工作面推进至850 m 时，背水坡侧坝脚处应力集中区域减小，其上覆的应力集中区域减小，总体呈现拉应力区域增大，压应力区域减小。迎水坡侧（干砌石段）应力集中区域增大，应力值由-0.693 MPa 增大至-1.13 MPa，坝脚位置（坝基土体）应力集中区域减小，应力值由0.01 MPa 增大至0.80 MPa，坝顶位置出现小范围的应力集中区域，最大应力值为-0.87 MPa。其它区域的应力值逐渐减小，均降低至-0.2 MPa 以下。

6）当工作面推进至900 m 时，背水坡侧坝脚处应力集中区域增大，其上覆的应力集中区域减小，应力值由拉应力转化为压应力。迎水坡侧（干砌石段）应力集中区域继续增大，最大应力值增大至1.7 MPa，坝脚位置的应力集中区域也逐渐增大。其它区域的应力值继续减小，且为压应力。

7）当工作面推进至910 m 时，背水坡侧坝脚处应力集中区域基本保持不变。迎水坡侧（干砌石段）应力集中区域基本不变，处于拉应力状态，坝脚位置的应力集中区域和应力值均在减小，总体呈现为拉应力状态。坝顶位置的应力集中区域增大，应力值由-0.87 MPa 增大至1.1 MPa。坝体中心区域出现了压应力转化为拉应力变化，最大应力值由-0.95 MPa转为0.059 MPa。其它区域的应力值继续减小，且为压应力。

8）当工作面推进至920 m 时，背水坡侧坝脚的应力集中区域消失，应力值降低为-0.021 MPa，上覆的应力集中区域消失，保持在压应力状态。迎水坡侧（干砌石段）应力集中区域减小，应力值增大至1.9 MPa，坝脚位置的应力集中区域减小，应力值也在减小，转变为压应力状态。坝顶位置的应力集中区域和应力值基本保持不变。坝体中部位置拉应力区域扩展至坝体顶部，且随着高度的增大，应力值逐渐减小，几乎为0。其它区域的应力值也逐渐减小，且为压应力。

9）当工作面从920 m 推进至2 000 m 的过程中，背水坡侧、迎水坡侧和坝体中部位置土体的应力值呈现1 个缓慢增大的过程，应力值均为负值（压应力），最大应力值为-0.25 MPa，说明土体在经历压缩变形过程。

由上述坝体各部位的应力变化可以得出：坝体经历了从拉伸变形破坏到压缩-还原的交替变化过程[17-18]，只是各部位出现拉伸-压缩转变过程的时间节点不同，受采动影响开始阶段背水坡侧坝脚先经历拉伸变形破坏，其次为迎水坡侧坝体，最后为坝体中部位置。经历由拉伸转变为压缩的过程也是由背水坡侧坝体开始，其次为迎水坡侧坝体，最后为坝体中部位置。

2.2.3 塑性破坏分析

根据坝体应力变化云图对应的塑性破坏云图如图5，由此分析坝体的变形破坏规律及特征。

由图5 可以看出：

图5 坝体中心位置的塑性破坏云图Fig.5 Diagrams of plastic failure at the center of the dam

2）当工作面推进600 m 时，背水坡侧坝脚出现小范围内的拉伸破坏区域（长1.7 m，高1 m），说明此区域内土体受下部采空区侧拉伸变形较大（应力值由压应力转换为拉应力）。坝体迎水坡侧和坝体中部位置均未出现变形破坏。

3）当工作面推进至700 m 时，坝体背水坡侧变形破坏区域增大，说明此处坝体受到下部表土层的下沉，对其产生的拉应力增大所致，拉伸破坏区域长2.5 m，高1.3 m。坝体迎水坡和坝体中部位置均未出现变形破坏。

4）当工作面推进至800 m 时，坝体背水坡侧产生拉伸破坏的区域继续增大（长4.8 m，高2.9 m）。迎水坡侧坝脚出现小范围的变形破坏区域（长1.6 m，高1 m），说明此处坝体随着背水坡侧坝脚位置的拉伸破坏，坝体中部位置和迎水坡侧坝脚向背水坡侧倾斜变形，迎水坡侧坝脚上覆没有限制其移动变形的实体结构，其受到的拉应力大于土体的自身强度，产生拉伸破坏。坝体中部位置未出现变形破坏区域。

5）当工作面推进至850 m 时，坝体背水坡侧的拉伸破坏区域继续增大（长17.8 m，高4.5 m），此处坝体受到下部表土层的拉伸作用继续增大所致。迎水坡侧坝体变形破坏区域出现急剧增大现象，说明此处坝体受到背水坡侧下沉变形，拉伸幅度加大，迎水坡侧土体受到的拉应力也在增大，拉伸破坏区域增大。

6）当工作面推进至900 m 时，坝体背水坡侧的拉伸破坏区域增大速度减缓，是由于背水坡侧坝体受到下部表土层拉应力减弱所致，背水坡侧坝体由拉应力转化为压应力。迎水坡侧坝脚受到背水坡侧的拉应力减弱，变形破坏区域增大缓慢。而坝体中部位置开始受到下部表土层的拉应力，出现拉伸破坏，此时坝体中部位置（底部）仅有1 m 左右的土体没有被拉伸破坏区域影响。

7）当工作面推进至910 m 时，坝体背水坡侧坝体继续保持压应力状态，坝体迎水坡侧受到的下部表土层拉应力减弱，拉伸破坏区域增长缓慢。坝体中部位置拉伸破坏区域增大，且逐渐向上发展，此时坝体的塑性破坏区域已贯穿坝基，变形破坏区域很有可能成为渗流的主要通道。

8）当工作面推进至920 m 时，坝体迎水坡侧坝脚位置开始由拉伸变形破坏转变为压缩变形破坏。坝体中部位置的拉伸破坏区域继续增大，直至工作面推进至930 m 时，整个坝体处于拉伸破坏区域。坝体中部位置也开始由拉伸变形转向压缩变形状态。

随着工作面开挖结束，坝体一直保持着拉伸破坏状态。虽然土体发生拉伸破坏后，会出现压缩、还原的过程，但软件只能记录出现过变形破坏，塑性破坏区域无法显示还原过程，因此，在工作面继续向前推进的过程中，坝体整体结构均出现了拉伸破坏，仅仅依靠坝体结构的自重使坝体恢复到原始结构状态，需要经历很漫长的历史过程。

青辰想说些什么安慰，却又不知道说什么才好。他能够想象到当时的战况有多惨烈，那日在天葬场，他已经见识了唐飞霄的手段。

3 坝体移动变形观测

利用topconGPS 水准仪和全站仪对工作面开采过程中坝体的移动变形进行观测[19-20]，从坝体中心位置向库区方向在坝顶布设1 条观测线，设置20 个测点，测点间距为30 m，测线全长600 m，坝下开采坝体移动变形测点布置如图6。设计工作面推进至600 m 之前，工作面每推进50 m 观测1 次坝体的移动变形值；工作面推进600 m 至1 800 m 过程中，每推进10 m 观测1 次；工作面推进1 800 ～2 000 m过程中，每推进50 m 观测1 次。

图6 坝下开采坝体移动变形测点布置图Fig.6 Layout of measuring points of movement and deformation of mining dam body under the dam

选取坝体中部1#、2#和3#测点的观测结果，分析坝体受采动影响的移动变形规律，工作面推进过程中3 个测点的下沉值见表2。

表2 工作面推进过程中3 个测点的下沉值Table 2 Subsidence values of 3 measuring points during the advancement of working face

依据表2 中工作面不同推进长度下坝体1#、2#和3#测点的下沉值，绘制的坝体上3 个测点受采动影响下的下沉变化曲线如图7。

由图7 中测点下沉变化曲线可知：1#、2#和3#测点受采动影响的下沉变化趋势基本相同，受采动的影响起始点在600 m 左右，在工作面推进至1 500 m 以后各测点处坝体的下沉基本平稳；工作面推至1 800 m 左右时，其中1#测点和2#测点处坝体达到最大下沉，最大下沉值为7.077 m，3#测点处坝体的最大下沉值为7.032 m。

图7 受采动影响下坝体上3 个测点的下沉变化曲线Fig.7 Subsidence variation curves of three measuring points on the dam body under the influence of mining

实测结果与数值模拟结果对比如图8。

图8 实测结果与数值模拟结果对比图Fig.8 Comparison between measured results and numerical simulation results

从图8 坝体中部位置的下沉情况对比可知：模拟得到坝体中心位置的下沉过程变化曲线与实测结果中1#测点变化曲线基本吻合，模拟得到坝体最终的下沉值为7.17 m，与坝体真实的最大下沉值7.077 m 相比，绝对误差为0.093 m，误差较小。同时，上述坝体的移动变形过程经历了实测工作面采动后地表下沉曲线的各个位置，符合地表移动变形特征，也证明了局部方法分析方法具有一定科学性和合理性。