5 矿岩自然崩落模拟技术研究

1、 五、矿岩自然崩落模拟技术研究5.1矿岩崩落研究的技术路线自然崩落法是一种利用地压进行采矿的方法，在拉底与削帮工程完成后，矿岩中会产生相应的次生应力场，岩体在次生应力场的作用下发生破坏直至崩落，这是人们普遍认可的矿岩自然崩落的机理。但是，对于不同的岩体，在次生应力场的作用下，其破坏方式一般都不相同。对于坚硬脆性岩体，崩落可能是由于岩体的脆性断裂；而对于软弱复杂岩体，矿岩的破坏可能会以大变形的方式导致崩落产生。我们知道，岩体的破坏与时间关系最为密切的因素是岩石的流变性质和岩石中微裂隙的亚临界扩展速率。岩石的流变性主要表现在蠕变和松弛两方面。所谓蠕变是指在应力不变的情况下，变形随时间的持续而逐渐增

2、长的现象。松弛是指在变形保持不变的情况下，应力随时间的持续而逐渐减小的现象。特别是蠕变现象在一些蠕变效应比较明显的岩石中更为常见。在这些岩石中，岩体工程的毁坏，往往不是因为围岩的强度不够，而是由于蠕变而产生过大的变形，致使岩体工程发生毁坏。因此，对这类岩石中的岩体工程进行设计时，必须考虑蠕变的影响。故很有必要研究岩石的蠕变特性。另一方面，岩石作为一种脆性材料，其中含育大量的微观裂隙，当裂隙尖端的强度因子接近岩石的断裂韧性时，裂纹发生亚临界扩展，经过一定的时间后，裂纹就会发生高速扩展从而导致断裂破坏。研究表明，裂纹从发生亚临界扩展开始到高速扩展所需的时间是一个与作用在裂纹上的应力有关的量。对于脆

3、性岩体，崩落速率的研究就可以从考虑裂隙的扩展与时间的关系入手，但对于金川软弱岩体，裂隙脆性断裂的可能性不大，与破坏有关的时间因素则主要是由岩体的流变性所决定，因此，我们可以通过取样进行流变试验和断裂力学试验来确定金川三矿区岩石力学基本参数。在此基础上，确定从流变的角度来预测崩落速率或是采用岩石的亚临界扩展速率作为判断条件来预测崩落速率。研究自然崩落法的崩落规律一般通过这样两种方法进行，一是采用做相似材料模拟试验的方法来直观地认识崩落的规律，二是采用数值分析的方法来模拟自然崩落的规律。由于技术上的原因，在对自然崩落法的可崩性进行三维分析时，相似材料模拟难以实现多因素的相似，甚至一些主要因素要完全

4、满足相似条件都会存在困难。数值分析法模拟崩落规律有其方便的地方，那就是可以进行多因素、多模型的分析和比较。在这方面，我们通过自行编制三维粘弹塑性有限元程序和采用FLAC3D相结合的办法，对金川三矿区矿体崩落特性进行研究。研究的目的是对各种条件下不同拉底方案以及同一拉底面积下不同的割帮高度和不同的原岩应力场进行模拟计算与分析，通过数值分析的方法将这些规律性的东西进行全面的反映，得出在这些条件下岩体应力的变化和岩体的破坏规律。这一技术路线用框图表示如图5-1所示。5.2模拟程序GEOVEP3D的实现概述随着微型计算机的飞速发展和普及，计算机在各学科中应用得越来越广，岩土学科亦不例外。计算机硬件发展

5、很快，不但性能越来越好，而且价格越来越便宜，加之系统软件的发展，如PowerStation4.0Fortran90编译系统，计算速度成千上万倍地加快，用粘弹塑性有限元法分析岩土实际问题，单元数和计算时间几乎不受限制，因此，粘弹塑性有限元法在岩土工程中的实际应用将越来越普及。岩石流变试验研究岩石断裂与扩展速率试验国内外研究现状调研图5-1矿岩崩落研究技术路线图有限元方法是分析岩土工程问题最有效的方法，经过几十年的发展和应用，已经比较成熟,且积累了大量的程序资料。但大多数有限元程序都是针对结构问题开发出来的，用于岩土工程问题的分析还存在着许多困难，尤其重要的是没有引入很适合于岩土材料的粘弹塑性模型

6、；其次有限元前后处理程序不完善，有限元程序的数据输入量大，很难排除数据输入错误，这些问题使得大家很难使用有限元程序分析岩土工程问题。因此，有必要对岩土粘弹塑性有限元程序的前后处理及岩土材料模型的引入进行研究。国外的有限元软件发展较快，功能不断地增强，在我国的应用领域也在不断地扩大，如SAP95在实际应用中确实发挥了较大的作用。但国内的有限元软件发展较慢，还没有一个有限元软件能真正大量地推广使用，大多是研究人员根据自己特定的需要编制出相应的有限元分析软件，并未达到商品化的程度。长此下去，使用国外有限元软件必然会带来一系列的问题，如运行结果的可靠性难以保证，很容易受国外公司人为地控制等。另一方面国

7、外还没有适合岩土工程分析的粘弹塑性有限元程序。因此，有必要开发适合中国国情的岩土（粘）弹塑性有限元程序国外比较成熟的BOPACE3D程序，主要用于结构的热弹塑蠕变分析，是一个功能比较完善的三维非线性等参单元有限元程序，该程序既考虑了材料非线性，又考虑了几何非线性。BOPACE原来用在IBM360/370和UNIVAC1108机工作站上，且有较完整的源程序。完全自己开发岩土粘弹塑性有限元程序必然会耗费大量的精力和财力，而利用BOPACE源程序进行改进可节约大量的时间，也是可行的。因此，本文作者经过一年多的努力，对BOPACE源程序的原理和结构进行消化和理解，在该程序的基础上开发出了适合于岩土工程

8、问题的粘弹塑性有限元程序GEOVEP3D,包括前后处理程序。BOPACE的源程序是用FORTRANIV语言编写的，将其转化成微机版本，其前后处理程序无法移植过来，必须重新编制。BOPACE程序的研究对象主要是金属材料，须将其转化成适合岩土材料的程序，且在转化过程中存在着许多语句错误，因此，必须读懂源程序，才能对错误进行必要的修改。修改源程序并不是一件很容易的事，但比重新开发岩土粘弹塑性有限元程序要容易得多。程序的求解方法及功能近年来，人们已经日益认识到在诸如土建、航空和机械工程中所遇到的各种类型的结构研究和设计中，使计算机具有有效的线性和非线性分析能力是一项很重要的和有价值的工作。首先，复杂结

9、构的线性和非线性分析完全可以利用微型计算机在实际结构的离散模型上进行，人们已经证明有限元法是一种非常有效的离散方法，这种方法既能用于分析结构和固体力学中的问题，又能用于分析热传导(及其他场)、液体流动、电磁等领域中的问题。本文将在BOPACE3D程序的基础上进行，对其消化、修改、补充后，研制适合于岩土工程及结构工程的粘弹塑性有限元程序。由于在许多情况下，只要求作线性分析，而且非线性分析总是先作线性分析，因此，本文设计了一个能非常有效地执行线性分析的程序，而且在线性分析之后，只要作相当少的输入改动，就能进行非线性分析。单元类型及材料种类程序所包含的单元类型有八节点和二十节点等参单元，非线性性质可

10、以由大位移、大应变和材料非线性性质引起。目前，可用的材料种类有：线性弹性、非线性弹性、曲线模式、混凝土模式、粘弹塑性等。求解特点该有限元软件是一个芯外存贮的解题程序，其平衡方程按块处理，它可以求解很大的有限元系统。另外，所有的结构矩阵按紧凑格式存贮，即只有非零元素才被处理，所以能得到最大的解题规模和求解效率。除了平衡方程的芯外存贮解法之外，实际上并不存在高速存贮器对所能使用的单元总数限制。为了获得程序最大的解题能力，根据其单元类型以及它们是线性单元还是非线性单元等不同情况按块处理。在求解过程中，低速存贮器用来贮存每块单元的全部信息。在非线性单元情况下，这些信息在时间积分过程中被修正，同时所需的

11、低速存贮容量，控制着所能处理的有限元系统规模。在非线性分析中，用平衡方程增量解法来计算有限元系统的响应，所采用的增量求解方案是一种加速的修正牛顿迭代法或BFGS法。为了提高求解效率和确保得到一个正确的解，用户能够指定形成新的有效刚度矩阵的求解步和指定进行平衡迭代的求解步。在进行增量求解之前，常系数结构矩阵(即线性有效刚度矩阵，可能用到的线性刚度、质量和阻尼矩阵)和载荷矢量，都已被组装，并贮存在低速存贮器上。而且，在逐步求解过程中，线性有效刚度矩阵只是因系统的非线性才作修正。方程的求解方法结构矩阵的有效存贮和平衡方程的有效求解是一个重要的问题。为了获得最大的解题能力，存贮方式需要优化。为了减少总

12、的求解费用，也需要有效的方程求解方法。在程序中，用了紧凑存贮方式，所有结构矩阵都按一维数组贮存，同时，只处理矩阵轮廓线以下的元素。作为一个例子，图5-2表示了在三角分解前后一个典型刚度矩阵中的元素图形。应该注意到在方程求解过程中，在轮廓线以内的零元素，一般不保持为零，所以这一些零元素要贮存，而所有轮廓图舸始刖度见阵的元索图52在按块存贮的刚度矩阵中典型的元素图线以外的元素都不必考虑。图5-2也表示在芯外存贮求解中，刚度和质量矩阵所用的块存贮方式。块数和每一块的列数是在总装过程之前由程序来计算，它取决于在矩阵总装和时间积分阶段中，全部可供利用的高速存贮容量的大小。应该注意到每块的列数都有变化，因

13、而供解题时使用的高速贮存器得到了最充分的利用。用芯外存贮的线性方程求解器Colsol来获得方程的解。这个方程基本上采用正定对称方程组的高斯消去法，但其系数矩阵按列运算。在求解过程中所完成的块运算如图5-3所示。在所有各类分析中也就是在线性、非线性、静力和动力分析中，都用同样的解法。同时，这种算法是由刚度矩阵(或有效刚度矩阵)的LDL分解、(有效)荷载矢量的消去和回代所组成。例如，在线性静力分析中，方程为KU=R，程序对U求解用：K=LDL(5-1)LV=R(5-2)DLtU=V(5-3)这里L和D分别是下三角矩阵和对角线矩阵。BAI|时ft快b靶磁的K-K*LNH1.图5-3刚度矩阵因子分解时

14、的块运算程序组织与框图在有限元前后处理部分中，将平面问题与空间问题分开编程，可单独运行。位移及应力分析部分适合于各种单元情况，通过数据文件输入和输出数据，该部分程序是在B0PACE3D程序的基础上开发而成的。在位移及应力分析程序的开发中，主要修改了源程序中许多语句错误，重编和新增了数据输入和输出部分的源程序，以适合于有限元前后处理程序的数据格式要求。另外，为满足矿岩自然崩落模拟的要求，自行开发增设了崩落判据程序。整个程序源代码近1万行，总框图如图54所示。模型关键控制点坐标输入第1超单元对应关键点号和柄料号等:.超理元釁遵爭吟1Yes!输出结果程序总控制卷数输夭材料雇性输灭鎚用元总邈U网格划努

15、数据输入有限元网格自动划外删除崩落单元信息第违寒吒剖沁迫剖归吃嵋因子网格的图形显示合理:单元剛度矩阵计篦形成总刚度矩阵载荷距阵计算迭代求解位移、应变和应力形成有限元甘析的节点及单元数据文件.蹩新单无編号、3:曲号图5-5网格自动划分程序框图图5-4程序流程总框图整个粘弹塑性有限元程序求解过程主要分为如下七个不同的阶段:有限元网格的自动形成有限元网格的自动划分程序主要针对岩土工程问题编写的。B0PACE3D源程序中的有限元前处理极不完善，是将所有节点和单元按一定的顺序进行一个一个的数据输入。考虑到工程实际计算模拟模型一般都较大，节点和单元少则几千，多则几万甚至几十万的数据输入量。如此繁琐而艰巨的

16、数据录入工作，给前期建模带来了极大的难度，且随着数据输入量的增大极易出错。有时尽管进行了反复的校核，但仍可能留下人为的错误而无法觉察到，这种错误的存在势必导致不正确的计算结果。如果在计算过程中发现计算结果不尽合理而要改变网格的划分，这时全部数据必须重新填写和输入，这样将导致工作量成倍的增加。鉴于BOPACE3D前处理的数据录入极不方便和难以确保数据的正确性，本文针对三维有限元划分网格的特点，自行开发了基于空间六面体8结点等参单元三维有限元网格自动划分的岩土工程通用程序，解决了BOPACE程序前处理难题，大大地提高了前期数据录入工作效率和数据的可靠度。网格自动划分程序的框图如图5-5所示。有限元

17、的图形显示有限元网格及边界条件的图形显示很重要，它可以直观显示有限元网格自动形成的结果，直观了解网格形成的正确性及合理性，指导使用人员修改有限元网格自动形成程序的输入参数。本文使用FORTRAN90加油站的图形功能及OpenGL三维图形函数库来实现图形显示，在FORTRAN90加油站编译系统上编译,有限元网格的图形显示程序框图如图5-6所示。尽管AutoCAD具有很强的绘图功能，能实现图形的显示和绘制，但得通过数据文件转换图形数据，然后使用Lisp语言或C语言绘图，如此处理存在许多问题。在FORTRAN语言内部通过OpenGL三维图形函数库绘图具有很大的方便性，不必传递图形数据，动态修改图形也

18、很方便。其次工程设计人员比较熟悉FORTRAN语言，使用简单、方便。OpenGL是硬件和图形的软件接口，实际上就是一个三维图形和模型库。由于它在三维真实感图形制作中具有优秀的性能，所以，诸如Microsoft、SGI、LBM、DEC、SUN、HP等在计算机市场占主导地位的大公司，都将它作为自己的图形标准，从而使之成为新一代的三维图形工业标准。OpenGL不仅仅是一个图形库，更是一个API(ApplicationProgrammingInterface)。它本身是一个与硬件无关的编程接口，可以在不同的硬件平台上得到实现，也正是因为如此，OpenGL中没有包含处理窗口和用户输入的命令。在OpenG

19、L中不提供三维造型的高级命令，虽然OpenGL也是通过基本的几何图元点、线、多边形来建立物体模型，但更确切地说它应该被称为新一代的三维图形开发标准。现在有很多优秀的三维图形软件，如3Dmax等可以方便地建立模型，但难以对其进行控制。把这些模型转化为用Fortran语言编写的OpenGL程序，就可以随心所欲地控制这些模型，制作CAD,制作三维动画，实现虚拟仿真，这些都使我们制作三维真实感图形更方便、更真实。OpenGL可以制作各种各样的三维图形，方便地实现三维图形的交互操作。OpenGL提供的图形操作有：绘制物体、变换、着色、光照、反走样、混合、雾、位图和图像、纹理映射和交互操作和动画。C语言及

20、FORTRAN语言可以很方便地调用OpenGL图形函数，实现图形显示，特别是OpenGL提供了许多坐标变换函数，绘制三维图形比较方便。因此，用FORTRAN语言绘制三维图形，选用OpenGL图形库是惟一较好的途径。而仅仅使用FORTRAN语言内部图形函数绘制平面图形是可以的，但绘制三维图形就非常困难。图5-6网格自动划分程序框图形成节点及单元数据文件网格自动生成后，对网格节点进行编号，形成有限元分析所需要的节点及单元数据文件。常数结构矩阵的装配在进行KUR求解之前，计算单元刚度矩阵，并组装成总刚度矩阵K，并贮存在硬盘上。此外，计算并贮存有效线性结构刚度矩阵。载荷矢量计算对每个时间步，计算外部施

21、加的荷载矩阵R并贮存到硬盘上。逐步求解采用迭代的方法计算位移、应变和应力程序，该程序是在BOPACE3D程序的基础上开发而成的，是有限元程序的重要组成部分，其功能强大、适用范围广。经过对源程序和理论文本的分初始化迭代变垦6P=KJ6Q计算总应变増量必总应变中分霸出譚.塑和蠕应变，并计算应力卜1耳畑zh1川3更新頤落数据Yes输出该増量步崩落数据F图5-7迭代求解位移、应变和应力程序框图图5-8崩落判据程序框图析和理解，得出位移、应力与应变有限元程序框图如图5-7所示。崩落判据程序框图崩落判据程序是根据对三矿区岩体特性进行研究后，经理论推导，得出相应公式，并在此基础上将其转化为程序的，崩落判据程

22、序框图如图5-8所示。5.3Flac3D原理简介5.3.1空间导数的有限差分近似三维快速拉格朗日法采用了混合离散方法，区域被划分为常应变六面体单元的集合体，而在计算过程中，程序内部又将每个六面体分为以六面体角点为角点的常应变四面体的集合体，变量均在四面体上进行计算，六面体单元的应力、应变取值为其内四面体的体积加权平均。如图5-9所示一四面体，节点编号为1到4第n面表示与节点n相对的面，设其内任一点的速率分量为U,i则可由高斯公式得：VA4图5-9四面体0udVi,jVndSj5-4)其中V为四面体的体积，S为四面体的外表面，n为外表面的单位法向向量分量。对于常i应变单元，u0为线性分布，n在每

23、个面上为常量，由式(5-4)可得：ii5-5)0140u=Xu1n(l)S(l)i,j3Vij1二i其中，上标1表示节点1的变量，(1)表示面1的变量。5.3.2运动方程三维快速拉格朗日法以节点为计算对象，将力和质量均集中在节点上，然后通过运动方程在时域内进行求解。节点运动方程可表示为如下形式：0(5-6)au1f1(t)i=ratm1其中F1(t)为在t时刻1节点的在i方向的不平衡力分量，可由虚功原理导出。m1为1节点i的集中质量，在分析动态问题时采用实际的集中质量，而在分析静态问题时则采用虚拟质量以保证数值稳定，对于每个四面体，其节点的虚拟质量为：TOC o 1-5 h zar1m1=Lmax(i)S(1)_2,i=1,3(5-7)9Vi其中a1=K+43G，K为体积模量，G为剪切模量。式(5-7)的前提是计算时步At=1。任一节点的虚拟质量为包含该节点的所有四面体对该节点的贡献之和。将上式左端用中心差分来近似，则可得到：5-8)U1(t+At)=U1(t-At)+F!(t)-ati2i2m1应变、应力及节点不平衡力三维快速拉格朗日法由速度来求某一时步的单元应变增量，如下式有了应变增量，100As=(u+u)Atij2i，jj，i即可由本构方程求出应力增量：Ao=H(o，As)+A/ijijijijij(5-9)5-10)其中H为已知的本构方程，A/为大变形情况下对应力所