覆岩破坏特征数值模拟研究现状及进展

王晓蕾, 王 晶

(1.吕梁学院矿业工程系， 吕梁 033000； 2. 煤矿机械装备维护与检测试验吕梁市重点实验室， 吕梁 033000)

中国是一个煤矿大国，煤炭在中国能源结构中占有重要地位[1]，对于中国经济发展具有重要的促进作用，是中国经济不可或缺的能源，虽然近几年中国大力发展新能源，但是效果不明显[2-4]。《能源发展战略行动计划(2014—2020年)》指出，到2020年，一次能源消费总量控制在48亿t标准煤左右，煤炭消费总量控制在42亿t左右，其中煤炭占一次能源比重为62%。由此可见未来相当长的时间内中国能源还是以煤炭为主[5-8]。

随着中国煤炭资源的大量回采，优势地质条件下煤炭资源已趋于殆尽，煤矿企业不得不开采劣势地质条件下煤矿资源[9]，同时，随着开采深度的不断增加，煤层瓦斯含量较大，煤层透气性较低，严重威胁煤矿安全生产[10]。

煤层开采前其应力处于原始平衡状态，当煤炭资源回采后，受采动和支撑应力的作用[11]，煤层上覆岩层必然发生移动，最终破断，自上而下形成有规律的三带[12]，煤层破坏高度对于煤矿水害治理和瓦斯防治具有重要意义[13]。煤层发育特征受地质特征、开采条件等多种因素影响，煤层破坏高度的有效确定对于煤矿安全生产具有重要意义[14]。对于煤层破坏高度的研究目前主要采用理论分析、模拟实验、现场测试的方法，其中数值模拟实验操作简便，模拟不同工况的覆岩破坏，是一种非常好的测试技术[15]。

现基于煤层开采覆岩破坏特征数值模拟研究，总结煤层开采覆岩破坏特征数值模拟研究现状，分析煤层开采覆岩破坏特征数值模拟研究的不足，针对目前煤层开采覆岩破坏特征数值模拟研究存在的问题，提出未来煤层开采覆岩破坏特征数值模拟研究发展趋势。

1 煤层开采覆岩破坏

矿井中巷道开掘之前不受任何人类影响的岩体称为原岩体[16]。存在于原岩中的力为原岩应力，巷道采掘前岩体处于弹性变形状态[17]。假设上覆岩层为均质连续介质，岩体为半无限体，当距地表一定高度时，其主要受上覆岩层的压力作用[18-19]。

当巷道开挖后，原始应力被破坏，应力重新分布，巷道周围发生应力集中现象[20]。临空煤岩体应力升高，同时，其受力状态由三向变为两向受力状态，此时承重主要是采空区上部的煤岩体，其承受上部岩层应力作用[21]。并将这种作用传递到工作面两端的煤壁上以及工作面前后的煤岩体上[22]。在采空区区域内为煤岩体卸压区。而在工作面两端煤壁上形成了应力集中区[23]，如图1所示。

图1 工作面支撑压力分布Fig.1 Working face support pressure distribution

采煤工作面前、两端煤壁、采空区冒落的矸石分别承受的为前支承压力[24]、侧向的支承压力以及承载支撑作用的后支撑压力[25]。由于采空区的卸压以及应力集中区的存在，必然使得煤层顶板发生破坏[26]。最终自上而下形成有规律的三带，如图2所示。

Ⅰ为垮落带； Ⅱ为裂隙带； Ⅲ为弯曲下沉带图2 覆岩三带Fig.2 Three zones of overburden

2 覆岩破坏数值模拟方法

随着科学技术的不断发展，中外学者根据覆岩破坏的特征及机理提出了众多数值模拟软件，主要包括UDEC、FLAC、RFPA、CDEM、PFC、ABQUES。这些模拟软件对于煤矿安全生产起到了非常重要的作用。

2.1 UDEC

UDEC是一款基于离散单元法理论的一款计算分析程序[27]。它是处理不连续介质的二维离散元程序，可以模拟非连续介质承受静载或者动载作用下的响应[28]。

该方法求解主要包括物理方程和运动方程两部分，计算时采用迭代进行，直到每块物体不再出现不平衡力和力矩为止[29-30]。其物理方程为

F=-KU

(1)

式(1)中：F代表法向应力；K代表节理的法向刚度系数；U代表块体间的叠合量。

由于运动方程符合牛顿第二运动定律，采用动态松弛法进行求解，其方程为

(2)

式(2)中：u代表块体形心的位移；c代表黏性阻尼系数；k代表刚度系数；m代表块体单元质量；f代表外载荷；t代表时间。

采用动态松弛法[式(2)]进行求解，是一种显式求解。计算过程中对于时间差要求小，而且要求合理确定阻尼系数。

王国锋[31]采用UDEC模拟软件对赵庄煤矿覆岩三带发育特征进行了研究，以1307工作面作为模拟面，其模型长度和高度分别是400 m和60 m，其模型如图3所示。

图3 UDEC数值模拟模型Fig.3 Simulation model of UDEC

以工作面煤以及覆岩岩层物理力学性质作为模型参数，其力学参数如表1所示。

表1 UDEC模拟力学参数Table 1 Mechanical parameters of UDEC

模拟过程中选取左侧边界100 m进行开挖，得出开挖过程中覆岩变化如图4所示。

图4 UDEC模拟覆岩破坏特征Fig.4 Fault characteristics of overburden of FLAC

由图4可知，当工作面回采20 m后，覆岩产生了新边界，但未发生垮落，随着工作面的继续推进，当回采到30 m时，煤体上方发生了垮落，垮落高度可达13 m，当回采到80 m时，采空区受到上覆岩体的作用，采空区被压实，冒落带达到最大值36 m，当工作面回采到150 m时，裂隙带发育距顶板48.5 m。该位置布置瓦斯抽放巷对于煤矿安全生产具有重要意义。

徐刚等[32]采用UDEC对崔家沟煤矿覆岩裂隙演化进行了研究，分析了覆岩三带裂隙演化规律，得出了垮落带和断裂西发育高度；黄清锋等[33]采用该模拟软件对某矿6303大倾角工作面进行了模拟，得出了工作面覆岩破坏呈V形分布，采场后方破坏范围最大，中部次之；许刚刚[34]采用该方法对矿井才懂裂隙演化规律进行了研究，分析采动裂隙随开采宽度变化特征，对于采动裂隙定量化提供了实践基础；汤玉兵等[35]采用数值模拟对薛沙河下煤层覆岩裂缝发育特征和发育高度进行了研究，建立了该矿的数值模拟模型，得出了导水裂缝带发育高度。

2.2 FLAC

FLAC是1985年Cundall博士提出来的一种工程力学计算的有限差分程序，主要用于分析岩土工程[36-37]，该软件具有多种结构单元，能够有效分析岩土及其他材料组成的复合结构物，其内置多种材料本构模型[38]，能够研究复杂工程的力学行为，得出各种云图及曲线图，具有很好的模拟效果[39]。

流固耦合机制是该模拟模型中比较复杂的内容，其主要采用单相达西流进行处理，其主要包括运移定律、平衡定律、本构法则、流体边界条件4部分[40]。

对于流体的运移主要采用达西定律进行描述，即

(3)

流体的质量遵循质量平衡方程，即

∂ξ/∂t=-∂qi/∂xy+q0

(4)

式(4)中：ξ为流体质量变化量；q0为初始渗流速度；xy为模型y轴方向长度。

孔隙流体的反应方程取决于饱和度的值，主要对完全饱和状态下的本构法则进行分析。

∂P/∂t=M(∂ξ/∂t-α∂ε/∂t)

(5)

式(5)中：M代表比奥模量；ε代表体积应变；α代表体积系数。

而流体的边界初始条件，孔隙压力为0，饱和度为1，因此透水边界条件为

qn=h(p1-p2)

(6)

式(6)中：qn代表渗流速度分量；h代表岩土透水系数；p1代表边界孔隙压力；p2代表透水区域孔隙压力。

张开弦[41]采用FLAC对祁东煤矿开采覆岩变形破坏规律进行了研究。该煤矿生产能力2.4 Mt/a，存在较多含水层，最大厚度可达59.1 m，是最大额突水层，是矿井开采主要威胁水源。

试验采用三维模拟模型，走向、倾向、高度分别是600、600、350 m。并在顶部施加一定载荷模拟覆岩自重。其力学参数如表2所示。

表2 FLAC数值模拟参数Table 2 Numerical simulation parameters of FLAC

根据力学参数和岩层厚度得出其计算模型如图5所示。

图5 FLAC数值模拟模型Fig.5 Simulation model of FLAC

共开挖400 m，分40步并在模型两边设置保护煤柱。得出开挖过程中覆岩位移特征如图6所示。由图6可知，当工作面开挖后，覆岩发生破坏，经历了垮落、开裂、离层等几个阶段。煤层开挖初始阶段弯曲程度较小，随着工作面的继续开挖，下沉继续，随着工作面距离的增加，覆岩下沉量越小，随着工作天的继续开挖，周期性进行垮落，当推进到200 m时，覆岩呈拱形，随着继续推移，当工作面完全回采结束，覆岩云图趋于完整。

图6 FLAC模拟覆岩破坏特征Fig.6 Fault characteristics of overburden of FLAC

冯超等[42]采用FLAC对小保当煤田煤层开采引起的破坏进行了研究，分析了导水裂缝带发育高度；杨高峰等[43]采用该软件对镇城底煤矿8号煤层顶板两带发育高度进行了研究，得出了覆岩导水裂缝带发育高度，为矿井防治水提供了技术支撑；郝利生等[44]采用该软件对斜沟煤矿18205工作面开采覆岩两带破坏特征进行了研究，得出了两带发育高度，并通过实测验证了模拟的准确性；孙润等[45]采用该软件对屯兰煤矿工作面覆岩导水裂缝带发育高度进行了研究，确定了最大破坏高度。

2.3 RFPA

RFPA软件主要应用弹性力学来分析煤岩体破坏应力特征[46-47]，在分析的过程中考虑了岩石本身的非均值性以及随机性缺陷[48]。同时，模型中的单个单元对整体影响有限，当模型中足够多基元数量力学行为的叠加对模型整体产生影响[49]，就利用大量基元力学行为的集体效应实现对岩石破坏失稳过程的全程模拟[50]。

由于基元赋值符合Weibull分布，因此建立细观与宏观介质力学联系，其计算公式为

φ(α)=ϑ/α0[(α/α0)ϑ-1e-(α/α0)ϑ]

(7)

式(7)中：α为煤岩体力学性质；ϑ为分布函数的性质参数；α0为力学参数平均值；φ(α)为煤岩体介质的力学统计分布密度。

假设模型的所有弹模平均值为E0，其分布函数积分为

(8)

式(8)中：φ(E)代表分布值；E代表弹性模量。由式(8)统计分布组成一个样本空间，由于性质参数的不同，得出的空间分布也不同。

段俭君[51]采用RFPA软件对任家庄煤矿工作面开采覆岩破坏高度进行了研究。该煤矿主要开采3、5号煤层，工作面顶板以泥岩、粉砂岩为主，采用走向长壁一次采全高采煤方法。建立400 m×200 m的模型，其力学参数如表3所示。

表3 RFPA模拟力学参数Table 3 Mechanical parameters of RFPA

根据煤岩体力学性质以及厚度进行建模，得出其模型如图7所示。

图7 RFPA数值模拟模型Fig.7 Simulation model of RFPA

模拟煤层从左向右依次开挖，开挖长度设置为120 m。共开挖12步，每步10 m。得出工作面开采过程中覆岩破坏特征如图8所示。

由图8可知，当工作面推进到30 m时，直接顶发生了垮落，随着工作面的继续推进，当推进到50 m时，裂隙发育高度继续向上部延展，并且老顶初次来压，当工作面推进到80 m时，顶板周期性来压，裂隙发育，裂隙带发育高度达到35 m处。

图8 RFPA模拟覆岩破坏特征Fig.8 Fault characteristics of overburden of RFPA

赵韶波[52]采用RFPA对寺河二号井覆岩破坏高度进行了研究，确定了覆岩三带发育高度，并用高位钻孔对覆岩破坏高度进行了验证；崔永青等[53]采用该软件对马堡煤矿工作面开采覆岩三带进行了研究，得出垮落带和裂隙带破坏高度，与实测结果一直，验证了高技术的有效性；孙显龙等[54]采用该软件对煤矿顶板破坏特征进行了研究，得出了工作面回采过程中覆岩破坏特征；颜林艳等[55]采用该技术对贵州煤矿开采覆岩破坏特征进行了分析，开采承重覆岩重力形成应力拱，构造的存在是覆岩破坏的途径。

2.4 PFC

PFC是颗粒流程序，它是由美国Itasca公司开发的一种计算模拟软件[56]，主要用于研究散粒体的分析，属于离散元软件范畴[57]。其理论基础主要为牛顿第二定律和力与位移的关系，通过微观角度来分析材料的宏观行为[58]。当颗粒所受到的不平衡力达到预设值时，模型则会停止计算，此时模拟结束[59]。

该颗粒流软件主要包括物理方程和运动方程两个方面[60]。

对于物理方程，主要分析颗粒之间的距离，其计算公式为

(9)

式(9)中：xA、xB代表两颗粒的圆心位置。

两颗粒之间的重叠量计算公式为

U=RA+RB-d

(10)

式(10)中：U代表颗粒之间重叠量；RA、RB代表两颗粒的半径。

对于运动方程，其加速度和角加速度计算公式为

ux(t0)=Fx/m

(11)

ωx(t0)=Mx/I

(12)

式中：ux代表计算度；Fx代表合力x方向的分量；m代表颗粒的质量；I代表转动惯量；ωx代表角加速度；Mx代表颗粒受到力矩在x方向的分量；t0代表某一时刻。

陆敬锋[61]采用PFC软件对新疆大南湖矿区采动覆岩运移规律进行了研究。建立模型长、高分别是180、100 m。模型顶板施加2.5 MPa的应力，其力学参数如表4所示。

表4 PFC模拟力学参数Table 4 Mechanical parameters of PFC

根据煤岩体的厚度以及力学参数建立数值模型如图9所示。

图9 PFC数值模拟模型Fig.9 Simulation model of PFC

将岩层分为20个岩层，煤层设置高度为3 m，每次开挖10 m，共开挖10次，开采长度为100 m。开挖过程中覆岩破坏特征如图10所示。

由图10可知，当工作面推进到10 m时，工作面顶板发生下沉现象，直接顶出现了裂缝，直接顶没有垮落；当工作面推进到50 m时，由于关键层的存在，能够承受较大的载荷，抑制顶板下沉，垮落带高度为6.56 m；当工作面推进到100 m时，裂隙向深部发展，含水层贯通，导水裂缝带发育到下部含水层，上部隔水层没有明显破坏，未出现贯通裂隙。

图10 PFC模拟覆岩破坏特征Fig.10 Fault characteristics of overburden of PFC

江成浩等[62]采用PFC软件对滕东煤矿综放工作面裂隙演化规律进行了研究，得出了覆岩孔隙率；钟江城等[63]采用该软件对浅埋煤层覆岩裂隙演化进行了研究，裂隙开度与倾角是影响矿井突水突砂的主要影响因素；罗怀廷等[64]采用该软件对固体充填采煤中的覆岩破坏特征进行了研究，破碎矸石在侧限压缩条件下所产生的侧压力与其所受轴压成正比，粒径越小颗粒转动越活跃，大颗粒自身具备运动惰性。

2.5 CDEM

CDEM软件是由中国科学院和北京极道成然科技有限公司联合开发的力学分析系列软件，该软件是基于连续介质的力学离散方法[65-66]，主要用于模拟材料的变形和破坏过程，软件单元之间通过虚拟弹簧传递相互作用力，在计算的过程中块体单元可以实现断裂[67]。模拟材料的破坏时通过界面的张开、滑移和运动速度等表述，非常适合描述地质体由连续向非连续的破坏过程[68]。

该软件采用刚度矩阵法进行求解，在每个迭代步通过式(13)计算单元自身的节点力，并将此节点力分配至单元对应的节点上。

Fi=KiuiC

(13)

式(13)中：Fi代表单元i的节点力向量；Ki代表单元i的节点位移向量；ui代表单元i的单元刚度矩阵；C代表黏聚力。

其界面的弹簧力计算公式为

(14)

式(14)中：Fn、Fs分别代表弹簧的法向、切向应力；Kn、Ks分别代表弹簧的法向、切向刚度；Δdn、Δds分别代表弹簧的法向、切向位移。

当破坏计算时，采用莫尔库伦准则进行修正，修正的计算式为

(15)

式(15)中：T代表抗拉强度；φ代表内摩擦角。

焦振华[69]采用CDEM软件对寺河二号井保护层开采覆岩破坏特征进行了研究。试验中采用模型长和高分别是200 m×72.2 m，其力学参数如表5所示。

表5 CDEM模拟力学参数Table 5 Mechanical parameters of CDEM

根据煤岩体厚度及力学性质建立相关模型如图11所示。

图11 CDEM模拟模型Fig.11 Simulation model of CDEM

为了消除边界效应，在模型的两侧留设50 m的煤柱，并在模型顶部施加7.5 MPa的岩层自重力，得出工作面回采中覆岩破坏特征如图12所示。

图12 CDEM数值模拟覆岩破坏特征Fig.12 Fault characteristics of overburden of CDEM

由图12可知，当工作面回采到20 m时，工作面顶板发生初次垮落，随着工作面的继续推进，直接顶随着采动而发生冒落，当工作面推进到30 m时，基本顶发生垮落。最终工作面导水裂缝带发育高度为30 m。

王振伟等[70]采用该方法对安家岭煤矿覆岩破坏进行了研究，得出了弯曲带与垮落带和裂隙带的厚度比为1∶2，初次垮落步距为115～165 m；郝应祥[71]采用该软件对孙家沟煤矿覆岩开采对地表的影响进行了研究，得出了顶板破坏规律和位移的变化特征，确定了地表沉降范围。

2.6 ABQUES

ABQUES是基于有限元的软件，其主要研究的是相对简单的线性以及众多复杂的非线性问题。可以很好地模拟各种形状的单元库以及众多类型材料的性能[72-73]。其适合冲击载荷和爆炸、地震这种短暂的瞬时的动态事件的分析与研究，对于高度非线性问题也非常有效，包括非线性大变形等问题的求解研究[74]。

其工作面开采覆岩变化的黏聚力模型服从Drucker-Prager屈服条件，即

(16)

式(16)中：I1、J2分别代表应力张量第一和第二不变量；α、k代表关于内摩擦角和黏聚力的参数。

砌体界面单元的本构关系为

t=Eε

(17)

式(17)中：t代表应力向量；E代表弹性模量；ε代表名义应变向量。

黏结单元指数计算公式为

(18)

式(18)中：D代表损伤变量；δ0代表初始的有效位移量；δmax代表加载过程中的最大位移量；δf代表失效时的有效位移量。

田森等[75-76]采用ABQUES对保护层开采覆岩裂隙演化规律进行了研究。该煤矿主采1号煤层，煤层埋深可达1 100 m，开采区域长度为90 m，其力学参数如表6所示。

表6 ABQUES数值模拟力学参数Table 6 Mechanical parameters of ABQUES

根据煤岩体的力学性质和厚度建立煤层采动有限元模型如图13所示。

图13 ABQUES数值模拟模型Fig.13 Simulation model of ABQUES

开采过程中，通过杀死单元的方式进行煤层的开挖，煤层开采共分9步进行，每步开挖90 m，共开挖90 m，覆岩变化特征如图14所示。

图14 ABQUES模拟覆岩破坏特征Fig.14 Fault characteristics of overburden of ABQUES

由图14可知，当工作面推进到30 m时，覆岩应力发生了明显变化，且靠近采空区中心应力变化最为显著，也是卸压最快的区域；当工作面回采到60 m时，覆岩应力继续发生着变化，同样，采空区附近应力变化最为明显，形成了驼峰状变化曲线；当工作面推进到90 m时，覆岩应力重新分布，卸压区范围逐渐变大，采空区附近是最大，在支撑应力区域内卸压效果较差。

王迁等[77]采用ABQUES软件对司马矿1112工作面煤层开采覆岩变化特征进行了研究，得出了下沉特征以及不同推进距离下与下沉量的关系和工作面覆岩应力变化特征。

3 存在问题及发展趋势

3.1 存在问题

(1)对于工作面回采覆岩破坏特征的数值模拟研究，多以平面模型为主。而采矿工程是一个三维模型，其工作面的开采覆岩变化特征是一个时空演变的复杂过程。单纯的平面模型难以准确地模拟覆岩破坏特征，难以反映开采覆岩破坏的实质问题。

(2)与简单的岩土工程对比，工作面回采过程中覆岩移动、损伤、破坏、失稳是必然要发生的。目前的模拟软件只能对理想的弹性、塑性以及损伤体进行有效的变形和受力进行分析。而采矿工程为更重要的是材料和结构破坏后的力学特性、煤岩体是如何破坏以及破坏稳定的变化过程。

(3)数值模拟软件在模拟过程中会对边界条件以及材料属性进行简化，多对模拟结果产生一定影响，并且模拟结构离散形式不同，模拟的结果和精度也不相同，随机性较大，造成模拟结果的可信度大大降低[78-80]。

3.2 发展趋势

(1)随着中国煤炭开采技术的不断发展以及环保和可持续发展理念的加强，提出了绿色开采和智能化开采，这些技术属于较前沿课题，缺乏一定的技术支持。未来数值模拟应该结合绿色开采和智能化开采技术研发新型的模拟模块。适合煤矿技术的发展。

(2)采矿覆岩破坏是一个时空变化过程，平面模型不足以模拟开采过程。未来应研发可视化的三维模拟模型，高强度还原采矿过程中覆岩破坏特征，提高模拟准确度。

(3)数值模拟仅仅考虑开采深度和采厚因素，对于其他地质、水文因素几乎无涉及，造成模拟准确度不高。未来数值模拟应综合覆岩破坏多影响因素，优化煤岩体力学参数合理确定方法，引入新的思维方法和数学工具，提高模拟结果的准确度。

(4)随着煤矿开采深度的不断增加，对于软岩在高地应力蠕变特性突出。在现有数值模拟软件中未考虑该问题。未来数值模拟应考虑高地应力下软岩蠕变特性，对于保证采矿工程的长期稳定性至关重要。