人造-天然复合岩层冒落机理数值模拟研究*

谭宝会，梁 博，贾凯跃，李明润

(西南科技大学 环境与资源学院，四川 绵阳 621010)

0 引言

下向分层胶结充填法因具有安全、高效、低贫损等优点而被广泛应用于贵重有色金属矿山和稀有金属矿山地下开采中，但当某些特定因素发生改变后(如矿石价格下跌、品位降低等)，矿山可能不再适合采用生产成本较高的胶结充填法进行采矿，转而考虑采用生产成本较低的无底柱分段崩落法。此时崩落法采场顶板以上通常至少包含了两种属性明显不同的材料，其中一种是人工构筑的胶结充填体，另一种是充填体上部的天然围岩，这两种材料相互叠加构成了一种特殊的人造-天然复合岩层。崩落法采场在首采分段采矿过程中需要诱导顶板复合岩层冒落，形成正常采矿所需的散体覆盖层。因此研究复合岩层的冒落机理及模式，对于制订崩落法采矿方案及安全管控措施有着重要意义。但崩落法采空区顶板的冒落通常发生在无法直接观测的崩落区内，这就增加了冒落机理、冒落模式及冒落发展过程的研究难度。

近年来，数值模拟方法的快速发展为顶板围岩的冒落机理研究提供了一条便捷途径。数值模拟方法可以很好地分析地下采矿过程中岩体的位移、应力及应变的分布模式[1-3]。因此有学者将连续介质数值模拟方法[如有限元方法(FEM)、有限差分方法(FDM)]应用于地下采矿引发的顶板冒落机理研究，但是地下采矿诱发顶板岩体冒落为松散岩块，冒落前后岩体力学特性将发生显著改变，因此利用FEM、FDM等方法模拟顶板冒落时，将岩体定义为一种连续的弹塑性材料并利用弹塑性理论将顶板的冒落归因于材料的屈服是有一定缺陷的。而离散元方法(DEM)、非连续变形分析(DDA)以及数值流形方法(NMM)在模拟岩体不连续的大变形、大位移方面具有显著优势，被广泛应用于模拟地下开挖引发的顶板冒落行为[4-7]。DDA方法是一种对不连续块体系统的静态和动态力学行为进行计算的数值模拟方法，其与DEM相比，共同之处在于模拟对象都是离散的块体单元，但由于计算收敛性等因素，DDA方法在模拟不连续块体的大位移、大变形等方面(如岩体冒落、边坡滑坡等)更具优势[8-9]。

综合比较当前可利用的数值模拟工具，采用DDA或者DEM等方法对岩体的冒落机理进行研究更为合适。本文以某矿山下向分层胶结充填法转无底柱分段崩落法为工程背景，采用DDA方法对崩落法采场顶板胶结充填体及上部围岩构成的复合岩层冒落机理、模式及过程展开研究，以期为崩落法采场覆盖层形成方案的设计和实施提供指导。

1 DDA方法的改进

1.1 DDA子块体开裂算法的优化

在以往的预离散子块体DDA开裂模拟方法中，通过采用虚拟节理对子块体进行黏结以模拟连续体的变形，进而根据虚拟节理面上的块间接触力判断沿预设节理面的拉伸或剪切破坏情况；该算法能够较好地模拟岩体的开裂路径和破坏形态，但受预设虚拟节理面方向对块间接触力大小的影响，得到的岩体开裂强度与实际可能存在较大差异。为解决这一问题，倪克松等[10]对DDA方法中的子块体开裂模拟算法进行了优化，在该算法中不再根据块间接触力判断开裂，而改进为根据邻近子块体的应力状态判断开裂，并仍假定裂纹沿虚拟节理面产生。研究表明，改进后的开裂算法具有较高的准确性，并减小了虚拟节理分布对开裂破坏路径模拟结果的影响。

1.2 DDA块体切角功能的实现

一般来说，在DDA模型中划分的子块体是最小的块体单元，无法二次开裂，因此当这些块体带有锐角时在冒落或移动过程中就极易与周边块体相互挤压而卡死，从而影响块体运动行为的发展。为解决这一问题，在DDA前处理程序中添加了块体的切角功能，假设模型中某个块体为锐角三角形，其中一条边长为L，通过设定切角比例k，便可在该条边的首尾两端各切去kL/2，从而消除了块体中的尖角(见图1)。

图1 改进的DDA块体切角功能

实验结果表明，切角后的块体系统发生不正常卡死的概率大大降低。此外，切角后的模型在一定程度上也更符合真实的岩体结构，但要确定合理的切角比例则需要进行大规模的岩体结构信息调查，这在实际操作中难度较大，但从减少块体系统卡死的概率来说，只要进行适当切角，其作用就会非常明显。

2 数值模型构建

2.1 模型的建立及数值控制参数设定

选取垂直矿体走向的典型剖面建立DDA模型，模型从下至上按照实际情况划分为4种不同属性的材料(不包含模型边框)，其中最底层为10 m厚的矿体，该层代表实际中崩落法采场的首采分段。考虑到实际中矿体的强度相对较大，对顶板冒落所产生的影响较小，因此为了尽可能地缩减模型中块体的数量，将矿体材料的网格简化为9个矩形刚性块体，模拟时对位于中部的7个块体从左至右依次开挖，以代表实际中沿进路方向的退采，每个开挖块体的宽度为10 m，大约等同于实际中的4～5个回采步距。

矿体之上为30 m厚的胶结充填体，共包含6个充填分层。胶结充填体中包含了两种不同强度的不连续面：一种为上下分层交界处的水平不连续面，其代表的是实际中胶结充填体上下分层的接触面，该不连续面的强度要小于胶结充填体的强度；另一种不连续面为胶结充填体内部的虚拟节理，其强度可取胶结充填体试样的实测强度。

本次建模采用了蜂巢状网格作为胶结充填体的虚拟节理，以便更真实地模拟胶结充填体的破坏形态，同时考虑到当前版本的DDA程序尚无法实现子块体的二次开裂，这在一定程度上限制了胶结充填体剪切滑移破坏的发展，为了弥补这一缺点，在胶结充填体中增加了一组垂直方向上的贯通虚拟节理。

胶结充填体之上为50 m厚的片麻岩，片麻岩中包含了两组优势结构面，其倾角及倾向分别为88°∠56°和3°∠52°；由于这两组节理均具有较好的连续性，因此在模型中将这两组节理均设定为贯通节理，同时为了缩减模型中的块体数量，在保持这两组节理走向、倾角及接触面力学参数贴近真实的基础上，将节理间距扩大10倍，这种做法虽然会在某些方面降低模拟结果的真实程度(如岩体开裂后的碎胀系数等)，但并不会从本质上影响岩层的失稳模式。

片麻岩之上为10 m厚的第四系表土，由于第四系表土强度极低，十分松散，极易发生破坏，因此为进一步缩减模型中的块体数量，仅在模型中设定了10 m厚的第四系表土为代表，同时采用2组互相正交的虚拟节理来对其强度进行赋值并表征破坏。最终模型中块体数量维持在4.0×103个左右，模型几何尺寸及测点分布情况如图2所示。

模型中各材料的关键物理力学参数是根据矿山提供的地质报告以及一系列相关的物理实验获得的。为了使模型中各材料的破坏更贴近真实情况，设定当模型中某个不连续面发生破坏后，该不连续面上的黏聚力、抗拉强度、摩擦角等强度参数均会发生衰减[11]。模型中各材料的物理力学参数见表1。

图2 DDA模型几何尺寸及测点布置

表1 模型中各岩层的关键物理力学参数

DDA方法所用到的数值控制参数是由用户自定义的，合理的数值控制参数通常要经过大量的前期试算才能得到。采用合理的数值控制参数在进行求解计算时，每一时步的计算会在最小次数的开-合迭代之后达到收敛[12]。本次模拟使用的数值控制参数见表2。

表2 模型求解所用的数值控制参数

2.2 模拟方案设定及初始应力平衡求解

为了使模拟更贴近现实，模型在开挖之前需要先在静态条件下求解初始应力平衡，即在模型中先实现地应力的还原再进行开挖。经现场地应力测试，该矿区地应力以水平应力为主，垂向压力基本为上覆岩体自重，各测点的水平应力与垂向应力之比(侧压比)λ大约维持在1～2，且最大水平应力方向基本垂直于矿体走向。为研究水平应力对复合岩层冒落发展的影响，本次数值模拟构建了两组模拟方案，即方案1和方案2，方案1侧压比λ=1，方案2侧压比λ=2。

图3以模型中1#测点为例展示了两组模拟方案的初始应力平衡求解过程。1#测点位置对应的实际埋深在130 m左右，由理论计算可知，该测点处的垂向压力σv=3 MPa。因方案1侧压比为1，故该测点的水平应力σh=3 MPa；而方案2侧压比为2，故该测点的水平应力σh=6 MPa。

图3 模型求解初始应力平衡

从图3可以看出，两组模型在物理时间t=20 s时(共计算20 000时步)达到了初始应力平衡。在方案1中，1#测点处的垂向应力达到3 MPa，水平应力也达到3 MPa。在方案2中，1#测点的垂向应力达到3 MPa，水平应力基本达到6 MPa。可见两组方案的应力条件基本达到预期状态，此时便可以从左至右依次移除底部所设定的7个块体，执行开挖后程序自动转为动态条件进行计算。模型中第1步开挖设定在20 000时步，第二步开挖设定在60 000时步(t=24 s)，之后每隔10万时步进行一次开挖，两次开挖之间保持较长时间间隔的目的是使每次开挖引发的顶板冒落都能够得到充分发展。

3 模拟方案及结果分析

图4给出了模拟方案1中复合岩层随着开挖逐渐发生失稳冒落的动态演化过程。由图4(a)可知，在第1步开挖后，形成了跨度为10 m的采空区，此时空区顶板围岩几乎未发生任何开裂和破坏。由图4(b)可知，在第2步开挖后，空区跨度达到了20 m，此时顶板胶结充填体开始发生破坏，胶结充填体沿空区两侧逐层向上破坏，很快形成贯通破坏裂隙，延伸至第五分层时才中止，此时上部片麻岩也发生了部分开裂，最终第一层胶结充填体岩梁发生了弯曲离层，但未冒落，而是在岩梁两端摩擦阻力及水平应力作用下形成了稳定的裂隙梁，在此过程中有个别散块冒落至空区底板。由图4(c)可知，在第3步开挖后，空区跨度达到了30 m，此时空区上方的胶结充填体、片麻岩和第四系表土瞬间完全垮塌冒落，以柱状形式整体滑落至空区底板，此阶段可视为空区顶板突发大规模冒落。由图4(d)可知，第4步开挖后，新开挖空区上部的胶结充填体及覆岩呈现为临空的悬臂梁形式，悬臂梁在自重作用下于支点处发生折断并冒落至空区。为节省计算时间，在第4步开挖引发的冒落完成后终止了计算，可以推断的是在后续退采中新开挖空区上方的复合岩层将周期性地发生悬臂梁折断冒落。

(a)开挖跨度10 m

图5为模拟方案2中复合岩层随着开挖逐渐发生失稳冒落的动态演化过程。

(a)开挖跨度10 m

由图5(a)可知，在第1步开挖后，空区跨度为10 m，此时空区顶板各岩层几乎未发生任何开裂和破坏。由图5(b)可知，在第2步开挖后，空区跨度达到了20 m，此时仅有空区顶板胶结充填体的第一分层产生了少许离层裂纹。由图5(c)可知，在第3步开挖后，空区跨度达到了30 m，此时第一分层胶结充填体冒落至空区底板，同时上部分层中也开始产生少量裂纹。由图5(d)可知，在第4步开挖后，空区跨度达到了40 m，此时第二分层胶结充填体发生了弯曲离层，第3-第5分层两端发生贯通式开裂，但并未滑落，而是形成了稳定的裂隙梁，同时空区上方对应的片麻岩和第四系表土中产生大量裂纹。由图5(e)可知，在第5步开挖后，空区跨度达到了50 m，此时胶结充填体和顶板围岩突发大规模整体滑落。由图5(f)可知，第6步开挖后，新开挖空区上部复合岩层以临空悬臂梁折断的形式发生破坏。在模型第6步开挖引发的冒落活动稳定后，终止了计算，可以推断的是在后续退采中新开挖空区顶板围岩将周期性地发生悬臂梁折断垮塌。

对比方案1和方案2的冒落过程可知，在该矿山胶结充填体和顶板围岩的强度和结构特征条件下，崩落法采场空区顶板复合岩层顶板很容易发生柱状整体陷落，从而引发直通地表的大规模冒落，而且水平应力显著提升了各岩层的抗破坏能力，增加了顶板的临界冒落跨度，但采空区的跨度越大，突发冒落的规模也就越大，更容易引起冒落冲击灾害。

为深入分析空区顶板的冒落机理，图6以方案1为例给出了模型计算过程中的最大剪应力分布云图，以及空区上方岩层整体冒落后模型中的破坏裂纹分布模式。

(a)第60 600时步最大剪应力云图

由图6可知，模型开挖后空区两侧发生了明显的剪应力集中，在初期时剪应力主要集中在胶结充填体第一分层岩梁两端[见图6(a)]，率先致使第一分层岩梁两端发生剪切破坏；随后，集中剪应力转移至第二分层，使第二分层岩梁两端发生剪切破坏，以此类推，岩梁两端的剪切破坏逐层向上发展[见图6(b)、图6(c)]。由于胶结充填体强度较低，剪切裂纹很快就贯通了整个胶结充填体，同时胶结充填体在自重作用下发生弯曲下沉，致使岩梁中部产生拉应力集中，从而发生拉伸破坏[见图6(d)]，严重破坏了胶结充填体的完整性。当岩梁两端所受的水平挤压力和摩擦力不足以阻止胶结充填体滑落时，胶结充填体便会发生冒落。失去胶结充填体的支撑后，上部片麻岩在重力作用下沿陡立的贯通节理发生整体剪切下滑，同时强度较低的片麻岩也随之陷落。总之，由于胶结充填体和顶板围岩强度均较低，空区顶板围岩往往会在较短的时间内发生整体垮塌，从而造成大规模冒落。

为了展现开挖引发的应力重分布规律，图7以模拟方案1为例给出了模型中1#-4#测点的垂向应力演化曲线。由图7可知，在模型连续开挖到40 m跨度的过程中，当测点位于空区一侧时，测点处将产生垂向应力集中现象，且距离空区越近，垂向应力集中现象越明显；而随着开挖的进行，一旦测点的位置由空区一侧转为空区顶板之上后，其垂向应力被立即释放，直至测点所在块体冒落至空区底板，上部冒落岩体逐渐堆积在该块体上，该测点的垂向压力才逐渐出现回升，此时该测点上的垂向应力即为上部散体产生的垂压。在整个开挖过程中，3#测点和4#测点始终位于空区一侧。由图7(c)可以看出，3#测点在前两步开挖后垂向应力基本未发生改变，当第三步开挖后，垂向应力由3.0 MPa增加到3.5 MPa左右，增加幅度为17%左右，可认为此时该测点受到了采动压力的影响；由图7(d)可以看出，4#测点的垂向应力自始至终基本不变，表明在开挖过程中该测点未受到采动压力的影响。

(a)1#测点

数值模拟结果表明，采空区顶板发生冒落的根本原因是开挖造成模型中应力重新分布；具体来说是开挖后空区上方垂向应力向空区两侧转移，致使空区两侧产生剪应力集中，胶结充填体在空区两端的集中剪应力作用下逐层向上发生剪切破坏，同时横梁弯曲下沉时在中部产生拉应力集中，引发拉伸破坏，顶板围岩在失去胶结充填体的支撑后，便在自重作用下沿软弱的陡立结构面发生柱状剪切滑移破坏。

由于该矿区矿体埋藏较浅，水平应力较大，同时胶结充填体强度较低而顶板围岩富含弱陡立节理，在这些因素的共同作用下，崩落法采场回采过程中顶板复合岩层易突发直通地表的大规模冒落。此时，气浪冲击危害和冲击振动危害将是最可能致灾的两大危险源。工程实践证明[13-15]，在崩落法首采分段退采过程中均匀扩展采空区，避免在空区内部留支撑矿柱，可促使空区顶板及时发生冒落；同时在空区底板预留足够厚度的矿石散体垫层，既可以有效缓冲顶板冒落的冲击振动作用，还可封堵出矿口以预防气浪冲击。

4 结论

a.数值模拟结果表明，随着采空区的不断扩大，空区上方垂向应力逐渐向空区两侧转移，致使空区两侧产生剪应力集中，胶结充填体在空区两端的集中剪应力作用下逐层向上发生剪切破坏；由于胶结充填体的强度较低，剪切破坏很快贯通整个胶结充填体，同时在自重作用下岩梁发生弯曲下沉，致使岩梁中部产生拉应力集中，从而发生拉伸破坏，严重破坏了胶结充填体的完整性；当岩梁两端所受的水平挤压力和摩擦力不足以阻止胶结充填体滑落时，胶结充填体便发生冒落。

b.由于矿山片麻岩中富含贯通的光滑陡立节理，当胶结充填体发生失稳下滑后，上部片麻岩便失去支撑，极易在重力作用下沿陡立节理发生整体剪切滑移破坏，同时由于第四系表土的抗剪强度极低，当其临空后也极易发生整体陷落。

c.研究结果表明，由于该矿区矿体埋藏较浅，水平应力较大，同时胶结充填体强度较低而顶板围岩又富含连续性较好的光滑陡立节理，因此在崩落法回采过程中采空区顶板复合岩层易突发直通地表的柱状大规模冒落，此时气浪冲击和冲击振动将是复合岩层冒落过程中最可能致灾的两大危险源，在实际生产过程中，建议采取均匀扩展采空区和预留矿石散体缓冲垫层等措施进行防控。