0 引 言
在近距离煤层群开采时，由于煤层间开采相互影响作用(尤其煤层间距较小时)，使得下层煤开采所引起的覆岩运移、采动应力演化规律有别于相距较远煤层。近些年来学者们对近距离煤层开采的覆岩裂隙演化规律及应力分布情况进行了大量富有成效的研究，且多集中在细观及宏观方面研究、实验室相似模拟研究，并取得了一系列重要成果。在细观层次方面：CHEN等[1]采用工业CT研究了保护层开采条件下卸载煤体损伤特性，将卸荷围压、固定偏应力和卸荷围压、固定轴向位移进行了比较，发现损伤变量和渗透率显著增加，但后者组合增加效果更加显著；谢和平[2]通过对大量岩石断口的微观观察，发现岩石的破坏主要是由局部损伤拉破坏造成的，并得到了岩石裂隙产生与扩展的判断依据；凌建明[3]采用电镜扫描研究了卸载岩石细观损伤特性，发现岩石试样的载荷-位移曲线与其损伤演化的各个阶段有对应关系，在荷载适当情况下，岩石细观裂纹的扩展具有显著的时效特征；钱鸣高等提出了关键层理论[4]及“O”形圈理论[5]；刘天泉等[6]提出了“竖三带”及“横三区”理论；康永华等[7]研究了覆岩性质对“两带”高度的影响，研究发现覆岩的移动和破坏是由下而上逐层发展的。实验室相似模拟方面研究：姜福兴等[8]研究了采动应力场与采场覆岩空间破裂的关系，发现采动将引起覆岩发生4种类型破裂形式，并基于此提出了由微震事件分布推断区段煤柱稳定性的方法；文献[9]通过实验室监测研究了裂隙的发育程度，揭示了覆岩采动裂隙的两阶段发展规律及“O”形圈分布特征；文献[10-13]研究了采动覆岩裂隙场的时空演化规律，并基于此提出顶底板卸压瓦斯抽采工程技术方法及瓦斯和煤与瓦斯共采工程技术体系。
随着计算机性能的快速发展，在地下煤矿开采的前期设计、支护方案选取及反演分析等方面，数值模拟技术是一种重要的研究手段。鉴于数值方法可将岩体、节理面及其他非均质结构面尽可能精细化考虑到模型中，从而研究其变形及力学性质，模拟结果可更接近于实际问题。但耗时巨大的计算量制约着精细化模拟的发展，如PFC及UDEC/3DEC中grain based model模拟手段[14-17]或是FEMDEM模拟方法目前都限于二维计算和实验室岩块尺度的模拟[18]。
针对地下采煤的大尺度开挖问题，为提高模拟结果的可靠性和合理减少模拟的计算时间，采用离散元DEM模拟软件3DEC及BBM(Bonded Block Model)黏连块体模型对近距离煤层开挖的采动裂隙场分布及演化进行了模拟研究，并探讨了该方法在煤矿开采实践中应用的可行性，且研究成果可用来加深对煤层开采产生的应力场和裂隙场的理解，也可对煤矿开采实践提供一定的参考，如顶底板破坏深度的估量和煤层群瓦斯抽放钻孔的布置等。
1 BBM-DEM模拟方法
BBM(Bonded Block Model)黏连块体模型是在商业化离散元(DEM)模拟软件3DEC中新开发的块体模型。BBM 块体单元可采用多面体几何形状来表征岩体单元，如Voronoi块体单元、四面体单元等。借助于BBM，可实现模拟裂隙在岩体中的萌生、发育及扩展过程。通过研究接触面的破坏(拉伸破坏、剪切破坏及拉-剪破坏)机理，可加深对岩体动态破坏演化的理解。同时，在煤层开挖模拟过程中，采动裂隙场在岩体中的分布和演化规律可实现更精确及直观的模拟。
与传统规则块体单元的建模过程相比，BBM 块体模型的建模顺序有所不同，结合3DEC及BBM模型概念，BBM-DEM的模拟方法如图1所示。
1)建立常规的3DEC计算模型。如图1a所示，将模型尺寸简化为2 m× 1 m× 3 m, 模拟3层岩体，分别建立各岩层，即直接底，煤层和直接顶，每层高度1 m。
2)在3DEC中，通过gen edge命令将各岩层划分为不同大小的四面体单元。如图1b所示，采用2种(0.1 m和0.2 m)不同的单元尺寸(四面体边长)将图1a模型离散成不同数量的四面体单元。
3)利用Fish把步骤2中的四面体单元的节点坐标及对应各面的拓扑结构导出，从而实现四面体单元成为独立的块体。
4)利用call命令将步骤3中生成的四面体块体重新读入3DEC，从而生成BBM模型。对应如图1b的2种单元尺寸，如图1d所示，分别有36 371个和4 666个四面体块体组成了计算模型。进而对四面体块体再次进行四面体单元的离散划分，从而进行模拟计算。
需要说明的是，当计算机计算能力足够大时，每个独立BBM块体单元的尺寸可以细化到细观尺度(如1 mm)，此时块体尺寸已小到可以模拟矿物颗粒尺度(类似grain based model)，从而精细化的揭示岩石或岩体的破坏变形规律。即提及的BBM-DEM精细化模拟方法，特指岩体块体单元及其接触面的精细化模拟。针对大尺度的煤层采动裂隙场的模拟，目前计算机性能还无法满足此精细化水平。因此，在BBM-DEM模拟方法中，四面体块体大小具有较大尺度选取范围(如微观、细观及宏观尺寸)，取决于处理的工程问题及计算机性能。为了说明此模拟方法的可行性和适用性，对基于近距离煤层保护层开采的模拟算例进行说明。
图1 BBM-DEM 模拟方法
Fig.1 Modeling approach of BBM-DEM
2 近距离煤层保护层开挖模拟
2.1 模拟区域地质概况
沙曲煤矿2号、3+4号及5号煤层为近距离煤层。北二采区钻孔柱状图及各煤层空间分布如图2[19]和图3所示。
图2 各煤层及工作面空间分布关系
Fig.2 Spatial location of coal seams and working faces
由图2及图3可见，2号煤层位于3+4号及5号煤层上方，埋深约500 m。经鉴定，2号煤层属于无突出危险性煤层，而3+4号、5号煤层均属于突出危险性煤层。故在实际开采中，2号煤层作为保护层煤层首先开采。22201工作面是2号煤层首个工作面，地质构造相对简单，平均倾角4°，平均厚度1.46 m。2号煤层保护层的开采，对下伏的3+4号，5号煤层起到了卸压保护作用[5,7,9-11]。24208工作面为3号和4号煤层合并开采的首个工作面，平均倾角4°，平均厚度4.37 m。各岩层厚度及岩性描述如图3所示。
2.2 BBM-DEM数值模型建立
采用BBM-DEM模拟方法对沙曲煤矿2号煤层保护层工作面的采动裂隙场进行数值分析，从而揭示采动裂隙场的演化规律及对下伏被保护层的卸压保护作用，说明该方法的可行性。
基于沙曲煤矿北二采区地质条件(图2和图3[19])，合理简化模型后(如不考虑岩层倾角的影响)，按照图1所示的建模方法，建立BBM-DEM计算模型。图4给出了数值计算模型的构建、边界条件、岩层情况、及四面体块体所对应的局部放大图。为提高计算效率，模型尺寸为220 m × 1 m × 200 m(图4)，在计算中，将模型所有节点沿Y方向的速度固定，从而将三维问题简化为二维平面应变问题。在建模过程中，对于模型顶和底部的中砂岩，其对应的四面体尺寸为4 m, 所研究煤层及其邻近岩层，对应尺寸1 m的四面体块体。模型共考虑了18层煤岩体，共有81 052个四面体块体。如将模型建为三维模型(如沿Y方向220 m)，则有约18 000 000个四面体块体，但目前还无法在普通的工作站或者个人计算机实现如此大的计算量。
图3 沙曲煤矿北二采区钻孔柱状
Fig.3 Borehole log of the northern second mining area of Shaqu Coal Mine
如图4所示，模型左、右边界为位移边界条件，限制各块体的水平方向运动。模型底部为固定边界条件。模型上部边界为应力边界条件，载荷为上覆岩层的自重。由于2号煤层埋深约500 m，模型顶部埋深约450 m，基于康红普等[20]关于山西煤矿矿井原岩应力场的研究，模型顶部施加的原岩应力：垂直主应力sv=10.3 MPa,最大水平主应力sH=7.73 MPa, 最小水平主应力sh=3.59 MPa。 数值分析的流程为：BBM-DEM模型构建、合理边界条件施加、块体及接触面变形及强度参数选取、原岩应力场施加、应力场初始平衡、开挖2号保护层煤层(开挖计算循环为10 m，左右预留宽50 m煤柱)，模拟结束，分析模拟计算结果。
2.3 岩块及接触面变形及强度参数选取
由图4可知, 各岩层在BBM-DEM模拟方法中都是由黏连在一起的四面体块体所组成，四面体块体及其接触面的变形及强度对整个煤层采动裂隙场演化模拟的准确性和有效性有重要的影响。
在上述分析流程中，煤岩体的变形及力学参数选取是非常重要的科学研究课题之一。目前，岩体变形及强度参数的估算主要有直接方法和间接方法。在直接方法中，小尺度实验室岩块力学试验和大尺度原位试验是主要手段。在间接方法中，数学解析解、反演分析、围岩经验分类系统及数值模拟法是其主要手段。每种估算方法都有其优劣势，如周期、费用、选取地点差异性、可行性设计阶段适用性、节理裂隙系统等，大量研究仍在持续进行中，在此不做详述。
图4 BBM-DEM数值模拟模型
Fig.4 Numerical model of BBM-DEM
所模拟的四面体块体尺寸近1 m，不能直接赋予实验室岩石变形及强度参数，需要将实验室得到的参数进行折减。一般情况下，岩体变形及力学参数在数值模拟中是需要校准和标定的。如第1节所述，若四面体块体大小能精细到实验室岩块尺度(如0.05 m)，则岩块的长期强度参数可直接赋予模型的块体。模型数值分析的目的只是为说明BBM-DEM模拟方法的可行性，确保估算的岩体参数在合理范围内，而不在于如何精确地估算岩体参数。故以实验室岩石物理力学参数为基础，岩石实验室试验数据参照阎峰等[21]的建矿地质报告，根据Hoek-Brown强度准侧和GSI岩体分类系统对各岩层的岩体变形及强度参数进行估算(利用Roclab软件)。表1给出了在模拟中采用的各层岩体的变形及强度参数。
在四面体块体的变形及强度参数确定之后，接触面变形及强度力学参数的估算是另一难点。在前人研究的基础上，并结合笔者数值分析的经验，采用以下方法估算接触面的变形及强度参数。
1)接触面法向刚度(JKn)：利用公式JKn=n×max(K+0.75G)/Δz进行预估算[20]，确保法向刚度的取值合理。式中：n取值范围为1～10；K和G分别为体积模量和剪切模量；Δz为接触面之间最小距离。为简化起见，在模拟中取n=5，Δz可通过块体大小计算。
2)接触面切向刚度(JKs)：基于文献[22-23]的研究，利用公式Ks/Kn =G/E进行估算，式中：E为岩体弹性模量；Ks为块体切向刚度；Kn为块体法向刚度。
3)接触面黏聚力(JC)和内摩擦角(Jf)：基于文献[24]利用UDEC的研究结果，煤岩块接触面JC和Jf约为煤岩体黏聚力和内摩擦角的120%和50%。为简化分析，考虑到模型中四面体块体的尺寸(1 m)，假定接触面JC为对应煤岩体黏聚力的1.5倍，假定接触面Jf为对应煤岩体内摩擦角的0.8倍。
4)接触面抗拉强度(JT)：假定尺寸4 m和1 m四面体块体接触面的抗拉强度约为岩石单轴抗拉强度的30%和50%。如实验室试验测得中粒砂岩的单轴抗拉强度均值为4.5 MPa，从而接触面编号1和接触面编号2的JT估算值分别为1.35、2.25 MPa。此外，接触面的残余黏聚力、残余抗拉强度、及残余摩擦角都简化为0。
表1 各层岩估算的变形及强度力学参数
Table 1 Estimation of deformation and mechanical strength parameters of rock masses
岩层层号密度/( kg·m-3)体积模量/GPa剪切模量/GPa黏聚力/MPa内摩擦角/(°)剪胀角/(°)抗拉强度/MPa中砂岩1,202 5007.14.705.242123.8中砂岩2,7,11,15,192 50014.29.4010.246185.6细粒砂岩3,82 6207.75.305.638203.6砂质泥岩4,142 6007.23.902.836163.2煤岩5,13,171 3502.51.151.828180.8粉砂岩6,102 2003.82.102.231203.4黑色泥岩9,12,162 6508.33.602.234162.8黏土泥岩182 6808.43.402.636183.8
在表2中，共有8种接触面参数，对应着表1中的不同岩层内部的接触面。同时，在2种岩层之间的接触面，则为其赋予较弱岩层的接触面参数，如和煤层(煤层所对应接触面编号为5)接触的直接顶和直接底的接触界面，对其赋予煤层接触面的变形及力学参数。需要指出的是，表1和表2的参数取值仅是为模型赋予合理的模拟计算参数。另外，值得提出的是，可通过现场监测的位移或应力来进行模型参数的反演分析，从而增加所取模型参数的可靠性。
基于表1及表2岩体及接触面的变形及强度估算参数，在BBM-DEM数值模拟中，四面体块体采用摩尔-库仑破坏准则，接触面采用面接触的库仑滑移破坏准则进行模拟。以下将给出保护层开挖岩层的运动规律及开挖卸荷效果的数值模拟初步结果。
3 保护层开挖模拟结果及讨论
3.1 覆岩运动及采动裂隙演化分析
在模拟2号保护层煤层开挖过程中，以开挖长度10 m为一计算循环，每次迭代60 000步。共分12次开挖。模型左右边界各预留宽50 m煤柱，开切眼在X=50 m处。图5显示了覆岩层随工作面推进的运动演化规律及采动裂隙场分布情况。需指出的是，为了更清晰观察采动裂隙场的情况，图5只显示了模型的局部范围(X:48～122 m;Z:80～130 m)。图5的每个推进距离，有3种形式的图组合在一起说明在该推进距离时的裂隙及岩层演化规律。最上方的图对应的是各个岩层运移图，中间图为对应虚线矩形框中放大的裂隙图，下方的图对应矩形框中放大的四面体块体单元和接触面破坏模式图。
表2 各接触面估算的变形及强度力学参数
Table 2 Estimation of deformation andmechanical strength parameters for contact properties
接触面编号JKn/GPaJKs/GPaJf/(°)JC/MPaJT/MPa119.47.833.67.81.352154.063.036.815.32.25386.035.030.48.41.604100.040.028.84.21.35524.010.022.42.70.41637.015.024.83.21.40798.041.027.23.31.20876.030.028.83.91.10
注：接触面编号自上而下依次对应表1的8种岩层。
由图5a可知，自工作面开切眼起，当工作面推至40 m时，直接顶产生大量裂隙，形成悬臂结构，且离层加剧。此时，基本顶产生离层及断裂，造成工作面的初次来压。同时，每次的模拟开挖都会对2号煤层的覆岩结构产生影响，特别是对直接顶的泥岩以及上方的细粒砂岩造成较大影响。随着直接顶悬露面积的增加，可以清晰观测到直接顶出现了4次倾斜破断，而非竖直方向断裂。这与模型选用的BBM四面体块体有关。通过对块体单元接触面的破坏模式分析，可发现沿着四面体块体接触面的拉破坏(蓝色十字)是主要的破坏模式，而剪切破坏主要发生在实体煤壁侧的接触面，主要是实体煤前方的应力集中所造成的。当工作面推进到70 m时，由图5b可知，沿开采方向20～50 m的直接顶发生大面积垮落，直接顶垮落高度接近煤层开挖的高度。同时，细粒砂岩和中砂岩层出现大量横向及斜向裂隙，出现下沉及断裂，此时对应着周期来压，周期来压步距约为70 m。同时，在直接顶大面积垮落和触底过程中，会出现局部非充分接触和全接触情况，如图5b中(2)所示。由于岩体中岩块的非均质性，在其破裂离层过程中可能会发生局部触底现象。利用BBM块体来模拟，能较为真实地反映直接顶触底的过程。另外，大量的拉破坏发生采空区上方的各岩层中，同时接触面的剪切滑移破坏出现在工作面前方的实体煤岩体中。当工作面推进至100 m时(图5c)，可发现直接顶泥岩(距开切眼20 ～ 80 m)大面积垮落到煤层底板，且上方细粒砂岩和中砂岩岩层产生大量裂隙，逐渐压实采空区的，产生弯曲下沉。图5d显示了当工作面推进到120 m时的覆岩及裂隙分布情况。可发现此时直接顶除两端之外几乎全部垮落，上覆岩层大量裂隙发育，有明显下沉。
图5 工作面推进中覆岩运动规律及裂隙分布情况
Fig.5 Overburden movement as advancing of working face and fracture distribution
由此可见，该BBM-DEM模拟方法可模拟在工作面推进过程中上覆岩层的运动规律，且可模覆岩的动态压实采空区过程。由接触面的破坏机理可知采动裂隙的分布及发育程度。
3.2 覆围岩位移场分析
当模型开挖到2号保护层煤120 m时，图6给出了其垂直位移场的分布。实际模拟的垂直位移分量最小值和最大值分别为 -1.38 m和0.71 m，未在图6中标出。为得到清晰的等值图分布，垂直位移所取范围为 -0.8～0.6 m，正值代表向上运动，负值代表下向运动。此外，为了监测应力及位移的变化，一系列的测点(10 m间隔)布置在模型的5条测线中，如图6所示。
1—测线1，位于2号煤层顶板；2—测线2，位于2号煤层底板；3—测线3，位于中砂岩；4—测线4，位于3+4号煤层；5—测线5，位于5号煤层
图6 2号煤层开挖120 m时垂直位移场分布
Fig.6 Vertical displacement contour when No.2 coal seam advanced 120 m
由图6可得，沿着2号煤层开挖方向，垂直位移在2号煤层上覆岩层中呈非均匀分布，直接顶的大面积垮落和基本顶的离层断裂造成了较大的垂直位移，且其分布与直接顶的垮落方式及地点紧密相关。
位移场分布基本反映出了采空区上方的垮落带、裂隙带及弯曲下沉带的特征。同时2号煤层的开挖对下伏的各个岩层，特别是被保护层3+4号煤层及5号煤层也有一定程度的影响。为进一步了解保护层开采对下伏岩层的影响，图7给出了在工作面推进120 m时，垂直位移分量在所选的5个测线上的分布。
由图7可得，保护层的开采对下伏各岩层的位移分布均有影响。直接顶(测线1)的垂直位移受2号保护煤层开挖的影响最大，部分测点的垂直位移超过1 m。而直接底(测线2)在局部区域也有较大的垂直位移，可达1 m左右，对应潜在的底鼓现象。由测线3—5可得，处于采空区中部压实区以下的垂直位移要大于开切眼和工作面附近的位移。另外，还可看出保护层以下的岩层，其垂直位移呈现由上至下减少的趋势。这意味着在近距离煤层中，保护层的开采对附近的煤岩层有着重要的影响作用。
图7 2号煤开挖120 m时5条测线上垂直位移分布
Fig.7 Vertical displacement distribution along five monitoring lines when No.2 coal seam advanced 120 m
3.3 围岩应力场及保护层开采卸压分析
由于煤层的开挖扰动，原岩应力会在开挖附近区域重新分布，形成应力集中区和降低区。为分析保护层煤层开挖对下伏被保护层煤层的卸压作用，图8给出了当2号煤层开挖至100 m时的垂直应力场。为清晰展现垂直应力的分布，图8中垂直应力的取值范围是 -40～2 MPa，而在模拟中，局部最大垂直应力可达到200 MPa, 未在图8中标出。另外，类似于图6，在不同的岩层中布置了5条测线用于检测垂直应力的分布情况。
1—测线1，位于2号煤层顶板；2—测线2，位于2号煤层底板；3—测线3，位于中砂岩；4—测线4，位于3+4号煤层；5—测线5，位于5号煤层
图8 2号煤层开挖100 m时垂直应力场分布
Fig.8 Vertical stress contour when No.2 coal seam advanced 100 m
由图8可知，在开切眼后方和工作面前方的实体煤岩体和顶板岩层中形成了应力集中区，而在采空区及上、下方的岩层区域形成了应力降低区。同时，在2号煤层采空区附近的区域大体都处在应力降低区。
但在局部采空区重新压实的区域，其应力可恢复到较高的水平(如10 MPa)。并且在采空区中段压实的下方各岩层中都出现了垂直应力的局部回升现象。这是因为在模拟中，BBM块体在力的作用下会发生张拉、滑移及旋转，能更加有效地模拟岩层断裂和动态失稳的过程。为进一步说明在近距离煤层中，保护层开采对下伏被保护层的卸压作用，图9给出了在工作面推进100 m时，垂直应力在所选5条测线上的分布情况。
图9 2号煤层开挖100 m时5条测线上垂直应力分布
Fig.9 Vertical stress distribution along five monitoring lines when No.2 coal seam advanced 100 m
由图9可得，保护层开采对近距离煤岩层的卸压作用。2号煤层未开挖时，测线1—5上的初始垂直应力水平依次为12.36、12.4、12.7、13.0、13.2 MPa，测线1—5对应的最大应力集中系数依次为3.4、3.8、1.9、2.6、2.6。由于2号煤层直接顶和基本顶的大面积垮落和应变能释放，在采空区区域附近的垂直应力除了偶尔几个测点的应力回升都发生了很大程度的卸压。
由测线1和2可得，开切眼后方(X<50 m)和工作面前方(X>150 m)都出现了不同程度的应力集中现象。而在采空区区域(50 m<X<150 m)内出现了大范围的卸压。其他3条测线(3—5)上的垂直应力具有同样的变化趋势。另外，齐庆新等[19]利用不同的理论及经验公式对2号煤层的底板破坏深度进行计算，发现其理论计算的最小值就已经接近5号煤的顶板。这进一步说明2号煤层保护层的开采对下伏煤层具有明显的卸压效果，也说明BBM-DEM可用来分析煤层的采动裂隙场。
4 结 论
1)基于BBM-DEM模拟手段，给出了针对煤矿开采的BBM-DEM建模过程及模拟方法。探讨了该方法的优势及局限性，以及在采矿工程中应用的可能性。
2)采用BBM-DEM模拟方法对沙曲煤矿2号煤层保护层工作面的开采进行数值模拟分析，以此说明此方法的可行性。从上覆岩层的运动及裂隙演化规律、围岩位移场分布和围岩应力场分布3个方面分析了保护层开采对下伏岩层的卸压保护作用。并且发现采空区应力在压实区非均匀分布，具有压实局部化现象。模拟的卸压效果与其他学者的研究结果具有一致性。
3)BBM-DEM模拟技术可用来模拟煤层开挖的过程，特别是对近距离高瓦斯煤层群的采动裂隙场模拟。随着计算机计算性能的快速发展，此方法具有可精细化的潜力。且模拟结果可加深对煤层开采产生的应力场和裂隙场的理解，也可对煤矿开采实践提供一定的参考，如顶底板破坏深度的估量和煤层群瓦斯抽放钻孔的布置等。
[1] CHEN Haidong, CHENG Yuanping, ZHOU Hongxing, et al. Damage and permeability development in coal during unloading[J]. Rock Mechanics and Rock Engineering, 2013, 46(6): 1377-1390.
[2] 谢和平. 岩石材料的局部损作拉破坏[J]. 岩石力学与工程学报, 1988, 7(2), 147-154.
XIE Heping. Tensile failure of local damage of rock material[J]. Chinese Journal Rock Mechanics and Engineering, 1988, 7(2), 147-154.
[3] 凌建明. 压缩荷载条件下岩石细观损伤特征的研究[J]. 同济大学学报, 1993, 21(2):219-226.
LING Jianming. Study of mesocrack damage of rock under compression condition[J]. Journal of Tongji University, 1993, 21(2):219-226.
[4] 钱鸣高, 缪协兴, 许家林，等. 岩层控制的关键层理论[M]. 徐州:中国矿业大学出版, 2000.
[5] 许家林, 钱鸣高. 覆岩采动裂隙分布特征的研究[J]. 矿山压力与顶板管理, 1997, 14(3/4): 210-212.
XU Jialin, QIAN Minggao. Study of feature of distribution of overlaying stratum fractures caused by mining[J]. Ground Pressure and Strata Control, 1997, 14(3/4): 210-212.
[6] 刘天泉, 仲维林, 焦传武，等. 煤矿地表移动与覆岩破坏规律及应用[M]. 北京:煤炭工业出版社, 1981.
[7] 康永华, 赵国玺. 覆岩性质对“两带”高度的影响[J]. 煤矿开采, 1998，3(1):52-54.
KANG Yonghua, ZHAO Guoxi. The effect of overburden strata characteristics on heights of the fractured and caved zones[J]. Coal Mining Technology, 1998,3(1):52- 54.
[8] 姜福兴, XUN Luo, 杨淑华. 采场覆岩空间破裂与采动应力场的微震探测研究[J]. 岩土工程学报,2003, 25(1):23-25.
JIANG Fuxing, XUN Luo, YANG Shuhua. Study on microseismic monitoring for spatial structure of overlying strata and mining pressure field in longwall face[J]. Chinese Journal of Geotechqiue Engineering, 2003, 25(1):23-25.
[9] 钱鸣高,许家林. 覆岩采动裂隙分布的“O”形圈特征研究[J].煤炭学报,1998,23(5):466-469.
QIAN Mingao,XU Jialin. Study on the “O-shape” circle distribution characters of mining-induced fractures in the overlaying strata[J]. Journal of China Coal Society, 1998,23(5):466-469.
[10] 李树刚, 石平五, 钱鸣高. 覆岩采动裂隙椭抛带动态分布特征研究[J]. 矿山压力与顶板管理, 1999,14(3):44-46.
LI Shugang, SHI Pingwu, QIAN Minggao. Study on the dynamic distribution features of mining fissure elliptic paraboloid zone[J]. Ground Pressure and Strata Control, 1999, 14(3): 44-46.
[11] 袁 亮. 低透高瓦斯煤层群安全开采关键技术研究[J]. 岩石力学与工程学报, 2008, 27(7):1370-7378.
YUAN Liang. Key technique of safe mining in low permeability and methane-rich seam group[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(7): 1370-7378.
[12] 袁 亮. 卸压开采抽采瓦斯理论及煤与瓦斯共采技术体系[J]. 煤炭学报, 2009, 34(1): 1-8.
YUAN Liang. Theory of pressure-relieved gas extraction and technique system of integrated coal production and gas extraction[J]. Journal of China Coal Society, 2009, 34(1): 1-8.
[13] 刘泽功. 卸压瓦斯储集与采场围岩裂隙演化关系研究[D].北京：中国科学技术大学, 2004.
[14] WANG Xin, CAI Ming. Modeling of brittle rock failure considering inter-and intra-grain contact failures[J].Computers and Geotechnics, 2018, 101：224-244.
[15] CHO N, MARTIN CD, SEGO D. A clumped particle model for rock[J]. Int J Rock Mech Min Sci，2007,44(7):997-1010.
[16] GAO F,STEAD D, ELMO D. Numerical simulation of microstructure of brittle rock using a grain-breakable distinct element grain-based model[J]. Comput Geotech，2016，78:203-17.
[17] POTYONDY D. A grain-based model for rock: approaching the true microstructure[C]//In:Proceedings of rock mechanics in the Nordic Countries,2010.
[18] XIANG J,MUNJIZA A, LATHAM JP,et al. On the validation of DEM and FEM/DEM models in 2D and 3D[J]. Eng Comput,2009,26:673-687.
[19] 齐庆新,季文博,元继宏,等. 底板贯穿型裂隙现场实测及其对瓦斯抽采的影[J]. 煤炭学报, 2014,39(8):1552 -1558.
QI Qingxin, JI Wenbo, YUAN Jihong,et al. In-situ observation on the floor-connected crack field and its effects on methane extraction[J].Journal of China Coal Society, 2014,39(8):1552-1558.
[20] 康红普，林 健，颜立新, 等.山西煤矿矿区井下地应力场分布特征研究[J].地球物理学报，2009, 52(7):1782-1792.
KANG Hongpu, LIN Jian, YAN Lixin,et al. Study on characteristics of underground in-situ stress distribution in Shanxi Coal mining fields[J].Chinese Journal of Geophysicals，2009, 52(7):1782-1792.
[21] 闫 峰，李建芳，李 刚，等.华晋焦煤有限责任公司沙曲二号矿井建矿地质报告[R].太原：山西省煤炭地质144勘查院,2015.
[22] FARAHMAND K, DIEDERICHS MS. A calibrated synthetic rock mass (SRM) model for simulating crack growth in granitic rock considering grain scale heterogeneity of polycrystalline rock[C]//The 49th US rock mechancis/geomechanics symposium,2015.
[23] KAZERANI T,YANG Z Y,ZHAO J.A discrete element code model for predicting shear strength and degradation of rock joint by using compressive and tensile test data[J].Rock mechanics and rock engineering，2012, 45(5),695-709.
[24] GAO Fuqiang. Simulation of failure mechanisms around underground coal mine opening using discrete element modeling [D].Vancouver:Simon Fraser University，2013.
