基于FLAC3D的大采高综放工作面推进矿压显现规律数值模拟.doc

采矿工程数值分析与应用 题目：基于FLAC3D的大采高综放工作面推进矿压显现规律数值模拟 采矿工程数值分析与应用1 关键问题根据煤层、顶板冲击倾向性鉴定结果和曾发生的动力现象，并考虑到1301工作面复杂的开采条件(深部、特厚煤层、高地压、强承压水、高温、厚表土层、构造发育等)，认为1302N工作面开采面临潜在的冲击地压、帽裂等动力灾害威胁。因此，本文提出应用FLAC3D研究工作面推进过程中矿山压力显现规律。根据该规律，减少冲击地压等的动力灾害的发生，保证工作人员的生命安全，增大生产效率，提高产量。2 工程背景2.1矿井基本概况山东新巨龙能源有限责任公司位于山东省菏泽市巨野县新巨龙镇，在巨野煤田中南部，东距巨野县城约20公里，西距菏泽市40公里，兖新铁路和327国道在井田上穿过，北临日东高速公路，东依京福高速公路、京沪铁路、京杭大运河，西靠 京九铁路、济广高速、德商高速，交通便利。公司占地面积 522808 平方米，矿井井田东起田桥断层，西至煤采地层底界露头，南起邢庄及刘庄断层，北至陈庙 断层及第一勘探线，地理坐标为北纬35053530，东经1154711618，南北长约12公里，东西宽约15公里，面积约180平方公里。交通位置详见图2-1，其中A为新巨龙能源公司所在位置。矿井地质储量16.83亿吨，可采储量5.1亿吨，设计生产能力600万吨/年，设计服务年限82年。图2-1 交通位置图巨野属黄淮流域，北临黄河，境内水系健全，水资源丰富。既有充足的地表水、地下水，又可常年引流黄河水，全县水资源总量3.76亿m3，可利用地表水1.3亿m3，可利用地下水2.47亿m3，人均水资源储量413.1m3。即将动工兴建的大野水库，库容达2.5万m3，可为工农业发展和城镇居民生活提供用水保障。巨野属暖温带大陆气候，四季分明，气候温和，雨水充沛。年均气温13.5C；年均降水量655mm；无霜期平均213天；年日照时数2329.2-2578.3小时。2.2矿井地质条件2.2.1地层与主采煤层特点（1）煤层发育情况地层区划属华北地层区鲁西地层分区，区内多为第四系覆盖。3煤层平均厚度为8.82m，煤层直接顶为厚度为2.35m的粉砂岩，其单向抗压强度为32.13 MPa，基本顶为厚度为12.42m的粉砂岩互层，其单向抗压强度为100.17 MPa，由于煤层厚度较大，采空区冒落高度相对较高，煤层顶板又较为坚硬，因此可能存在顶板大面积悬顶。随着工作面的继续推进顶板集聚足够的弹性能，突然断裂对工作面支架、煤壁造成冲击诱发采场、巷道冲击地压发生。因此，3煤层顶板是否形成大面积悬顶是工作面顶板明显动压发生的必要条件。根据1301N工作面、已施工的1302N上下平巷揭露资料，工作面第一联络巷以南部西部区域3煤厚度在2.8-3.3m，东部靠近1301N工作面上平巷附近厚4.3-8.0m；第一联络巷以北西部区域3煤厚9.7-10.9m，1301N工作面上头揭露3煤厚8.0-9.5m；煤层中间夹0-1.2m泥岩或炭质泥岩。3煤层属较稳定，结构复杂。（2）煤层类型及煤质3层煤可采指数为1，为较稳定煤层，煤层倾角06，平均3；煤种在二联巷处分界，煤质牌号为肥煤；宏观煤岩类型为半暗-半亮型煤，低灰低硫易洗选，发热量30.13MJ/kg，硫份0.55%，灰分12.95%。（3）煤层顶底板煤层顶底板条件主要见表2-1：表2-1 煤层顶底板情况表（根据东部L-8钻孔）顶、底板名称岩石名称厚度特 征基本顶中砂岩28.5灰白色，夹深灰色条带，中厚层状，薄层状，水平层理，含长石，绿泥石、白云母片及黄铁矿，分选中等，硬度较大。直接顶伪顶粉砂岩1.3深灰色，薄层状，水平层理，成分较纯，含白云母片及较少植物茎叶化石，硬度较小，平坦断口，参差状断口。 直接底泥 岩2.52深灰色，厚层状，含黄铁矿，煤粒及植物要化石，具斜纵向裂隙，半充填黄铁矿片晶，硬度小，较脆，棱角状断口。老底细砂岩8.5灰白绿色，夹有粉砂岩薄层，内有方解石薄膜及黄铁矿晶体。2.2.2地质条件综合评价（1）采场矿压显现剧烈主采煤层上覆薄基岩厚度平均不到200 m，但基岩上部有约650 m的巨厚冲积层（如图4-1所示），加之采用大采高综放开采，采场矿压显现剧烈。矿井井巷工程埋深较深，上覆巨厚冲积层薄基岩，顶板岩石多为泥质砂岩，遇水即涨，性脆，易破碎，支护困难；煤层较厚，工作面上下平巷沿煤层底板掘进，巷道支护困难。加之采用大采高综放开采，采场矿压显现剧烈。（2）水压高龙固煤矿直接充水含水层是主采煤层（3煤）的顶底板砂岩（“三砂”）和底板的太原组“三灰”，水压都高达8 MPa。在大采高综放开采条件下，3煤顶板砂岩含水层将以淋水的形式进入矿井；由于采煤过程中产生的底板次生裂隙，3煤底板砂岩含水层将以涌水的形式进入矿井，共同构成矿井涌水量的一部分。（3）地温高矿井第一水平为-810 m，由于深部开采，采区内原始岩温高达46，出水温度最高达50，采掘工作面空气温度更是超过35，湿度达到100%。另外，除了超高地温（出水点温度50）、本矿井还存在超大涌水（1700 m3/h），根据现场实测，-810 水平水温达4851，-709水平水温达46。且同水平水温比岩温高46.3，水的放热量传递给风流的速度快，水温高是龙固矿井高温热害的重要特点。（4）易自然发火龙固煤矿主采煤层为3（3上、3下）煤层，煤尘爆炸指数为35.3240.43%，3层煤自燃倾向性等级为类，自然发火期46天。但由于龙固矿井井下岩层温度（平均38）较高，使煤的自然发火期缩短到41天。（5）工作面有动力现象1301工作面上、下平巷掘进过程中曾多次发生煤炮，工作面回采过程中，曾于2010年2月3日发生采场支架压死现象。根据煤层、顶板冲击倾向性鉴定结果和曾发生的动力现象，并考虑到1301工作面复杂的开采条件(深部、特厚煤层、高地压、强承压水、高温、厚表土层、构造发育等)，认为1301工作面回采过程中面临潜在的冲击地压等动力灾害威胁。2.3采煤方法采用大采高一次采全厚开采，割煤高度3.5-4.0m。工作面推进长度1500m，倾斜长度220m。回采巷道沿煤层底板掘进，巷道尺寸45004200mm。图2.2 煤岩层综合柱状图3 软件介绍3.1 FLAC/FLAC3D简介 FLAC（Fast Lagrangian Analysis of Continua）是由Itasca公司研发推出的连续介质力学分析软件，是该公司旗下最知名的软件系统之一。FLAC目前已在全球七十多个国家得到广泛应用，在国际土木工程（尤其是岩土工程）学术界和工业界享有盛誉。 FLAC有二维和三维计算软件两个版本，即FLAC2D（1984）和FLAC3D（1994）。这里进行一下说明，本书在阐述软件系列时，以FLAC统一称谓FLAC2D和FLAC3D；分述FLAC2D和FLAC3D时，FLAC仅指代FLAC2D。FLAC V3.0以前的版本为DOS版本，V2.5版本仅仅能够使用计算机的基本内存（64K），因而求解的最大结点数仅限于2000个以内。1995年，FLAC升级为V3.3的版本，由于能够扩展内存，因此大大增加了计算规模。FLAC目前已发展到V5.0版本。FLAC3D作为FLAC的扩展程序，不仅包括了FLAC的所有功能，并且在其基础上进行了进一步开发，使之能够模拟计算三维岩、土体及其它介质中工程结构的受力与变形形态。FLAC3D目前已发展到V3.1版本。 3.2 FLAC/FLAC3D的主要特点 FLAC/FLAC3D界面简洁明了，特点鲜明，其使用特征和计算特征在众多数值模拟软件中别具一格。3.2.1 FLAC/FLAC3D的使用特征 FLAC/FLAC3D的使用特征主要表现为： 命令驱动模式 FLAC/FLAC3D有两种输入模式： 人机交互模式，即从键盘输入各种命令控制软件的运行； 命令驱动模式，即写成命令流文件，由文件来控制软件的运行。其中，命令驱动模式为FLAC/FLAC3D的主要输入模式，尽管这种驱动方式对于简单问题的分析过于繁杂，对软件初学者而言相对较困难，但对于那些从事大型复杂工程问题分析而言，因涉及多次参数修改、命令调试，这种方式无疑是最有效、最经济的（当然，由于二维建模相对简单，照顾不少用户的使用习惯，在FLAC中也可以采用界面操作模式即GIIC模式进行分析计算）。 专一性 FLAC/FLAC3D专为岩土工程力学分析开发，内置丰富的弹、塑性材料本构模型（其中FLAC内置11个，FLAC3D内置12个），有静力、动力、蠕变、渗流、温度5种计算模式，各种模式间可以相互耦合，以模拟各种复杂的工程力学行为。 FLAC/FLAC3D可以模拟多种结构形式，如岩体、土体或其它材料实体：梁、锚元、桩、壳以及人工结构，如支护、衬砌、锚索、岩栓、土工织物、摩擦桩、板桩等。通过设置界面单元，可以模拟节理、断层或虚拟的物理边界等。 借助其强大的绘图功能，用户能绘制各种图形和表格。用户可以通过绘制计算时步函数关系曲线来分析、判断体系何时达到平衡与破坏状态，并在瞬态计算或动态计算中进行量化监控，从而通过图形直观地进行各种分析。 开放性 FLAC/FLAC3D几乎是一个全开放的系统，为用户提供了广阔的研究平台。通过其独特的命令驱动模式，用户几乎参与了从网格模型的建立、边界条件的设置、参数的调试到计算结果输出等的全部求解过程，自然能更深刻理解分析的实现过程。 利用其内置程序语言FISH，用户可以定义新的变量或函数，以适应特殊分析的需要。例如，利用FISH，用户可以设计自己的材料本构模型；用户可以在数值试验中进行伺服控制；可以指定特殊的边界条件；自动进行参数分析；可以获得计算过程中节点、单元的参数，如坐标、位移、速度、材料参数、应力、应变和不平衡力等。 此外，用户还可以利用C+ +程序语言自定义新的本构模型，编译成DLL（动态链接库），在需要时载入FLAC/FLAC3D，且运行速度与内置模型相差不大；用户也可以利用有限元软件或其它专业建模工具建立复杂三维模型，导入FLAC3D，以弥补在建立三维复杂模型等方面的不足。 3.2.2 FLAC/FLAC3D的计算特征 作为有限差分软件，相对于其它有限元软件，在算法上，FLAC/FLAC3D有以下几个优点： 采用“混合离散法”（Marti and Cundall 1982）来模拟材料的塑性破坏和塑性流动。这种方法比有限元法中通常采用的“离散集成法”更为准确、合理。 即使模拟静态系统，也采用动态运动方程进行求解，这使得FLAC/FLAC3D模拟物理上的不稳定过程不存在数值上的障碍。 采用显式差分法求解微分方程。对显式法来说，非线性本构关系与线性本构关系并无算法上的差别，根据已知应变增量，可很方便地求得应力增量、不平衡力并跟踪系统的演化过程。此外，由于显式法不形成刚度矩阵，每一时步计算所需内存很小，因而使用较少的内存就可以模拟大量的单元，特别适于在微机上操作。在大变形问题的求解过程中，由于每一时步变形很小，因此可采用小变形本构关系，将各时步的变形叠加，得到大变形。这就避免了推导并应用大变形本构关系时所遇到的麻烦，也使得它的求解过程与小变形问题一样。 3.2.3 FLAC/FLAC3D的求解流程 采用FLAC/FLAC3D进行数值模拟时，有三个基本部分必须指定，即：有限差分网格；本构关系和材料特性；边界和初始条件。 网格用来定义分析模型的几何形状；本构关系和与之对应的材料特性用来表征模型在外力作用下的力学响应特性；边界和初始条件则用来定义模型的初始状态（即边界条件发生变化或者受到扰动之前，模型所处的状态）。 在定义完这些条件之后，即可进行求解获得模型的初始状态；接着，执行开挖或变更其它模拟条件，进而求解获得模型对模拟条件变更后作出的响应。图 1-1给出的是FLAC/FLAC3D的一般求解流程。 对于多单元模型复杂问题，如动力分析、多场耦合分析等的模拟，可以按这一求解流程，先采用简单模型（单元数较少的模型）观察类似模拟条件下的响应，接着再进行复杂问题的模拟以使之更有效率。 3.3 FLAC/FLAC3D的应用范围 尽管最初开发FLAC是用于岩土工程和采矿工程的力学分析，但由于该软件具有很强的解决复杂力学问题的能力，因此，FLAC及其扩展软件FLAC3D的应用范围现已拓展到土木建筑、交通、水利、地质、核废料处理、石油及环境工程等领域，成为这些专业领域进行分析和设计不可或缺的工具。其研究范围主要集中在以下几个方面： 岩、土体的渐近破坏和崩塌现象的研究； 岩体中断层结构的影响和加固系统（如喷锚支护、喷射混凝土等）的模拟研究； 岩、土体材料固结过程的模拟研究； 岩、土体材料流变现象的研究； 高放射性废料的地下存储效果的研究分析； 岩、土体材料的变形局部化剪切带的演化模拟研究； 岩、土体的动力稳定性分析、土与结构的相互作用分析以及液化现象的研究等。 图 1-1 FLAC/FLAC3D的一般求解流程3.4 FLAC/FLAC3D的不足 毋庸置疑，FLAC/FLAC3D是十分优秀的岩土工程数值模拟软件，其实用性和专业性得到了广泛证实。但不可否认，FLAC/FLAC3D尤其是FLAC3D也存在着诸多不足，主要集中在以下几个方面： 求解时间受网格尺寸的影响很大。 对于一般的弹塑性问题，FLAC的求解时间大致与N3/2（N为单元数目）成正比，FLAC3D求解时间大致与N4/3成正比。由此可以看出，FLAC/FLAC3D对网格尺寸十分敏感，同一模型采用不同尺寸的网格单元可能导致求解时间相差数倍之巨。 某些模式下的计算求解时间很长。 由于很多物理过程（如固结过程、长期动力影响、材料流变等）与时间相关，模拟时必须考虑时间效应。对于这些物理过程的时间效应，FLAC/FLAC3D均采用真实时间予以考虑，因而造成求解时间很长，在有些情况下计算时间甚至是无法令人接受的。 前处理功能较弱。 FLAC3D对于复杂三维模型的建立仍然十分困难。尽管FLAC3D软件为用户提供了12种初始单元模型，通过连接、组合匹配这些初始单元模型可方便快捷地建立规则的三维工程地质体模型；同时，也可通过内置语言FISH，编写命令来调整、构建特殊的计算模型，使之更符合工程实际。但是，由于FLAC3D在建立计算模型时采用的是键入数据/命令行文件的方式，加上FISH语言独特的源代码表达方式，直接扼杀了一般工程技术人员运用FLAC3D进行工程分析的想法。即使对于有相当数值模拟经验和能力的分析人员来说，建立较复杂的地质体模型，如地形起伏大的峡谷区地质模型，也是一件费时费力的苦差。这也是造成FLAC3D三维模拟计算周期长、难度大，制约其进一步推广应用的主要原因之一（胡斌，张倬元，等2002；廖秋林，曾钱帮，等2005）。 尽管如此，FLAC/FLAC3D的不足之处还是可以采取一定办法予以克服的。其计算时步受网格尺寸影响较大和某些模式下计算时间过长的问题，由于涉及到软件内核即算法和计算效率的问题，需从算法和计算机性能上予以改进（戴荣，李仲奎等2006），普通用户是难以解决的。但随着算法的不断改进和完善，以及高性能计算机的普及，这些不足之处有望得到改善。至于FLAC3D前处理较弱的问题，普通用户即使在现有条件下也完全可以通过借助其它软件予以弥补和完善，对于这部分内容本书将单独成章予以详述。4 模型建立由于该地区煤层平均倾角为3度，比较小，可以近似的认为是水平地层，因此，本文将模型简化为水平地层的模型的开挖问题。以1302N工作面地质条件为基础，建立数值计算模型如图2所示。模型尺寸为400m*500m*77.24m，共划分为22720个单元和25704个节点，模型x轴方向为工作面倾向，y轴为工作面的推进方向（走向），z轴方向为铅锤向上，模型四周为固定边界条件，上边界为自由边界条件，施加16.6Mpa压应力载荷，沿模型x、y方向施加16.6Mpa的压应力载荷模拟水平应力。图2 模型5参数确定所涉及的岩体的具体参数见表5-1表5-1 所用岩体的力学参数岩层密度(kg/m3)体积模量（GPa）剪切模量（GPa）粘结力（MPa）抗拉强度（MPa）摩擦角（）细砂岩126609.971161244泥岩24905.544140 3煤14002.31.55.46342粉砂岩26003.93.021.740细砂岩2270010.07.56.01241粗砂岩25606.04.54.27306 计算过程计算过程如图6-1所示：图6-1 计算过程图7结果分析7.1 开挖后工作面的位移当工作面沿推进方向推进25米后，在20米处其垂向应力分布如图7-1所示，其位移云图上，在工作面中底板处位移值中最小，基本无变化，而工作面顶板中部变形最大，并向两端逐渐减小，呈现一个下凹的弧形，。图7-1沿推进方向20米处的剖面上垂直位移云图当工作面沿推进方向推进50米后，在45米处其垂向应力分布如图7-2所示，其位移云图上，在工作面中底板处位移值中最小，基本无变化，而工作面顶板中部变形最大，并向两端逐渐减小，呈现一个下凹的弧形，中间顶部出现最大位移,并且比图7-1的位移大。图7-2沿推进方向45米处的剖面上垂直位移云图当工作面沿推进方向推进75米后，在70米处其垂向应力分布如图7-3所示，其位移云图上，在工作面中底板处位移值中最小，基本无变化，而工作面顶板中部变形最大，并向两端逐渐减小，呈现一个下凹的弧形,中间顶部出现最大位移，较图7-2的位移更大。图7-3沿推进方向70米处的剖面上垂直位移云图当工作面沿推进方向推进100米后，在95米处其垂向应力分布如图7-4所示，其位移云图上，在工作面中底板处位移值中最小，基本无变化，而工作面顶板中部变形最大，并向两端逐渐减小，呈现一个下凹的弧形,中间顶部出现最大位移，较图7-3的位移更大。图7-4沿推进方向95米处的剖面上垂直位移云图7.2 开挖后工作面的垂直应力的分布当工作面沿推进方向推进25米后，在20米处其垂向应力