厚覆盖岩层下采矿引起地表变形的三维有限元分析

第22卷第10期岩石力学与工程学报22(10)：1629～16352003年10月ChineseJournalofRockMechanicsandEngineeringOct.，20032002年5月30日收到初稿，2002年7月4日收到修改稿。作者杨家岭简介：男，59岁，1968年毕业于武汉大学测绘学院，现任研究员，主要从事数值计算方面的研究工作。厚覆盖岩层下采矿引起地表变形的三维有限元分析杨家岭丰定祥葛修润(中国科学院武汉岩土力学研究所岩土力学重点实验室武汉430071)摘要介绍了在厚覆盖岩层条件下，确定地下采矿引起地表变形特征的三维有限元弹塑性计算分析工作。关键词岩石力学，厚覆盖岩层，地下采矿，地表变形，三维有限元法分类号P642.25，P642.26文献标识码A文章编号1000-6915(2003)10-1629-073DFEMANALYSISONGROUNDSURFACEDEFORMATIONCAUSEDBYMININGUNDERTHICKROCKOVERBURDENYangJialing，FengDingxiang，GeXiurun(KeyLaboratoryofRockandSoilMechanics，InstituteofRockandSoilMechanics，TheChineseAcademyofSciences，Wuhan430071China)AbstractTheanalysisof3Delasto-plasticFEMisintroducedforcalculatingthedeformationofthegroundsurfacecausedbyundergroundminingunderthickrockoverburden.Keywordsrockmechanics，thickrockoverburden，undergroundmining，settlementofgroundsurface，3DFEM1引言武汉钢铁公司金山店铁矿东区地下采场之上是居民稠密的张福山村和铁路线及车站，面积约0.9km2，地下矿体的开采所引起的地表变形极有可能危及这些建筑物的安全。为选择合适的开采方案，开展了地下开采引起地表变形特征及其控制对策问题的研究[1～5]。该采场有其自身的特点：在这里未进行过露天开采，不存在露天坑与地下采空区相连的问题；起采高程在地表以下200余米的深处，即地下采场上面有200余米的厚覆盖岩层。研究工作认为：厚覆盖岩层下的地下采空区可以看成是某种形态的地下空间，形成这样的地下空间时必然在围岩中形成破损区，具有这种破损区的地下空间视破损区的不同状态可能是稳定的，也可能是不稳定的；若地下空间是稳定的，则破损区常呈椭球状(在垂直于矿体走向的截面上呈椭圆形)，且破损区的顶部距地表有足够的距离。这样的地下空间及其破损区可以达到自身平衡稳定，因此不宜采用“渐进崩落”的连续破坏模式来影响地表的移动变形。为此，采用了数值计算分析和物理模型试验两种模拟方法，包括二维分析和三维分析，并且进行了多种可能的假定条件下的结果比较，在综合分析的基础上得到比较稳妥的地表变形结果。下面仅介绍三维有限元的计算结果及分析。2地质模型及力学参数金山店铁矿东区从25勘探线至44勘探线，东西长1000余米，南北宽900余米。采区范围内平行排列的3条断层F1，F3和F4近东西走向，南倾，倾角为87°～55°，其中，F4规模较大，横贯东西。因此，采区的地质结构大致可以分为3部分：以F4•1630•岩石力学与工程学报2003年断层为界面的上盘以沉积变质岩(如角页岩)为主；以F1断层为界面的下盘以石英闪长岩为主；F1～F4之间的断裂带，其中，包含较破碎的角页岩、岩浆岩、大理岩等，I，II号铁矿体即分布其中。I，II号矿体是该矿区的主要矿体，分布纵贯东西，两者南北相距15～200m，较为集中赋存的标高为－200～－500m，亦是拟将要开采的部分。采区地表地势西北高，最高高程海拨约90m，东南部高程降至30m。32#勘探线以西，分布着居民稠密的张福山村；38#和39#勘探线附近有铁路线和车站。各类岩石计算力学参数见表1。表1岩体计算力学参数Table1Mechanicalpamametersforcomputation岩性(或断层)容重γ/kN·m－3内摩擦角ϕ/(°)粘聚力c/MPa弹性模量E/GPa泊松比μ抗压强度/MPa抗拉强度/MPaA1：角页岩C1：岩浆岩27.0502.050.2575～800.4A2：角页岩(差)B：大理岩26.0451.050.30500.2C2：岩浆岩(差)27.0501.550.25650.3D：矿石34.0380.740.3550～400.1断层F1，F3，F4260.18法向刚度Kn/MN·m－3180.0切向刚度Ks/MN·m－360.00.05注：残余强度：岩石c残=0.5c，ϕ不变；断层c残=0.0，ϕ不变。3计算模型与计算方法3.1计算几何模型计算模型范围东西长1660m，南北宽1250m，垂直方向自地表至－800m，高约860m。三维有限元网格以8节点等参数单元模拟各类岩石和矿体，以8节点等参数节理单元模拟F1，F3，F4三条断层。全部网格包含15126个节点，11909个岩石单元(其中有703个矿体单元)和1536个断层节理单元(图1～3)。3.2初始地应力场(1)自重应力场图1三维有限元网格立体图Fig.1Meshfor3DFEM⎪⎭⎪⎬⎫−===zyxzhσμμσσγσ1(1)式中：γ为容重，h为埋深，μ为泊松比。(2)自重应力场+水平构造应力场⎪⎪⎭⎪⎪⎬⎫===hKhKhyyxxzγσγσγσ(2)式中：xK，yK分别为x，y方向水平应力的侧压力系数，取0.1=xK，5.1=yK。xσ，yσ作为边界的面分布荷载加在边界面上，与自重应力荷载共同形成计算体的初始应力平衡方程的荷载向量，解平衡方程得初始应力场。3.3开挖卸荷荷载在初始地应力场存在的条件下，模拟地下开采过程的荷载仅为开挖卸荷力。每个开挖单元的卸荷力由下式计算：vqNvBRvvd}{][d}{][}{0Te0Tee∫∫−=σ(3)式中：[B]为单元几何矩阵，}{0σ为单元初应力向量，[N]为单元形函数矩阵，}{0q为单元体积力向量。3.4三维弹塑性有限元计算方法弹塑性有限元按初应力法进行增量迭代混合求解。把每一开挖步的卸荷荷载分为若干级逐级施加，在每一级荷载增量中迭代求解，当前后两次迭代求得的位移之差小于预先规定的精度(10－4cm)后，再施加下一级荷载。岩石单元按理想塑性体对待，采用Drucker-Prager屈服准则：第22卷第10期杨家岭等.厚覆盖岩层下采矿引起地表变形的三维有限元分析•1631•图2地表网格及剖面示意图Fig.2Surfacemeshandsections图3－410m高程平剖面网格图Fig.3Planemeshatthelevelof－410m注：图中数字为单元号•1632•岩石力学与工程学报2003年KJJF−′+=21α(4)式中：ϕϕα2tan129/tan+=，c为粘聚力，ϕ为内摩擦角；ϕ2tan129/3+=cK；1J为第一应力张量不变量；2J′为第二偏应力张量不变量。当21JJ′+α＞K时，岩石单元被认为剪坏，屈服后单元的强度参数取残余强度，超出的应力转换为等效节点力转移到周围各节点上，继续迭代。当然，在剪屈服判别之前，先要进行拉破坏判别处理，被拉坏的单元不再进行抗剪处理。节理单元采用摩尔-库仑强度准则=F|τ|+−ϕαtanc，类似的处理方法见文[3]。4采矿方法的模拟及计算方案金山店铁矿东区地下开采拟采用无底柱分段崩落法。开采范围为－200～－550m中70m的一个阶段，共分为5个开采阶段，每个开采阶段分上40m和下30m两个开采步距，每个步距自西向东开采，先采I号矿体，后采II号矿体，I号矿体比II号矿体超前一个步距开采，在矿体垂直方向上开采顺序自上而下。实际计算中，在矿体走向方向上不再细分开采步，一次把一个步距的矿体全部采完。因此，整个开采模拟计算共分11个开采步。考虑到地表保护的重点对象是居民稠密的张福山村，设计在其相关部位的地下开采需保留必要的永久矿柱，包括部分连续矿柱和间隔矿柱(60m间隔)。即在32#勘探线以西、－340m高程以上设连续矿柱；32#勘探线以西、－340～－550m高程设间隔矿柱。计算中没有考虑模拟充填过程，开采过程的模拟只是分步把矿体单元去掉，形成一个采空区，并且视这种采空区为永久的地下空间。三维有限元计算将描述其采空效应，即地应力扰动释放调整的非线性过程，以及所产生的开采区围岩至地表的变形(位移)结果。矿体被开采后，采空区四周的围岩会出现一个松动区，这个松动区范围内的岩体对矿区顶板将失去或者部分地失去承载作用。在有限元模拟计算中，这个松动区的确定是依赖岩石单元的强度特性及当时的应力状况来圈定的，即当开采引起地应力二次分布后，若某岩石单元承受的应力强度超过了其按屈服准则计算的单元屈服强度，则认为该单元屈服或破损。将屈服单元所超载的应力转移到周围其他单元上，再进行平衡迭代，直到不再产生新的屈服单元。在整个矿体开采过程的非线性迭代计算过程中，采空区围岩的屈服松动由发生到扩展，有的情况在厚覆盖岩层下，只是局部的，且最后不再扩展，趋于稳定；有的情况屈服松动区不断扩展，直到地表，这种情况必会导致严重的地表陷落。三维有限元弹塑性计算做了若干个方案，主要方案的计算条件如表2所示。表2各计算方案条件表Table2Computationalschemesandconditions方案简称开采起止高程/m初始地应力场是否留永久矿柱屈服破损单元处理方式基本方案1～200～－550自重应力留保留部分承载能力基本方案2－270～－550自重应力留保留部分承载能力对比方案1－200～－550自重应力不留保留部分承载能力对比方案2－270～－550自重应力不留保留部分承载能力对比方案3－200～－550自重应力留全部丧失承载能力对比方案4－270～－550自重应力留全部丧失承载能力5计算结果分析上述各计算方案的主要计算结果对比见表3。基本方案1计算的地表最大沉陷区和地表最大水平位移区位置见图4。基本方案1计算的34#勘探线纵剖面、38#勘探线纵剖面以及－410m高程平剖面的采空区围岩屈服破损区见图5～7。计算结果表明：(1)基本方案1(－200m起采至－550m，张福山村部位下留永久矿柱，自重应力场)，最大地表沉降量达44.3cm，地表最大沉降区位于38#勘探线附近，靠近铁路线和火车站，最大沉降中心距张福山村约300m。该方案计算的最大地表倾斜1.1=Tmm/m，最大地表水平变形mm/m8.0=ε。根据《矿山测量试行规程》，对于一般砖石结构的建筑物，其临界值为：mm/m3=T，mm/m2=ε。计算结果对张福山村的安全性影响不大。第22卷第10期杨家岭等.厚覆盖岩层下采矿引起地表变形的三维有限元分析•1633•表3三维有限元各计算方案结果Table33DFEMcalculationresultsofdifferentschemes计算方案地表最大沉降/cm地表最大水平位移/cm最大地表倾斜T/mm·m－1最大地表水平变形ε/mm·m－1沉陷中心屈服破损区基本方案144.315.81.10.8明显偏东，38#勘探线附近主要分布在－100m高程以下，32#勘探线以东基本方案238.413.50.90.8明显偏东，38#勘探线附近主要分布在－100m高程以下，32#勘探线以东对比方案156.822.31.21.1西移160m有所扩大对比方案246.017.31.10.9西移160m有所扩大对比方案399.330.55.52.7扩展到地表对比方案483.224.63.42.7扩展到地表图4基本方案1地表最大沉陷区及地表最大水平位移位置图Fig.4Locationofthemaximumsubsidenceareamaximumhorizontaldisplacementofgro