随着我国能源需求量和开采强度的逐年增加，煤矿正逐步向深部开采状态发展[1-2]。由于深部煤矿巷道围岩处于高地应力、高地温、高岩溶水压的“三高”地质环境，且岩体本身力学性质普遍较为软弱，在开挖扰动与卸荷作用下极易诱发围岩破裂碎胀和挤压变形，从而严重威胁煤巷的稳定与安全[3]。
深部煤巷围岩破裂碎胀挤压大变形的力学本质为从峰前损伤扩容连续变形(细观裂纹萌生、扩展、汇聚)到峰后破裂碎胀非连续变形(宏观裂隙贯通以及块体错动、翻转和分离)的渐进破坏演化过程[3-4]。就目前而言，无论是理论解析方法、室内试验还是原位试验，尚无法对深部围岩这一复杂变形破坏行为及规律进行定量描述与精准预测[5]。
随着计算机科学技术的高速发展，数值仿真逐渐成为岩土工程领域的主要辅助研究手段。就目前而言，关于岩石破裂过程模拟分析的常用数值方法主要可分为3类:基于连续的数值方法，基于非连续的数值方法，以及连续-非连续耦合的数值方法。基于连续的数值方法，常用的如有限单元法(Finite Element Method,FEM)[6]、边界单元法(Boundary Element Method,BEM)[7]和有限差分法(Finite Difference Method，FDEM)[8]等，主要通过构建特定本构模型以及嵌入相应力学准则进行岩石体力学性质和裂纹萌生扩展演化过程的模拟仿真，然而上述方法难以反映煤巷围岩的宏细观非连续性，更无法有效捕捉大量复杂裂纹扩展贯通所导致的破碎岩石块体的挤压、摩擦和移动、转动等大位移以及破裂碎胀大变形力学行为特征。基于非连续的数值方法，常用的如离散单元法(Discrete Element Method,DEM)[9]、非连续变形分析方法(Discontinuous Deformation Analysis,DDA)[10]和刚体弹簧元法(Rigid Body Spring Method,RBSM)[11]等，在模拟岩石体峰后破裂碎胀和挤压大变形方面具有较大优势，然而上述方法预先对求解域进行离散处理，难以真实反映岩石体峰前阶段损伤扩容连续变形。此外，在块体自身应力应变场精度计算、块体自身开裂模拟等方面也存在较多缺陷[12]。
由石根华提出的数值流形方法(Numerical Manifold Method，以下简称NMM)[13]是一种典型的连续-非连续耦合数值方法，其采用了两套覆盖系统——数学覆盖和物理覆盖，能够在同一计算架构下灵活而有效地处理连续与非连续问题;该方法所采用的接触算法，能够很好地模拟非连续面之间的挤压、摩擦、张开等相互作用以及块体的平移、转动等大变形特征[14]。笔者在原始NMM程序基础上进行改进，发展了一种深部煤巷围岩破裂碎胀大变形模拟预测分析方法，主要研究工作如下:
(1)在NMM前处理程序中植入Voronoi随机多边形自动生成算法，生成能够较好表征煤巷围岩体非连续结构面的Voronoi随机多边形块体环路集合体;基于该几何表征发展了能够考虑岩石体非连续面黏聚力和张拉强度的界面接触计算模型，构建了基于界面接触模型的破裂算法准则。
(2)发展了基于显式积分格式的数值流形方法，通过Voronoi随机多边形环路自动生成算法和界面接触计算模型的嵌入，搭建了基于Voronoi多边形环路构型的显式数值流形方法。
(3)在改进的NMM方法平台基础上，发展了能够表征巷道掌子面纵向支撑效应的开挖模拟算法，从而完成深部煤巷隧道围岩破裂碎胀大变形模拟预测分析方法的构建。
(4)基于上述模拟预测分析方法，对深部煤巷开挖卸荷作用下围岩渐进破坏演化规律与碎胀大变形力学行为特征，以及最终开挖损伤区的形成与分布特征等进行了详细研究。
1 显式数值流形方法
1.1 有限覆盖系统与位移近似函数
NMM的覆盖系统由数学覆盖(Mathematical Cover，MC)和物理覆盖(Physical Cover，PC)构成，物理覆盖是物理网格对数学覆盖的再剖分，物理覆盖的重叠部分定义为流形单元(Manifold Element)[13-14]。在原始NMM中一般采用常规三角形网格作为数学网格，形成的数学覆盖为正六边形数学覆盖(图1中M1,M2和M3)，通过与物理网格(图1中红色实线表示)的切割形成物理覆盖，其中每3个物理覆盖相互重叠构成一个流形单元(图2)。
图1 三角形数学网格与物理网格
Fig.1 Triangular mathematical mesh and physical mesh
图2 物理覆盖与流形单元
Fig.2 Physical covers and manifold elements
NMM在每个物理覆盖上构建局部位移近似函数，通常采用多项式进行局部位移函数
近似:
(1)
其中，x为位置矢量矩阵;PT为多项式的基函数矩阵;d为位移自由度列阵;i为对应物理覆盖的编号。通过权函数对局部位移近似函数进行加权平均，可以得到整个覆盖系统上的总体位移函数uh为
(2)
其中，npc为构成流形元的物理覆盖数量，在原始NMM中值为3。将式(2)代入式(1)，则全局位移函数可表示为
(3)
其中，Ni(x)为形函数;wi(x)为覆盖i对应的权函数，且满足以下条件:
(4)
1.2 接触算法与开闭迭代
接触算法是非连续数值分析方法的主要功能模块，用于模拟块体之间的挤压、摩擦、分离等相互作用关系，继而计算块体的变形和大位移。原始NMM继承了DDA的接触检索和接触计算方法，基于块体间无嵌入、无张拉的原则，通过开闭迭代进行接触计算的收敛判断。
在NMM中以块体环路为基本个体，采用矩形框接触检索方法对所有环路进行接触搜索。NMM中首先给定一个最大位移阈值d0，即在每个时间步长内，所有流形单元节点位移的最大值不超过该位移阈值。在接触检索过程中，对每一个块体环路定义矩形包络圈(图3)，使该环路恰好被该矩形框完全包络。在每一时步的开始，通过两两块体环路上环路节点的循环遍历，计算其所在矩形包络圈的最小距离，此时如果该距离小于预先给定的接触判断阈值(一般取2d0)，则该两个块体环路被判定为潜在接触对，并对相应环路编号和环路节点编号进行存储以便后续的接触计算。
图3 潜在接触环路对与非潜在接触环路对
Fig.3 Potential and non-potential contact loop pairs
在接触检索计算后，NMM对每一组潜在接触环路进行接触类型判断，从而找出存在的接触对。在NMM中可能存在的接触对类型可以分为以下3类:角-角接触，角-边接触，以及边-边接触(图4)。其中边-边接触可以分解为两组角-边接触对。NMM的接触计算是基于时步进行的，上一时步末的所有几何物理参数将传递至下一时步初。为了使接触计算结果更为准确，NMM在每一时步内反复求解控制方程并进行接触状态的判断，其中伴随着罚弹簧的加设和移开过程即为开闭迭代[14]。
图4 二维NMM中的3种接触类型
Fig.4 Three contact types in 2-dimension NMM
1.3 NMM显式求解算法与控制方程
原始NMM是一种基于隐式积分格式的数值方法，由于每一时步都需要重新装配刚度矩阵并对系统方程组进行超松弛迭代求解，造成计算效率偏低，特别是涉及大规模非连续问题求解时，该问题尤为突出[15]。本文根据DEM中采用的显式差分格式与牛顿第二运动定律，构建了基于显式时间积分格式的数值流形方法(ENMM)[16]。
在原始NMM中，系统的运动方程可由最小势能原理推导得到:
(5)
其中，{Un}为第n时步PC的位移增量列阵;{Fn}为作用在PC上的荷载列阵;[M],[K]和[C]分别为质量矩阵、刚度矩阵和阻尼项[15]。假定每个时步内加速度为常量，基于直接积分方法求解第n+1时步的速度项
和位移项{Un+1}:
(6)
其中，Δt为时间步长。基于位移和速度的初始条件，并将式(6)代入式(5)得到NMM系统控制方程:
(7)
式中
(8)
本文采用类似于DEM中的显式积分方法[17]，基于中心差分格式计算每半个时步系统中各个PC的速度和位移增量:
(9)
将式(9)代入式(5)即可得到ENMM的控制方程为
(10)
其中，[M]为质量矩阵，可采用集中质量矩阵的方法，通过每行求和技术进行对角化而得
为荷载列阵:
(11)
式中，
为内力项;
为阻尼项。
2 Voronoi块体环路与界面接触模型
在宏观层面上，深部煤巷围岩体是由岩石基质和各级非连续结构面(如断层、层理、节理等)构成[18];在细观层面上，岩石由大量矿物颗粒、颗粒间胶接界面以及细观非连续(如微缺陷、微孔洞、夹杂等)构成[19]。在进行深部煤巷围岩力学行为模拟计算和评价分析时，关键问题在于构建能够合理表征围岩体宏细观随机非连续面分布特征与空间镶嵌状态的几何结构模型。采用Voronoi随机多边形进行岩石体随机非连续面的几何近似，是岩体力学性质与变形破坏行为研究的新方法。该方法的主要优势在于能够较为合理地表征岩石体宏细观非连续结构面的空间几何分布与镶嵌结构特征，且在很大程度上弱化了岩石体开裂路径对裂隙网格分布的依赖性[20]。基于上述优势，Voronoi随机多边形表征方法已被广泛运用于岩石体的力学性质模拟分析，并取得了较好的效果[18-22]。
本文通过Voronoi随机多边形自动生成算法在NMM前处理程序中的嵌入，生成满足要求的Voronoi随机多边形环路，并在该环路构型基础上发展了能够考虑随机非连续面黏聚力和张拉强度的界面接触算法模型。
2.1 Voronoi多边形环路生成算法
本文在NMM前处理程序基础上开发了一种Voronoi随机多边形自动生成算法，具体计算流程如图5所示。首先，在某一给定计算域生成随机点(图5(a));基于生成随机点的逆时针有序连接产生Delauny三角形网格(图5(b))，并计算网格中每个三角形的重心(图5(c));最后再次通过逆时针有序连接三角形重心，获得计算域的Voronoi多边形集合体(图5(d))。在该算法中，随机点的分布能够通过计算迭代次数进行控制，且计算迭代次数越多，随机点的分布(即各随机点之间的距离)越均匀，从而最终生成的Voronoi多边形集合体几何形状越均匀;计算域中随机点的总数量由人为预先确定，从而决定所生成Voronoi随机多边形的平均尺寸。通过人为控制计算域内随机点数目和迭代次数，从而能够生成满足计算要求的Voronoi随机多边形集合体。
图5 Voronoi多边形块体自动生成算法示意
Fig.5 Generation procedures of Voronoi polygons
图6 Voronoi多边形接触环路节点更新示意
Fig.6 Updating of points on Voronoi polygon loops
基于Voronoi随机多边形自动生成算法计算得到的Voronoi多边形集合体将作为NMM前处理中的物理网格，通过对常规三角形数学网格的再剖分，产生NMM主程序计算所需的接触环路、流形单元以及相应物理覆盖等基本数据信息。
2.2 Voronoi环路节点更新算法
在原始NMM程序的接触检索和计算过程中，对潜在接触环路对上的所有环路节点(包括环路角点和环路边节点)进行遍历判断(图6(a))。然而在引进Voronoi随机多边形集合体近似岩石细观结构或岩体非连续结构面后，计算模型中将不可避免地出现大量接触环路，此时如果仍对每个环路上所有节点进行循环遍历接触检索，将直接导致接触计算模块耗时过长，从而导致程序整体计算效率偏低。
为了提高程序计算效率，本文采用了一种基于Voronoi环路构型的节点更新算法，仅保留多边形块体角点并删除环路边界其余节点，从而大幅削减参与接触检索的环路节点数目。以模型中任意环路i(图6(a))为例，基于所有环路节点生成首尾逆时针相连的有序向量，通过计算首尾相接两向量的夹角余弦值，从而判别中间点是否为多边形环路的角点。如图6(a)所示，环路i中任意3个相邻节点k，k+1和k+2构成了两个首尾相接的向量nk,k+1和nk+1,k+2，当两向量满足
(12)
此时可判断两向量为非共线向量，则节点k+1为环路角点。通过该方法对所有环路进行角点检索，同时更新并压缩存储角点编号、坐标为新环路节点编号、坐标，最后得到一套只含环路角点的节点信息(图6(b))。
2.3 基于Voronoi环路的界面接触模型
原始NMM采用了基于罚函数的接触计算方法，通过侵入点(或边)的法向嵌入量和切向滑移量计算各接触对之间的接触力。然而在原始接触计算方法中，忽略了非连续面之间的黏聚力和抗拉强度，导致模拟计算结果与实际情况存在一定偏差;此外，模型中不可避免地存在大量接触对，接触检索与计算过程将造成巨大的计算量，导致计算效率偏低。本节基于环路节点更新后的Voronoi多边形构型，提出了一种改进的界面接触计算模型，能够很好地反映非连续面间的黏聚力和抗拉强度;此外，由于在接触胶结界面破坏前对应Voronoi多边形环路的边不参与接触检索与计算，从而能够有效提高模型的计算效率。
由于所生成的Voronoi随机多边形块体环路在受荷发生变形破坏之前是紧密贴合且无间隙的，因此假定环路间的相邻边存在一个无厚度的晶间胶结界面(图6(b))，在该胶结界面的构型基础上搭建考虑界面黏聚力和抗拉强度的界面接触模型。基于胶结界面可将岩石的破坏分为纯张拉破坏(I型)、压剪破坏(II型)和拉剪混合破坏(I-II型)3种主要模式[23](图7)。
图7 基于Voronoi块体环路构型的胶结界面接触模型 本构关系
Fig.7 Constitutive law of the interface contact model imbedded into the Voronoi polygon loops
2.3.1 I型破坏模式
当接触界面受拉时，界面受到法向黏结应力σ并在发生张开初期遵循线弹性本构行为(图7(a));一旦其开度(O)达到某一临界值Op界面随即发生屈服，此时该临界值称为接触界面的屈服开度:
Op=2Lft/Pfn
(13)
其中，L为接触胶结界面的长度;ft为岩石材料抗拉强度;Pfn为人为定义的接触界面法向刚度。在接触界面发生屈服后呈线性软化状态，当接触面开度达到残余开度Or时即发生I型张拉破坏。在线性软化阶段(Op<O≤Or)，接触界面法向黏结应力的降低由张拉损伤变量Dt决定:
Dt=(O-Op)/(Or-Op)
(14)
其中残余开度Or可由材料I型断裂能和抗拉强度共同决定:
(15)
I型破坏模式下接触界面的本构关系可表示为
(16)
2.3.2 II型破坏模式
当接触界面处于压剪状态时，界面受到一个切向的应力分量τ，该切向应力分量由接触界面受到的法向应力σ以及切向滑移量S决定(图7(b))。与I型张拉破坏模式类似，在切向滑移初期界面遵循线弹性行为直到切向滑移量达到临界值Sp，该临界滑移量由岩石材料的抗剪强度fs决定:
fs=c+σtan φ
(17)
其中，c和φ分别为岩石材料的黏聚力和内摩擦角。类似于I型破坏模式，在峰后阶段，接触界面遵循线性软化本构行为，且剪应力由剪切损伤变量(Ds)决定:
Ds=(S-Sr)/(Sp-Sr)
(18)
随着滑移量S的进一步增大，切向应力逐渐降低至残余摩擦力fr:
fr=σtan φ
(19)
II型破坏模式下接触界面的本构关系可表示为
(20)
临界滑移量Sp和残余滑移量Sr可根据下式计算得到:
(21)
式中，Pfs为人为定义的界面切向刚度;
为岩石材料的II型断裂能。
2.3.3 I-II型破坏模式
对于张-剪混合模式，接触界面的破坏由界面间的开度和滑移量共同决定(图7(c)):
(22)
晶间胶结界面在发生破坏前并不参与接触检索，而待其失效破坏后对相应块体环路进行标记，该胶结界面对应的两个环路边参与下一时步的接触检索和接触计算，这一改进能够在很大程度上提高程序的计算效率。通过基于胶结界面接触模型的破坏，实现对岩体受荷过程中裂纹的萌生、扩展和贯通全过程进行模拟;此外，Voronoi随机多边形环路集合体构型为捕捉胶结界面破裂后产生的非连续面间的挤压摩擦等相互接触，以及块体的运动翻转等大变形现象提供了可能。
3 开挖卸荷模拟算法
本节在前述章节算法平台的基础上，发展了一种基于软化单元法和删除单元法的煤巷开挖模拟方法。在初始地应力平衡之后，通过模型开挖核心筒内单元的软化和删除，最终完成开挖过程的模拟仿真。
3.1 软化单元法
在NMM方法的前处理过程中，能够较为容易地完成开挖核心筒内流形单元及Voronoi环路和非开挖区流形单元及Voronoi环路对应材料号的区别赋值(图8)，为软化单元的实现提供了较大便利。假设核心筒内单元材料的原始弹性模量为E0，通过预先设定的软化系数α对其进行依时步折减以计算得到当前单元材料弹性模量Ec:
Ec=αE0
(23)
图8 模型区域单元材料赋值及开挖步序
Fig.8 Material parameter assignment and excavation sequence of the numerical model
将折减后的当前弹性模量Ec代入下一时步初，实现开挖核心筒内单元弹模的逐渐软化，进而对煤巷开挖作用下掘进工作面的纵向支撑效应进行考虑。这一简化方法目前已有部分学者在隧洞准静态开挖二维模拟中进行尝试，并取得了较好的模拟效果[20,23]。
3.2 删除单元法
在上述软化单元步序完成之后，需要对隧道开挖面以内单元及环路进行挖除，从而最终完成开挖的模拟。如图8所示，在NMM前处理工作中，可以根据隧洞直径及中心点位置确定开挖边界，通过线段求交算法确定被开挖环路和单元并进行标记，最后完成开挖的计算模拟。具体原理及方法步骤如下:
(1)确定开挖边界。根据计算模型的开挖几何参数(隧洞开挖半径及中心位置)在前处理过程中预先确定开挖边界，即由多线段首尾相接形成开挖环路，并初步获得整个开挖计算模型的单元、节点及覆盖系统等相应几何数据文件。
(2)确定被挖除单元和环路。通过线段求交算法将开挖部分(开挖边界内)与非开挖部分(开挖边界外)定义为不同的Group，并对处于不同Group的环路、流形单元和物理覆盖分别进行标记。
(3)对开挖与非开挖环路单元赋予不同材料属性。在输入模型物理力学参数时，对开挖与非开挖Voronoi块体环路及单元的材料属性进行分组，在此基础上能够较为方便地实现核心筒内流形单元的软化，进行巷道开挖卸荷作用的模拟。
(4)删除开挖环路、单元、覆盖以及信息更新。在核心筒内单元软化到预定值后，进行环路和单元的挖除计算。在被开挖环路和单元完成标记后，通过对原始前处理生成文件进行单元节点编号、物理覆盖编号、环路节点及编号等信息的压缩更新来删除被标记的环路和单元，从而实现挖除模拟。
基于软化单元法和删除单元法的深部煤巷开挖模拟算法流程如图9所示。
图9 深部煤巷NMM开挖算法具体流程
Fig.9 Flow chart of NMM excavation algorithm in deep roadway
4 算例分析
4.1 模型参数
通过采用第3章节中的开挖模拟方法，对某一深部高应力马蹄巷软弱围岩在开挖卸荷作用下的破裂碎胀演化规律以及收敛变形特征进行模拟研究:煤巷围岩宏观力学性质见表1。该计算巷道为马蹄形断面，巷道宽为5.8 m，边墙高度为2.0 m;模型采用0.3 m的Voronoi多边形块体环路进行围岩体随机结构面分布的近似，几何边界取为60 m×60 m(约为10倍洞径)，从而在保证计算效率的同时尽量消除边界效应的影响;模型施加地应力为水平应力25.8 MPa、竖向应力21.9 MPa;该模型主要微观力学参数由单轴压缩和巴西劈裂数值试验标定得到(表2)。
4.2 模拟结果分析
图10给出了该马蹄巷模型在不同开挖时步软弱围岩碎胀大变形裂纹扩展演化规律。由图10可知，巷道拱肩与底角部位由于开挖卸荷作用产生应力集中，导致首先发生剪切裂纹的萌生与扩展(图10(a))，且随着开挖截面单元的逐步软化，裂纹逐渐向远场扩展贯通(图10(b),(c));由于在该应力分布条件下，巷道拱顶与底板应力集中系数较高，导致裂隙发育程度较高从而最终造成围岩高度破碎(图10(d))。
表1 完整岩石及岩体宏观力学性质
Table 1 Macro-mechanical parameters of the intact rock and jointed rock mass
完整岩石力学性质σc/MPaEr/GPaRQD岩体力学性质σcm/MPaEm/GPa抗拉强度/MPa20.2～28.810.4～17.949～7514.532.491.1524.5(平均)14.2(平均)62(平均)
表2 标定所得模型微观力学参数
Table 2 Calibrated micro-mechanical parameters of the NMM model for coal roadway excavation simulation
Voronoi块体微观参数密度ρ/(kg·m-3)弹性模量Em/GPa泊松比νm胶结界面微观参数法向刚度Pfn/GPa切向刚度Pfs/GPa摩擦角φj/(°)黏聚力cj/MPa抗拉强度ftj/MPa2 3702.50.310025233.51.2
图11给出了该模型在不同开挖时步围岩收敛变形矢量分布演化规律。计算结果显示，巷道拱顶与拱肩由于应力集中程度较高导致岩体被裂隙切割贯通并发生极为严重的碎胀变形，所产生的收敛量较大，其中拱肩最大收敛变形达到0.43 m，右侧拱肩最大收敛变形达到0.47 m;底板围岩碎胀程度较拱顶与拱肩较低，最大隆起部位发生于底板左侧，最大隆起量为0.29 m;由于巷道两帮裂纹发育相对较少、两侧边墙处围岩完整性相对较高，因此对应收敛位移也相对较小，分别为左侧边墙0.21 m、右侧边墙0.17 m。经计算，巷道左右边墙及底板的收敛应变分别为7.2%，5.9%和10.0%，处于5%～10%，表明围岩发生了较为严重的挤压变形[24];将该模型预测数据与淮南朱集煤矿某-885 m水平高地应力巷道(最大主应力可达31.2 MPa)围岩变形监测数据进行对比，其演化规律具有较高的一致性[25]。总体而言，上述模拟结果较好地揭示了深部高应力煤巷软弱围岩破裂碎胀挤压大变形的力学行为特征。
图10 深部高应力煤巷围岩碎胀大变形演化规律
Fig.10 Fracturing expansion evolution process of surrounding rock mass in deep coal roadway
图11 深部高应力煤巷围岩收敛变形演化规律
Fig.11 Convergence deformation of the surrounding rock mass in the deep coal roadway
图12给出了该模型在开挖核心筒单元软化系数α=0时(即表征开挖掘进工作面远离该截面导致纵向支撑效应完全消失)对应的开挖损伤区分布特征。由计算结果可得，在马蹄形巷道四周形成了一个11.9 m×10.4 m的开挖损伤区，且该损伤区的宽度和高度均超过了巷道设计宽度和拱顶高度的2倍;巷道两帮损伤区深度大致相等且近似呈对称分布，这一规律与YANG[26]及ROY[22]等的研究结果较为一致;巷道拱顶部位损伤区深度较大，其损伤深度达到了3.8 m，而地板损伤区深度相对较小，为1.7 m。较大范围的损伤破裂导致围岩极易发生挤压变形，应及时进行衬砌支护，从而保证巷道的稳定性。然而，受限于目前的监测探测技术手段，无法获取真实的深部巷道围岩裂纹演化分布规律及开挖损伤区分布特征，因此难以实现上述计算结果的直接对比验证。
图12 深部高应力煤巷围岩开挖损伤区分布特征
Fig.12 Distribution of the excavation damaged zone in surro- unding rock mass of the deep coal roadway
5 结 论
(1)Voronoi随机多边形块体环路的引入，为较好反映深部煤巷围岩体宏观随机非连续面几何分布特征与空间镶嵌结构特性提供了新的思路方法;同时，界面接触计算模型在Voronoi随机多边形环路中的嵌入为复杂多裂纹扩展提供了更多潜在路径，能够在同一计算方法架构下对围岩峰前损伤扩容连续变形以及峰后破裂碎胀非连续变形进行模拟仿真。
(2)在改进的NMM算法平台上发展的基于软化单元法和删除单元法的开挖模拟算法，能够较好地反映深部煤巷开挖卸荷作用和掌子面的空间支撑效应;通过该方法能够较好地对深部煤巷开挖卸荷作用下围岩破裂碎胀行为特征、围岩径向收敛变形特征以及最终形成的EDZ分布特征等进行模拟分析，具有一定的实用性。另外需要指出的是，基于软化单元法的巷道围岩纵向支撑效应模拟是一种简化的二维模拟方法，其中空间约束效应以及卸荷速率与巷道开挖进尺间的定量关系仍有待进一步研究。
(3)由于该方法目前仍处于初始阶段，在进行深部煤巷围岩碎胀大变形模拟计算中，尚未考虑复杂地质条件(例如断层破碎带，复合地层，高地应力、高地温、高渗透压THM耦合作用等)的影响。该改进算法的提出，为复杂工况、地质条件下深部煤巷围岩破裂碎胀挤压大变形的模拟提供了良好基础。
[1] 刘泉声,黄兴,时凯,等.煤矿超千米深部全断面岩石巷道掘进机的提出及关键岩石力学问题[J].煤炭学报,2012,37(12):2006-2013.
LIU Quansheng,HUANG Xing,SHI Kai,et al.Utilization of full face roadway boring machine in coal mines deeper than 1 000 m and the key rock mechanics problems[J].Journal of China Coal Society,2012,37(12):2006-2013.
[2] WANG Jinhua.Development and prospect on fully mechanized mining in Chinese coal mines[J].International Journal of Coal Science & Technology 2014,1(3):253-260.
[3] 卢兴利,刘泉声,苏培芳.考虑扩容碎胀特性的岩石本构模型研究与验证[J].岩石力学与工程学报,2013,32(9):1886-1893.
LU Xingli,LIU Quansheng,SU Peifang.Constitutive model of rocks considering dilatancy-bulking behavior and its calibration[J].Chinese Journal of Rock Mechanics and Engineering,2013,32(9):1886-1893.
[4] 黄兴,刘泉声,刘恺德,等.深部软弱地层TBM掘进围岩变形破坏特性室内试验研究[J].岩石力学与工程学报,2015,34(1):76-92.
HUANG Xing,LIU Quansheng,LIU Kaide,et al.Laboratory study of deformation and failure of soft rock for deep ground tunneling with TBM[J].Chinese Journal of Rock Mechanics and Engineering,2015,34(1):76-92.
[5] 何军.基于数值流形方法的裂隙扩展模拟及其在岩土工程中的应用[D].武汉:武汉大学,2016.
HE Jun.Crack propagation analysis based on numerical manifold method and its application on geotechnical engineering[D].Wuhan:Wuhan University,2016.
[6] FRIES T P,BELYTSCHKO T.The extended/generalized finite element method:An overview of the method and its applications[J].International Journal for Numerical Methods in Engineering,2010,84(3):253-304.
[7] 黄茂松,张治国,王卫东.基于位移控制边界单元法盾构隧道开挖引起分层土体变形分析[J].岩石力学与工程学报,2009,28(12):2544-2553.
HUANG Maosong,ZHANG Zhiguo,WANG Weidong.Response analysis of shield tunneling-induced layered soil deformation based on deformation controlled boundary element method[J].Chinese Journal of Rock Mechanics and Engineering,2009,28(12):2544-2553.
[8] HUANG X,LIU Q S,PENG X X,et al.Mechanism and forecasting model for shield jamming during TBM tunnelling through deep soft ground[J].European Journal of Environmental and Civil Engineering,2017:1-34.
[9] DAMJANAC B,BOARD M,LIN M,et al.Mechanical degradation of emplacement drifts at Yucca Mountain-A modeling case study part II:Lithophysal rock[J].International Journal of Rock Mechanics and Mining Sciences,2007,44(3):368-399.
[10] 江贝,李术才,王琦,等.基于非连续变形分析方法的深部沿空掘巷围岩变形破坏及控制机制对比研究[J].岩土力学,2014,35(8):2353-2360.
JIANG Bei,LI Shucai,WANG Qi,et al.Comparative research on surrounding rock failure deformation and control mechanism in deep roadway driving along goaf based on discontinuous deformation analysis[J].Rock and Soil Mechanics,2014,35(8):2353-2360.
[11] 姚池，姜清辉，邵建富，等.一种模拟岩石破裂的细观数值计算模型[J].岩石力学与工程学报，2013，32(S2):3146-3153.
YAO Chi,JIANG Qinghui,SHAO Jianfu,et al.A mesoscopic numerical model for simulation of rock fracturing[J].Chinese Journal of Rock Mechanics and Engineering,2013，32(S2):3146-3153.
[12] LISJAK A,GRASSELLI G.A review of discrete modeling techniques for fracturing processes in discontinuous rock masses[J].Journal of Rock Mechanics and Geotechnical Engineering,2014,6(4):301-314.
[13] SHI Genghua.Modelling rock joints and blocks by manifold method[A].Proc.of the 32nd U.S.Symp.on Rock Mechanics[C].Sanata Fe.New Mexico,1992:639-648.
[14] 刘泉声,何军,吴月秀.数值流形方法接触检索算法的改进[J].岩石力学与工程学报,2016,35(1):40-49.
LIU Quansheng,HE Jun,WU Yuexiu.An improvement of contact detection algorithm of numerical manifold method[J].Chinese Journal of Rock Mechanics and Engineering,2016,35(1):40-49.
[15] QU X L,WANG Y,FU G Y,et al.Efficiency and accuracy verification of the explicit numerical manifold method for dynamic problems[J].Rock Mechanics and Rock Engineering,2015,48(3):1131-1142.
[16] QU X L,FU G Y,MA G W.An explicit time integration scheme of numerical manifold method[J].Engineering Analysis with Boundary Elements,2014,48:53-62.
[17] CUNDALL P A,STRACK O D L.A discrete numerical model for granular assemblies[J].Géotechnique,1979,29(1):47-65.
[18] 周创兵,陈益峰,姜清辉.岩体表征单元体与岩体力学参数[J].岩土工程学报,2007,29(8):1135-1142.
ZHOU Chuangbing,CHEN Yifeng,JIANG Qinghui.Representative elementary volume and mechanical parameters of fractured rock mass[J].Chinese Journal of Geotechnical Engineering,2007,29(8):1135-1142.
[19] NOROUZI S,BAGHBANAN A,KHANI A.Investigation of grain size effects on micro/macro-mechanical properties of intact rock using voronoi element discrete element method approach[J].Particulate Science and Technology,2013,31(5):507-514.
[20] BAI S Q,TU S H,ZHANG C,et al.Discrete element modeling of progressive failure in a wide coal roadway from water-rich roofs[J].International Journal of Coal Geology,2016,167(Supplement C):215-229.
[21] LAN H X,MARTIN C D,HU B.Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading[J].Journal of Geophysical Research.Solid Earth,2010,115(14).
[22] ROY N,SARKAR R,BHARTI S D.Prediction model for performance evaluation of tunnel excavation in blocky rock mass[J].Int J Geomech,2018,18(1):04017125.
[23] LISJAK A,FIGI D,GRASSELLI G.Fracture development around deep underground excavations:Insights from FDEM modelling[J].Journal of Rock Mechanics and Geotechnical Engineering,2014,6(6):493-505.
[24] HOEK E,MARINOS P.Predicting tunnel squeezing problems in weak heterogeneous rock masses[J].Tunn Tunn Int.,2000,5:45-51.
[25] 黄兴,刘泉声,乔正.朱集矿深井软岩巷道大变形机制及其控制研究[J].岩土力学,2012,33(3):827-834.
HUANG Xing,LIU Quansheng,QIAO Zheng.Research on large deformation mechanism and control method of deep soft roadway in Zhuji coal mine[J].Rock and Soil Mechanics,2012,33(3):827-834.
[26] YANG S Q,CHEN M,JING H W,et al.A case study on large deformation failure mechanism of deep soft rock roadway in Xin’an coal mine,China[J].Engineering Geology,2017,217:89-101.
