采动应力与采空区压实承载耦合分析方法研究

蒋力帅，武泉森，李小裕，丁 楠

(1.山东科技大学 矿山灾害预防与控制国家重点实验室培育基地，山东 青岛 266590; 2.山东科技大学 矿业与安全工程学院，山东 青岛 266590)

采动应力与采空区压实承载耦合分析方法研究

蒋力帅1,2，武泉森2，李小裕2，丁 楠2

(1.山东科技大学 矿山灾害预防与控制国家重点实验室培育基地，山东 青岛 266590; 2.山东科技大学 矿业与安全工程学院，山东 青岛 266590)

采空区垮落带在覆岩沉降作用下的应力恢复对采动应力场有显著影响。基于围岩应力与采空区压实承载的耦合作用，研究了垮落带岩体压实承载的力学特性及参数反演方法，分析了裂隙带岩体的拉伸劣化力学特性，提出了采动应力场与采空区压实承载耦合分析方法，通过实例仿真模拟分析进行了验证。研究表明：垮落带渐进压实承载的力学特性是指数型表达式，压缩刚度随压缩应变呈非线性增长。由双屈服本构模型反演得到与理论解拟合度较好的垮落带力学参数，实现了垮落带压实承载特性的仿真模拟。裂隙岩体劣化系数表达式反映了模量劣化特性，应变软化本构模型增加模量劣化算法组成岩体劣化本构模型，实现了应变软化和裂隙岩体模量劣化的动态分析，较真实严谨的反映裂隙带岩体的力学特性。基于垮落带与裂隙带特征的仿真模拟，实现了采动应力场与采空区覆岩“两带”特征的耦合分析。通过工作面支承压力场演化特征的实例仿真分析，验证了采动应力场与采空区压实承载耦合分析方法的可行性和正确性，提高了采动应力场和围岩采动响应研究的严谨性和可靠性。

采动应力;垮落带压实;裂隙带劣化;数值模拟;耦合分析

长壁工作面垮落法管理顶板是最为普遍的采煤方法。随着工作面推进，采场围岩应力分布不断调整及演化，采空区直接顶垮落、基本顶及其上部岩层破断运移、高位岩层缓慢沉降等破坏及移动现象随之出现，垮落的松散岩体充填采空区形成垮落带，上覆岩层破断运移形成裂隙带，高位岩层直至地表形成弯曲下沉带。工作面采空区冒落岩体在上覆岩层沉降作用下，由松散状态逐渐压实，支撑能力不断提高，最终形成支撑上覆岩层的承载体。这一动态力学过程实质上是围岩应力与覆岩运移及垮落带逐渐压实的耦合作用过程，围岩应力集中引起覆岩运移，覆岩运移又导致围岩应力分布演化。其中，支承压力分布是采动应力场的核心，是沿空巷道布置、回采巷道围岩控制、冲击地压防治、煤与瓦斯突出防治等实际问题的理论基础和重要依据[1-3]。

国内外学者对采动应力场的分布特征都十分关注，基于连续介质的数值模拟方法被广泛应用。在数值模拟计算中，通常以开挖方式处理工作面采空区，并未考虑垮落带的渐进压实特性及裂隙带特征的影响;或者将采出空间或垮落带赋予软弱岩性，并不符合垮落带压实过程的力学特性;或者选择适合垮落带特性的本构模型，但力学参数选取缺少依据。PAPPAS和C.MARK[4]指出采空区破碎松散的冒落岩体压实力学特性对探究长壁工作面开采引起的围岩扰动具有重要意义，并通过试验研究了破碎松散岩体的压实力学特性;张振南、马占国等[5-7]通过松散岩块的压实试验，分析了变形模量与轴向应力、轴向应变的关系;YAVUZ[8]以英国某煤矿的开采沉陷为背景，研究了采空区内岩体压实与应力恢复主要与采高、埋深等参数有关，并采用FLAC数值计算软件反演了采空区压实力学特性;白庆升、屠世浩等[9]基于采空区压实理论进行了FLAC3D反演分析。近年来，一些学者采用连续介质数值模拟方法，对采空区进行等效模拟或近似等效处理。目前，相关研究主要集中在采空区内冒落岩体的力学特性方面，数值模拟中未能正确反映垮落带岩体的压实力学特性，而采动影响下裂隙带岩体的力学特征及其对采动应力场的影响却鲜有探索。

垮落带松散岩体的压实与承载力学特性，是采空区应力恢复与分布规律的基础，明显影响着采动应力场特别是支承压力分布。覆岩裂隙带范围较大、裂隙纵横，裂隙带岩层的采动劣化对采动应力场的影响是不应忽略的。因此，在前人有关研究成果的基础上，基于采动应力场与覆岩破坏特征的耦合分析，综合研究采动影响下垮落带冒落岩体的压实承载特性和裂隙带岩体的拉伸劣化特性，通过对FLAC3D中本构模型的二次开发，研究其在连续介质数值模拟中的耦合再现技术，进而提出较真实严谨的采动应力场仿真研究方法，并结合工程实例验证该方法的可行性及效果。研究工作对支承压力分布、区段煤柱稳定性、沿空巷道布置、回采巷道采动影响与控制设计、冲击地压危险评估、煤与瓦斯突出防治、岩层与地表移动等方面，具有重要的科学意义及应用价值。

1 垮落带压实承载特性及参数反演方法

1.1 采动覆岩破坏高度

长壁工作面开采过程中，上覆岩层将不断发生垮落及破断运移，在采空区形成覆岩“三带”结构，其中垮落带和裂隙带直接影响着采动应力分布和覆岩运移，“两带”高度及其力学特性是采动应力场仿真中模拟极其重要的组成部分，对水体下采煤及下解放层开采等具有重要意义[10]。

采用钻孔双端堵水器、钻孔冲洗液漏失量等方法，通过现场实测可得到具体工程条件下的覆岩“两带”高度。在大量现场实测分析的基础上，国内外学者研究提出了工作面垮落带和裂隙带高度的计算式。Palchik[11]结合乌克兰顿涅茨克煤田的实践经验与前人理论分析成果，认为垮落带高度与采高及碎胀系数有关，即：

式中，Hc为垮落带高度;h为采高;b为碎胀系数。

BAI,KENDORSKI等[12]在中、美两国大量现场实测数据的基础上，通过统计回归分析得到垮落带和裂隙带高度的回归计算式，见式(2)。该计算方法包含不同岩体强度对“两带”高度的影响，精确性较高，因而得到广泛应用。

式中，σc为单轴抗压强度;Hf为裂隙带高度;c1,c2,c3,c4为岩层强度系数，见表1。

表1垮落带和裂隙带岩层强度系数
Table1Coefficientsforheightofcavingandfracturedzone

岩石类别σc/MPac1c2c3c4坚硬>402116122中硬20～4047191636软弱<206232315

我国“三下”采煤规程中提供了覆岩“两带”高度的计算方法。金太等[13]综合分析了兖州矿区的现场实测资料，得到缓倾斜厚煤层分层开采和综放开采“两带”高度的变化规律与计算式。根据现场实测结果，或针对覆岩条件采用上述有关计算式，可得到比较可靠的“两带”高度，为建立数值计算模型提供依据。

1.2 垮落带压实承载过程的力学特性

随着工作面推进，直接顶板发生冒落及碎胀，充填采空区的垮落带岩体最初呈现松散状态，在裂隙带及弯曲下沉带的沉降荷载作用下，逐渐压实同时支撑上覆运移岩层。垮落带岩体的渐进压实过程及其承载力学特性，对采空区与煤体应力状态有着极其重要的影响。对松散破碎岩体压实演变过程中的岩体力学响应特性，国内外学者进行了理论和试验研究[5-9,14]。由SALAMON[14]提出的破碎岩体压缩过程的应力-应变关系见式(3)，被国内外研究人员广泛认可并应用[2,8,15-17]。

式中，σ为垮落带岩体所受垂直应力;E0为岩体初始正切模量;ε为垂直应力作用下垮落带岩体的应变量;εm为岩体碎胀后最大应变。

其中参数E0主要受岩体碎胀系数b及σc的影响。YAVUZ[8]基于前人大量不同岩石单轴压缩试验数据，通过三维回归分析得到E0的表达式：

岩体碎胀后达到的最大应变εm可由碎胀系数b得到

联立式(1),(2)可得碎胀系数与岩体强度参数及采高的关系：

由式(3),(4),(5)得到垮落带岩体压实特性的表达式：

表达式(7)反映了垮落带岩体渐进压实-支撑过程的力学特性。由式(7)及式(3)～(5)可知，垮落带岩体压实过程中的应力-应变关系由碎胀系数b和单轴抗压强度σc确定，其中碎胀系数对压实力学特性的影响更为显著。

岩体单轴抗压强度为30 MPa时碎胀系数b对垮落带岩体刚度的影响如图1所示。可见，垮落带岩体的压缩刚度随着压缩应变的增长而明显提高，岩体碎胀系数对碎胀后的最大压缩应变及压缩刚度影响很大。高碎胀系数岩体的初期压缩刚度上升趋势平缓，后期上升幅度很大，压缩应变较大;低碎胀系数岩体的压缩刚度上升较快，压缩应变较小。

图1 碎胀系数对垮落带岩体压缩刚度的影响Fig.1 Effect of bulking factor on the compression stiffness of caved rock

碎胀系数为1.3和1.5时单轴抗压强度对垮落带岩体刚度的影响如图2所示。可见，在碎胀系数相同的前提下，高强度垮落岩体的压缩刚度明显大于低强度岩体，相同应力状态下高强度垮落岩体的应变总是小于低强度岩体。需要注意的是，垮落岩体的体积永远无法压实至初始完整体积，仍保留残余碎胀。

图2 单轴抗压强度对垮落带岩体压缩刚度的影响Fig.2 Effect of uniaxial compressive strength on the compression stiffness of caved rock

上述研究得到的垮落带岩体应力-应变关系式，揭示了垮落带岩体压实承载的力学特性，能够可靠地反映出从松散低刚度状态到压实高刚度状态的力学特性演化，表征了垮落带岩体在覆岩压力作用下逐渐压紧密实的工程实际过程，定量表达了整体刚度随体积减小到最终压实的指数型增长规律，可表达垮落带岩体压缩、承载及支撑过程的力学行为。

1.3 垮落带模型力学特性参数的反演方法

在数值仿真模拟研究中，正确再现采空区垮落岩体压实-支撑过程的力学响应，是真实模拟工作面开采过程中支承压力分布、煤柱应力状态、巷道采动影响及其围岩稳定性的重要基础。针对采空区垮落带的数值模拟，国内外学者进行了大量探索，主要方法是采用弹性模型或双屈服本构模型(double-yield model)[8-9]充填采空区和垮落带，其中双屈服模型准确体现出垮落带岩体压实过程中的力学特性，近年来被科研人员所认可和运用。

将岩石力学试验所得到的参数代入式(3)～(7)，可得出该工作面顶板条件下采空区垮落带岩体碎胀后的最大应变εm、碎胀系数b和压实过程中的应力-应变关系。建立FLAC3D垮落碎胀岩体精细数值计算模型，采用双屈服本构模型进行压缩模拟试验，监测模型的应力-应变关系，与垮落带应力-应变关系的理论解对比，采用试错反演法可以得到拟合度较好的采空区模型力学参数。

2 裂隙带岩体的拉伸劣化力学特性

岩体抗拉强度很低，且远小于抗压强度及抗剪强度，在拉应力作用下极易使节理与裂隙扩展甚至贯通，工作面围岩特别是顶底板岩层中会产生拉应力及拉伸破裂。随着工作面推进及垮落带碎胀岩体充填采空区，上覆悬露岩层在弯曲拉应力作用下发生挠曲下沉及离层，在拉应力超过岩层抗拉强度的部位，拉伸裂隙扩展甚至产生拉破裂，破断运移后得到垮落带碎胀岩体的支撑，形成拉伸裂隙及离层裂隙发育的裂隙带。

国内外学者纷纷指出岩石完整性、裂隙产状和发育程度直接影响到岩体的力学性能，并分别提出了岩体力学参数(尤其是杨氏模量)与裂隙发育程度的关系[18-19]。由于采动裂隙带的复杂性和隐蔽性，现场实测其裂隙发育程度难度较大、成本较高，对裂隙带内裂隙发育程度的准确、定量描述非常困难，笔者基于国际普遍认可和广泛应用的地质强度指标GSI(Geological Strength Index)体系，研究得到岩体拉伸破坏后的残余杨氏模量与拉应力作用下的裂隙发育程度GSIt的关系[20]，即工程岩体拉伸劣化系数的表达式：

式中，Er为岩体发生拉伸破坏后的残余杨氏模量;Em为岩体初始杨氏模量;GSIt为岩体在拉应力作用下的裂隙发育程度。

由式(8)可见，随着岩体裂隙发育程度的增加，杨氏模量呈现快速下降的变化趋势，具有显著的模量劣化特性。裂隙带岩体的力学特性进入了峰后力学性质的劣化阶段，特别是裂隙发育导致的杨氏模量劣化。

工作面围岩破坏的主要形式是剪破裂和拉伸破裂，围岩峰后特性软化不仅是剪破裂状态下的应变软化，还有拉破裂状态下的杨氏模量劣化。基于Mohr-Coulomb强度准则的应变软化本构模型，通过定义力学参数与塑性变形的负相关关系，较好地反映出岩体峰后力学性能软化、并呈现残余强度的力学行为。而拉破坏分析仅判别出单元的破坏状态与破坏类型，尚未实现杨氏模量的劣化。由于FLAC3D中内置的Mohr-Coulomb模型和应变软化模型中杨氏模量均维持恒定，不能准确严谨的反映裂隙发育、纵横的裂隙带岩体力学特性。为此，基于连续介质有限差分数值模拟软件FLAC3D，采用FISH语言对内置的应变软化本构模型进行二次开发，增加拉伸裂隙模量劣化算法，形成工程岩体劣化模型[20]。

工程岩体劣化模型的算法核心由岩体破坏识别算法和破裂劣化算法组成，该模型在常规显式有限差分计算中，通过每间隔一定运算时步后遍历三维模型的全部单元体，逐步渐进地动态识别每个剪切破裂和拉伸破坏的单元，并对其进行相应力学参数劣化，其中黏聚力、内摩擦角弱化遵循应变软化模型，杨氏模量劣化遵循拉伸劣化模型，单元力学参数更新后继续迭代运算，循环上述过程直至模型达到平衡状态[20]。

工程岩体劣化模型在动态识别围岩破坏性态的基础上，根据破坏类型进行应变软化和杨氏模量劣化，通过力学特性等效的方法实现了岩体拉伸破裂后根据裂隙发育程度对杨氏模量的劣化，从而较真实严谨的反映裂隙岩体的力学特性，对覆岩裂隙带模拟与采动应力场演化研究具有重要意义。

3 采动应力场与覆岩两带耦合分析

3.1 采动应力场与覆岩两带耦合方法

基于上述采动覆岩“两带”破坏特征及其力学机制，考虑采空区垮落带压实承载力学特性与裂隙带岩体破断劣化特性的动态演化力学过程，研究提出采动应力场与采动覆岩的耦合分析方法，如图3所示。

图3 采动应力场与覆岩两带耦合分析方法Fig.3 Diagram of numerical simulation method on mining- induced stress field

根据工作面地质条件和岩石物理力学性质，计算或实测得到垮落带范围，依据垮落带岩体压实过程的力学特性，建立垮落岩体精细数值模型，采用试错反演得到符合垮落带压实过程应力-应变关系的模型力学参数，通过以上步骤实现采空区垮落带压实过程的力学特性模拟。计算或实测得到裂隙带范围，根据裂隙带岩体的拉伸劣化特性，实现上覆裂隙带中裂隙发育过程的模拟。在连续介质的数值模拟中，实现垮落带压实和裂隙带力学特性的仿真，使采动应力场分析与采空区覆岩力学特征相耦合，从而可以提高采动影响研究过程的严谨性和结果的可靠性。

3.2 采空区覆岩力学特性参数分析

为进行采动应力场与采空区覆岩“两带”力学特性的耦合分析，验证上述方法的可行性，以河南能化集团赵固二矿工程地质条件为例，运用上述方法分析确定采空区覆岩力学特性参数。

赵固二矿处于焦作煤田东部，11030工作面开采二1煤层，煤层倾角0°～11°，平均厚度6.16 m，煤层结构简单，层位稳定，属近水平稳定厚煤层。工作面采用大采高综合机械化开采工艺，后退式全部垮落法管理顶板，开采深度平均652 m。直接顶以泥岩、砂质泥岩为主，基本顶以粉砂岩、细粒砂岩为主，底板以砂质泥岩为主。二1煤层和顶板的岩石试样采自煤层顶板20 m和底板15 m范围，采用RMT-150型伺服试验机进行岩石力学试验，试验测试结果见表2。

表2煤岩物理力学参数
Table2Physicalandmechanicalpropertiesofrocksandcoal

岩石层位岩性视密度/(kg·m-3)抗拉强度/MPa抗压强度/MPa杨氏模量/GPa泊松比砂岩25521091836316023顶板砂质泥岩2591347562163025泥岩277723438295029二1煤煤层1435072204281030底板砂质泥岩262928141899027石灰岩27541352996826022

11030工作面采高6 m，根据岩石力学参数试验结果表2，由式(2)和表1，计算得到工作面采空区的垮落带高度为12.7 m，裂隙带高度为45.5 m。将表2中有关岩石力学参数代入计算式(5),(6)，得到二1煤层工作面采空区垮落带岩体碎胀后的最大应变εm为0.32，碎胀系数b为1.47。将εm及b两个参数代入计算式(3),(4)，得到垮落带岩体压实过程的应力-应变关系如图4所示。据此，采用FLAC3D数值计算模型对垮落带岩体进行压缩模拟试验，采用试错反演法得到与理论解拟合较好的模型应力-应变关系如图4所示，数值模拟反演得到的垮落带双屈服本构模型参数见表3。

图4 双屈服模型力学特性的数值模拟反演Fig.4 Back analysis of stress-strain characteristics of double-yield model

K/MPaG/MPaγ/(kg·m-3)c/MPaφ/(°)σt/MPa19900100017000001300

由岩石力学参数表2，采用Hoek-Brown强度准则计算得到岩体力学参数见表4。其中K为体积模量，G为剪切模量，c为黏聚力，φ为内摩擦角，cr为峰后残余内聚力，εp为岩体强度变为残余值时的塑性应变。数值计算模型的煤层和岩层采用表4的力学参数。工作面覆岩属于中硬覆岩条件，通过分析研究[21]，裂隙带岩体GSIt取值50。

表4模型岩体力学参数
Table4Rockmassmechanicalproperties

岩石层位岩性K/GPaG/GPac/MPaφ/(°)cr/MPaεp/%砂岩91593945039001顶板砂质泥岩52313240032001泥岩24112135021001二1煤煤13061431014001底板砂质泥岩72403437034001石灰岩96654247042001

4 应用实例

为验证上述方法的正确性和可靠性，以赵固二矿11030工作面工程为背景建立模型算例，应用上述模拟方法研究工作面采动应力场演化规律，并与理论分析结果进行对比印证。

4.1 数值模型的建立

根据对称性原则，以11030工作面中线为对称轴，建立三维数值模型如图5(a)所示。模型长度350 m，其中工作面推进长度250 m，前后各留50 m边界;模型宽度180 m，工作面二分之一的长度90 m，侧向煤柱宽度90 m;模型高120 m。模型顶部施加15 MPa的垂直应力，X,Y方向施加的水平应力分别为垂直应力的0.8和1.2倍[22]，模型四周和底部采用位移限定边界。随着工作面推进，垮落带冒落岩体由双屈服本构模型填充，力学参数见表3。各岩层采用工程岩体劣化模型作为力学判据。

图5 采动应力场演化FLAC3D模型Fig.5 FLAC3D for mining induced stress field simulation

在煤层直接顶布置如图5(b)中虚线所示的4条应力测线，其中测线1,2,3沿推进方向布置，分别位于工作面对称轴、距工作面边界45 m处和工作面边界;测线4垂直于推进方向布置，从采空区边界贯穿边界煤柱，工作面推进200 m后记录工作面前后和煤柱的垂直应力分布。

4.2 工作面前后支承压力演化规律

一次采高为3 m和6 m的工作面前后垂直应力分布如图6所示。可见，工作面开采促使围岩内部应力重新调整和分布，在工作面前方煤壁前方产生超前支承压力区，采空区垮落带岩体在上覆岩层作用下逐渐压实，出现卸压与应力恢复现象，随着远离工作面，采空区内垂直应力逐渐增大。

图6 工作面垂直应力分布平面图Fig.6 Plan view of vertical stress distribution induced by mining

工作面中部、距边界45 m及工作面边界位置(即测线1,2,3)的垂直应力分布如图7所示。可见，在工作面不同部位，前方煤壁和后方采空区的垂直应力分布存在着一定差异，沿着工作面长度范围(x=0～90 m)并非均匀分布。在整个工作面范围内，煤壁和采空区的垂直应力峰值点均位于工作面中部，由中部向边界逐渐降低。

图7 工作面前后支承压力分布Fig.7 Distribution of mining induced abutment pressure

由图6,7可见，采高增大对前支承压力分布有显著影响，主要体现为煤壁应力降低区扩大、支承压力峰值降低、峰值影响区域前移、影响范围扩大。采高增大对工作面中部前方煤壁影响最显著，采高由3 m增加到6 m时，前支承压力峰值由32.3 MPa降低至30.9 MPa，峰值位置前移8.4 m;在距边界45 m处，前支承压力峰值由30.6 MPa降低至29.4 MPa，峰值位置前移4.8 m;在工作面边界，前支承压力峰值仅由27.9 MPa降低至27.7 MPa，比中部峰值分别下降了4.4 MPa和3.2 MPa，峰值位置前移0.6 m。

由图6,7反映出工作面后方采空区冒落带压实造成的应力恢复。采空区应力恢复在工作面长度方向有着显著差异，采空区中部的应力恢复较高，距离采空区边界越近应力恢复越低。采高的增大也有明显影响，在采空区中部工作面后方110 m处，采高3 m时后支承压力恢复到原岩应力的88.7%，而采高6 m时仅恢复到79.3%，应力恢复速度较慢。这是由于垮落带和裂隙带的高度与采高成正比(见式(1),(2))，大采高开采扩大了覆岩破坏范围，使采动覆岩需要更长的时间才能经过应力重新调整达到稳定。

在距边界45 m的采空区部位，采高3和6 m时的采空区应力恢复差异不大，后支承压力升高速度基本相近，在工作面后方110 m处后支承压力分别恢复到原岩应力的73.7%和67.4%。而在采空区边界位置，应力恢复特征与采空区中部相反，采高为3 m时，后支承压力在工作面后方110 m处恢复到原岩应力的20.1%，而采高为6 m时恢复到30.0%，这一现象与其他学者的研究结果一致[8,23]。这是由于采空区边界位于侧向基本顶结构下方，围岩应力调整过程中应力主要向采空区中部压实区和煤柱转移，因此采空区边界应力恢复程度显著低于采空区中部。由于大采高工作面垮落带及裂隙带高度的增加，采空区中部大范围的覆岩破坏减缓了应力恢复速度，使围岩应力调整时部分应力由采空区中间向两侧转移;同时采高的增大使侧向基本顶结构高度增加，使采空区边界位置的覆岩压力随之增大。

4.3 侧向煤柱集中应力与力学状态演化规律

工作面开采除了在前方煤壁形成前支承压力，同时也在采空区侧向煤柱上形成垂直应力集中，煤柱边缘区域在集中应力作用下发生塑性破坏而卸载，集中应力逐渐向煤柱深部转移，塑性破坏区或应力极限平衡区也逐步向煤柱深部扩展直至弹性应力区边界[24-25]。

一次采高为3和6 m时侧向煤柱垂直应力分布如图8(a)所示，走向中部煤柱垂直应力分布如图8(b)所示。可见，采高对侧向煤柱的应力分布影响非常显著。采高为3 m时，侧向煤柱应力峰值位置进入煤柱内部12.9 m处，即极限平衡区宽度12.9 m，应力峰值34.1 MPa;采高为6 m时，侧向煤柱上应力集中区整体向煤柱深部移动，应力极限平衡区宽度24.2 m，且应力集中范围扩大，应力峰值为34.7 MPa。一次采出厚度由3 m增大至6 m时，侧向煤柱中应力极限平衡区宽度由12.9 m扩大到24.2 m，扩大系数为1.88。因此，采高大小对侧向煤柱应力极限平衡区宽度、集中应力程度有着显著影响，是沿空巷道合理煤柱宽度、稳定性分析与支护设计的重要因素。

图8 侧向煤柱垂直应力分布Fig.8 Distribution of vertical stress in pillar

煤柱应力极限平衡区宽度是设计煤柱尺寸、计算顶板侧向断裂位置和确定沿空掘巷合理位置等工作的重要指标。应力极限平衡区宽度可由式(9)得到[26]：

式中，x0为应力极限平衡区宽度;λ为侧压系数;φ0为煤层界面的内摩擦角;K为最大垂直应力集中系数;γ为上覆岩层体积力;H为埋深;c0为煤层界面黏聚力;Px为支架对煤帮的支护强度。

计算式(9)根据表3得到的参数取值，其中煤层与顶底板交界面的力学参数c0和φ0低于煤体的c和φ值，c0和φ0分别取0.7 MPa和26°。由此计算得到采高3 m时应力极限平衡区宽度x0为12.2 m，采高6 m时x0为24.4 m。上述数值模拟仿真结果与理论计算基本吻合。

上述工作面前后及侧向支承压力分布及演化规律，与采空区垮落带渐进压实承载力学特性、裂隙带覆岩运移特征、采空区及煤柱形态特征具有耦合关系。由此表明，采用研究得到的垮落带力学特性及其双屈服模型反演参数、裂隙带岩体劣化模型，实现了采动应力场与采空区覆岩特性的耦合，采动应力演化的仿真模拟是可靠的，该方法对区段煤柱稳定性、回采巷道采动影响与控制设计、岩层与地表移动、保护层开采等工程问题，具有重要的科学意义及应用价值。

5 结 论

(1)在综合分析覆岩“两带”高度计算及影响因素的基础上，研究得到垮落带渐进压实承载过程的力学特性表达式，垮落带岩体压缩刚度随压缩应变呈现指数型增长，碎胀系数和单轴抗压强度对压缩刚度及应变有显著影响，阐述了垮落带模型力学特性参数的反演方法，正确反映出采空区垮落带压实—支撑过程的力学响应。

(2)基于采空区垮落带压实承载过程的力学特性与裂隙带岩体的破断劣化特性，采用反演力学参数后的双屈服模型和二次开发得到的工程岩体劣化模型，等效描述垮落带和裂隙带的力学特性，提出了采动应力场与覆岩两带耦合方法，实现了采动应力场与采空区覆岩“两带”特征耦合分析。

(3)通过赵固二矿11030工作面应用实例分析，采用采动应力场与覆岩两带耦合方法，研究了工作面采动应力场演化规律，结果表明采高对超前支承压力分布、采空区应力恢复和侧向煤柱集中应力分布均有显著影响。将计算结果与理论分析进行对比，验证了采动应力场与覆岩两带耦合方法的可行性和正确性，提高了采动应力场和围岩采动响应研究的严谨性和可靠性，对区段煤柱稳定性、回采巷道采动影响与控制设计、岩层与地表移动、保护层开采等工程问题，具有重要的科学意义及应用价值。

[1] 谢和平,周宏伟,王金安,等.FLAC 在煤矿开采沉陷预测中的应用及对比分析[J].岩石力学与工程学报,1999,18(4):397-401. XIE Heping,ZHOU Hongwei,WANG Jin’an,et al.Application of FLAC to predict ground surface displacements due to coal extraction and its comparative analysis[J].Chinese Journal of Rock Mechanics and Engineering,1999,18(4):397-401

[2] BADR S A.Numerical analysis of coal yield pillars at deep longwall mines[D].Colorado:Colorado School of Mines,2004.

[3] 浦海，缪协兴.综放采场覆岩冒落与围岩支承压力动态分布规律的数值模拟[J].岩石力学与工程学报,2004,23(7):1122-1126. PU Hai,MIAO Xiexing.Numerical simulation of dynamic falling of overlying rocks and abutment pressure in surrounding rocks for fully-mechanized top-coal caving stope[J].Chinese Journal of Rock Mechanics and Engineering,1999,23(7):1122-1126.

[4] PAPPAS D M,MARK C.Behavior of simulated long-wall gob material[R].Washington:United States Department of the Interior Bureau of Mines,1993.

[5] 张振南,缪协兴,葛修润.松散岩块压实破碎规律的试验研究[J].岩石力学与工程学报,2005,24(3):451-455. ZHANG Zhennan,MIAO Xiexing,GE Xiurun.Testing study on compaction breakage of loose rock blocks[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(3):451-455.

[6] 张振南,茅献彪,郭广礼.松散岩块压实变形模量的试验研究[J].岩石力学与工程学报,2003,22(4):578-581. ZHANG Zhennan,MAO Xianbiao,GUO Guangli.Experimental study on deformational modulus of friable rock during compactionp[J].Chinese Journal of Rock Mechanics and Engineering,2003,22(4):578-581.

[7] 马占国.采空区破碎岩体压实和渗流特性研究[M].徐州:中国矿业大学出版社,2009.

[8] YAVUZ H.An estimation method for cover pressure re-establishment distance and pressure distribution in the goaf of longwall coal mines[J].International Journal of Rock Mechanics and Mining Sciences,2004,41(2):193-205.

[9] 白庆升,屠世浩,袁永，等.基于采空区压实理论的采动响应反演[J].中国矿业大学学报,2013,42(3):355-361,369. BAI Qingsheng,TU Shihao,YUAN Yong,et al.Back analysis of mining induced responses on the basis of goaf compaction theory[J].Journal of China University of Mining & Technology,2013,42(3):355-361,369.

[10] 钱鸣高,石平五.矿山压力与岩层控制[M].徐州:中国矿业大学出版社,2003.

[11] PALCHIK V.Influence of physical characteristics of weak rock mass on height of caved zone over abandoned subsurface coal mines[J].Environmental Geology,2002,42(1):92-101.

[12] BAI M,KENDORSKI F S,VAN ROOSENDAAL D J.Chinese and North American high-extraction underground coal mining strata behavior and water protection experience and guidelines[A].Proceedings of the 14th International Conference on Ground Control in Mining[C].West Virginia Univ.,Morgantown,WV (United States),1995.

[13] 金太,席京德,蒋金泉,等.缓倾斜厚煤层矿山压力与岩层控制[M].徐州:中国矿业大学出版社,2000.

[14] SALAMON M.Mechanism of caving in longwall coal mining[A].Rock mechanics contribution and challenges[C].Proceedings of the 31st US symposium of rock mechanics,Golden,Colorado,1990:161-168.

[15] WHITTLES D N,LOWNDES I S,KINGMAN S W,et al.Influence of geotechnical factors on gas flow experienced in a UK longwall coal mine panel[J].International Journal of Rock Mechanics and Mining Sciences,2006,43(3):369-387.

[16] JIANG Lishuai,ZHANG Peipeng,CHEN lianjun,et al.Numerical approach for goaf-side entry layout and yield pillar design in fractured ground conditions[J].Rock Mechanics and Rock Engineering,DOI:10.1007/s00603-017-1277-0.

[17] XUE Junhua,WANG Hanpeng,ZHOU Wei,et al.Experimental research on overlying strata movement and fracture evolution in pillarless stress-relief mining[J].International Journal of Coal Science & Technology,2015,2(1):38-45.

[18] HOEK E,CARRANZA-TORRES C,CORKUM B.Hoek-Brown failure criterion-2002 edition[J].Proceedings of NARMS-Tac,2002,1:267-273.

[19] MITRI H S,EDRISSI R,HENNING J G.Finite-element modeling of cable-bolted stopes in hard-rock underground mines[J].Transactions-Society For Mining Metallurgy and Exploration Incorporated,1995,298:1897-1902.

[20] JIANG Lishuai,SAINOKI Atsushi,MITRI Hani S,et al.Influence of fracture-induced weakening on coal mine gateroad stability[J].International Journal of Rock Mechanics and Mining Sciences,2016,88:307-317.

[21] 蒋力帅.工程岩体劣化与大采高沿空巷道围岩控制原理研究[D].北京:中国矿业大学(北京),2016.

[22] YAN S,BAI J,WANG X,et al.An innovative approach for gateroad layout in highly gassy longwall top coal caving[J].International Journal of Rock Mechanics and Mining Sciences,2013,59:33-41.

[23] LI W,BAI J,PENG S,et al.Numerical modeling for yield pillar design:a case study[J].Rock Mechanics and Rock Engineering,2015,48(1):305-318.

[24] HU Zhenqi,CHEN Chao,WU Xiao,et al.Surface movement and deformation characteristics due to high-intensive coal mining in the windy and sandy region[J].International Journal of Coal Science & Technology,2016,3(3):339-348.

[25] HU Haifeng,LIAN Xugang.Subsidence rules of underground coal mines for different soil layer thickness:Lu’an Coal Base as an example,China[J].International Journal of Coal Science & Technology,2015,2(3):178-185.

[26] 侯朝炯.巷道围岩控制[M].徐州：中国矿业大学出版社，2013：186-194.

Numericalsimulationoncouplingmethodbetweenmining-inducedstressandgoafcompression

JIANG Lishuai1,2，WU Quansen2，LI Xiaoyu2，DING Nan2

(1.StateKeyLaboratoryofMiningDisasterPreventionandControlCo-foundedbyShandongProvinceandtheMinistryofScienceandTechnology,ShandongUniversityofScienceandTechnology,Qingdao266590,China; 2.CollegeofMiningandSafetyEngineering,ShandongUniversityofScienceandTechnology,Qingdao266590,China)

The goaf stress recovery under the effect of overlying strata subsidence has notable effect on mining-induced stress field.Based on the coupling effect between ground stress and goaf compression,the mechanical characteristics and properties back-analysis method of caved rock are studied,the tension-weakening characteristics of rock mass in fractured zone are analyzed,and hereby the coupling simulation method between mining-induced stress and goaf compression is proposed and verified using a simulation case study.The results show that the mechanical behavior of the compressed caved zone follows an exponential growth,the compression stiffness nonlinearly increases with the compression strain.The mechanical properties of caved zone obtained from back-analysis agree with the analytical results,which contribute to the simulation of the mechanical behavior of goaf.The tension-weakening model,which is developed from strain-softening model by adding the tension-weakening algorithm,can describe the weakening of Young’s modulus,hereby achieve to describe the mechanical behavior of fractured rock mass rigorously.Based on the mechanical behavior simulation of caved and fractured zone,the proposed method can conduct the coupling numerical simulation of mining-induced stress field and the overlying “two-zone”.

mining-induced stress;caved zone compression;fractured zone weakening;numerical simulation;coupling analysis

10.13225/j.cnki.jccs.2016.1717

TD323;TD823

:A

:0253-9993(2017)08-1951-09

国家自然科学基金资助项目(51704182);山东省自然科学基金资助项目(ZR2017BEE050);山东省高等学校科技计划资助项目(J17KA212)

蒋力帅(1989—)，男，江苏南通人，讲师。E-mail:jlsh1989@126.com

蒋力帅，武泉森，李小裕，等.采动应力与采空区压实承载耦合分析方法研究[J].煤炭学报,2017,42(8):1951-1959.

JIANG Lishuai，WU Quansen，LI Xiaoyu,et al.Numerical simulation on coupling method between mining-induced stress and goaf compression[J].Journal of China Coal Society,2017,42(8):1951-1959.doi:10.13225/j.cnki.jccs.2016.1717