第九章 岩石力学数值模拟分析

1、第九讲第九讲 岩石力学的数值模拟分析岩石力学的数值模拟分析 岩石力学在上个世纪的突出进展在于数值计算方法，即计算机的计算方法在岩土工程中的应用获得了巨大的进展。近数十年来，我国已开展了大量的土木岩土工程。有些工程已属于世界第一流的水平，如水电工程中已建立了大量的混凝土高坝和堆石坝，如世界第三高拱坝二滩拱坝、龙羊峡高坝。 现今可用于对岩体工程结构进行力学分析的数值方法多种多样，每一种方法有其针对性和特点，对一个具体的问题用数值模拟方法进行分析时，应选择一种最适合该问题的方法进行研究。 数值模拟方法的选择，取决于研究对象即岩体工程结构的岩石力学性质和数值模拟的目的。 严格地讲，岩石除具有弹性性质外

2、，还具有塑性性质和粘性性质，只不过在特定情况下，某种性质占主导地位而已。在岩体工程实践中，硬岩及应力水平不甚高的中硬岩，其力学性质主要呈现为弹性或弹塑性；高应力环境下的软岩，其力学性质主要呈现为塑性或粘塑性；对于服务时间较长的地下工程，岩石极软或软且应力水平很高，则在计算分析中不能忽视岩石的流变性质。 为了达到了解整个岩体工程系统的应力及变形规律的目的，各种数值方法均可采用，但以弹、塑性有限单元法或拉格朗日元法最为适宜。这两种方法的单元划分灵活，计算所需参数较少且易获得，软件也易于得到，成本较低。 局部工程结构的应力及变形分析，若岩石中硬以上，则各种方法均可采用；若岩石软弱，则宜采用能进行大变

3、形分析的拉格朗日元法；若岩体可能发生非连续破坏，则宜采用离散单元法、非连续变形法或流形元法或大变形分析的拉格朗日元法。 要实现了解岩体结构的破坏特征及动态破坏过程的目的，则只能采用离散单元法、非连续变形法或流形元法，因为这些方法就是针对岩体介质的非连续性而提出的。 此外，对于一个具体问题，是进行平面分析或是进行三维分析，也需做出恰当选择。严格地讲，所有的问题都是三维的，但如果采用平面分析既能达到目的，计算结果误差也不大，为了降低费用和快速方便起见，则以采用平面分析为宜。反之，则应采用三维程序代码进行计算分析。 表1给出了几种数值模拟方法所依据的基本原理、求解方式、离散化方法及其适用条件，可供选

4、择模拟方法时参考。数值模拟数值模拟方法方法基本原理基本原理 求解方法求解方法离散方式离散方式适用条件适用条件有限单元有限单元法法最小势能原理解方程组全区域划分单元岩石中硬以上，小变形，岩体不会发生非连续性破坏如滑动、转动、分离等边界单元边界单元法法Betti互等定理解方程组边界上划分单元同上离散单元离散单元法法牛顿运动定律显式差分按结构弱面分布特征划分单元岩石中硬以上，低应力水平，大变形，岩体沿弱面发生非连续性破坏非连续变非连续变形法形法最小势能原理解方程组按主要结构弱面实际情况划分单元大变形，岩体发生非连续性破坏数值流形数值流形法法最小势能原理解方程组全区域划分单元中硬以上岩体的连续或非连续

5、变形拉各朗日拉各朗日元法元法牛顿运动定律显式差分全区域划分单元岩石软弱，大变形，岩体的破坏以变形为主表1 几种数值模拟方法的对比近30年计算机发展很快，也推动了有限元的发展，并已应用于实际工程问题。但在遇到实际工程问题时，就有了困难，主要是实际工程量的DOF（自由度）太大。第一，工程对象几何形状复杂，如地下厂房的构造，尾水管内墩子，拱坝和地基往往也涉及很大的范围。第二，在求解一些角度点，如镇墩、排水孔结构等，自由度往往都很大。第三，在考虑岩体的开裂、岩体中的裂隙等，则要求更为庞大的自由度。这些都对计算能力提出了很高的要求，虽然现在的计算机的计算能力有了很大的提高，但是要完成更加复杂的计算任务，

6、为了提高计算的规模和减少计算时间，就需要进行并行计算。几乎所有的大型计算机都是并行机。并行计算对于超大型的数值模型运算已经是一个必不可少的工具。并行计算在国外得到了很大的发展，自从美国国家航天局（NASA）的A.K.Noor在1975年发表第一篇关于有限元并行性计算的文章以来，有限元并行处理技术几乎与并行计算机同步发展，并且由算法研究发展到了算法、软件和硬件相结合的研究。在国内，并行计算还处于发展阶段，而且硬件相对比较落后。在硬件上主要基于向量机、分布式并行机和共享存储式并行机；在内容上，涉及范围比较广，但在系统性和深度上有待进一步发展，软件的开发距离实际应用和商品化还很远；专门针对有限元并行

7、计算的硬件研究很少。（一）（一）并行计算的基本概念并行计算的基本概念并行计算机并行计算机并行计算是指将计算任务分为n份，并将其分配给n个计算节点（计算节点数对网络机群而言指工作站台数；对大型并行机而言指CPU个数）同时进行计算。并行实际上包含了两个概念：同时性（Simultaneity）和并发性（Concurrency）。同时性是指两个或者多个事件在同一时刻发生；并发性是指两个或者多个事件在同一时间间隔内发生。对于并行计算机，按照实际应用，可以分为如下的4类1：大规模并行机（Massively Parallel Processors，MPP）对称多处理器系统（Symmetric MultiPr

8、ocessors，SMP）分布存储系统（Distributed Shared Memory system，DSM）网络工作站机群（Cluster Of Workstations，COW）在上面4种并行机中，COW和MPP具有类似的地方，实际上，当前MPP和COW之间的界限正在逐渐模糊。例如IBM SP2被视为MPP，但它却具有机群结构。由于COW的性能价格比远优于MPP，系统扩充方便，网络连接形式多种多样，所以COW是发展可扩放并行计算的主流趋势。并行算法及并行程序开发并行算法及并行程序开发并行算法是区别于串行算法的另外一大类型的算法，它是适合于各种并行计算机上求解问题和处理数据的算法。并行程

9、序开发是对给定算法构造并行程序的活动，它要求算法设计者和计算机系统的体系结构的设计者进行广泛的交互。并行程序的设计不仅编程困难，而且调试和分析更加困难。目前，并行程序设计在各方面都处于一个初级阶段，它不支持一个成熟、稳定和通用的并行程序开发环境。并行开发环境指的是软件环境，它对并行计算的影响比一般串行机要大得多。软件对计算性能的影响的差别可达几个数量级。并行开发环境包括两部分内容：操作系统、通信平台、编译和调试工具以及性能测试软件等。其中最主要的是操作系统和通信平台。可以用于并行计算的操作系统主要有Windows（NT或者2000）和Unix，其中与Unix内核基本相同的自由软件Linux在并

10、行计算中占优绝对的统治地位，而比较流行的通信平台有P4、PVM、Express、PARMACS和MPI等。（二）有限元法的并行计算（二）有限元法的并行计算一般来说，有限元计算包括如下的步骤：（1）数据准备阶段，也就是前处理阶段，包括结构的离散化、单元信息和节点坐标的生成。（2）单元分析阶段，包括计算单元刚度矩阵、单元荷载向量等。对于动力分析，还需要计算单元质量矩阵。（3）生成总体刚度矩阵和总荷载向量。（4）约束处理。（5）方程组求解，得到全部节点的位移。（6）后处理。在上面的计算中，（1）、（4）和（6）阶段的计算花费的时间很短，而（2）、（3）和（5）阶段则是有限元分析的关键所在。尤其是（5

11、），将占整个有限元分析的70的计算量。因此，并行计算主要是针对这3个步骤进行的。用于多机并行系统的有限元并行分析策略用于多机并行系统的有限元并行分析策略对于三维有限元的分析，如果在自带存储器的多机并行系统上，可以采用如下的两种策略：SBS（Substruct By Substruct）策略 对要求解的结构采用某种剖分策略，比如，类似于子结构的方法，使得各个机器所承担的部分刚度/质量矩阵可以进行“孤立”自由度（内部自由度）和“公共”自由度（界面自由度）分块。“孤立”自由度将直接在“本地”进行“消去”，凝聚矩阵提交给“中心机”进行装配并完成求解。然后由各个“节点机”完成消去未知量的换算。所采用的“

12、剖分策略”，一般以孤立自由度数目的总和最大为优化目标，兼顾及“同时完成任务”。采用这种策略，并行化计算是建立在子结构一级的水平上。EBE（Element By Element）策略 EBE策略的基本思想就是将一个总体矩阵的向量积转化到一组单个矩阵的向量积的计算。在有限元计算中，总体矩阵就是由多个单元矩阵迭加而成，因此，非常适合使用EBE策略来进行求解。EBE策略是一种既能节省存储量、又能适宜于并行计算而且容易实现的方法。静力分析静力分析在有限元的静力分析中，主要是采用两大类算法：直接并行算法和迭代并行算法。对于直接并行算法，一种是采用子结构直接并行算法，其基本思想就是SBS策略的基本思想。另外

13、一种直接解法是仅对占很大计算量的有限元方程组进行并行计算，其他计算步骤仍然采用串行算法。例如，并行求解方程组可以采用并行Gauss消去法、并行Cholesky分解等方法。周树荃，邓绍忠对直接并行算法进行了很深入的研究，并在YHI向量机上进行了大量的数值试验。对于迭代并行算法，一般都采用同步运算或者基本完全异步控制，因此，迭代并行算法易于发挥并行机的向量处理功能。目前常用的算法有基于EBE策略的并行预处理共轭梯度法（Preconditioning Conjugate Gradient，PCG）、并行GMRES方法以及基于多色理论的超松弛迭代（Successive Over-Relaxation

14、Method，SOR）法等。动力分析动力分析TLDL动力响应分析 在这方面内容中，研究最多的是直接积分法，并取得了很大的成果。例如，高昀采用 分解方法，在Transputer及网络机群环境下，实现了基于Newmark方法的的动力响应分析。结构的动力特性研究 在这一方面，相对于动力的响应分析来说，研究的相对较少。T.Hwang和I. D. Parsonsp网格法在MIMD并行计算机上求解特征值问题。国内在这一方面也做了大量的工作，主要是采用Lanczos方法对其进行求解。另外，徐甲同采用QIF分解方法对多处理机上矩阵特征值的并行计算进行了研究。发展方向发展方向采用可以达到完全并行的EBE策略，可

15、以完全突破大规模计算对串行计算机的内存和计算速度的限制。理论上而言，只要参与计算的节点足够多，就可以对任意规模的数值计算问题进行求解。虽然采用EBE策略的预处理共轭梯度法是解决有限元并行计算很有效的方法，也是有限元并行计算发展的一个重要方向。尤其是在非线性问题中，EBE方法的效率是很高的。但是目前的各种预处理方法，不仅仅是EBE策略的预处理，就是对整体矩阵进行的预处理，都很难从理论上说明那一种更有效。如何给出一个有效的预处理方法仍然是一个需要解决的问题。另外，EBE方法在非线性问题中的收敛性也是需要进行研究的问题。 通常状态岩体是被节理和断层等构造切割成的不连续介质，这些不连续介质影响着岩体结

16、构的滑移和破坏，不连续介质数值模型即是为了解决这种地质问题而产生的。 非连续变形分析,即DDA(Discontinuous Deformation Analysis)，是由石根华博士提出的分析块体系统运动和变形的一种数值方法。它以自然存在的节理面或断层切割岩体形成不同的块体单元，以各块体的位移作为未知量，通过块体间的接触和几何约束形成一个块体系统。在块体运动的过程中，严格满足块体间不侵入。将边界条件和接触条件等一同施加到总体平衡方程。总体平衡方程是由系统的最小势能原理求得。求解方程组即可得到块体当前时步的位移场、应力场、应变场及块体间的作用力。反复形成和求解总体平衡方程式，即可得到多个时步后以

17、至最终平衡时的位移场及应力场等情况以及运动过程中各块体的相对位置及接触关系。因此，DDA方法可以模拟出岩石块体的移动、转动、张开、闭合等全部过程并据此判断岩体的破坏程度、破坏范围，从而对岩体的整体和局部稳定性作出正确的评价。DDA方法在三峡船闸高边坡稳定分析以及锦屏的稳定分析中得到应用。1、近期的发展、近期的发展 块体系统非连续变形分析是近年来发展起来的一种崭新的数值计算理论。在满足弹性理论的基本方程条件下能够反映出岩体变形的不连续性，它既具有有限元理论基础的严密性，又具有离散元法可计算块体大位移的特点，是一种很有发展前途的数值计算方法。它以新颖的理论思想，严密科学的构思，广泛应用前景相融为一

18、体。作者在三维理论的建立、完善的基础上，结合三峡船闸高边坡的稳定性分析进行了更深入细致的研究，全面发展完善和丰富了块体系统非连续变形分析理论，取得了如下主要研究成果。（1）提出三维非连续变形分析理论，并对其理论进行完善。并从理论上论证和分析三维非连续变形分析理论的正确性，以12个变形参数来描述块体系统的刚体运动和变形。建立总体平衡方程和各类子矩阵。（2）块体间高效而正确的接触判定是实现非连续变形分析计算的关键所在。提出利用以参考面法、参考线法来判定二维和三维块体间的接触。（3）参考面（线）是块体边界上的某些平面（线），不必重新生成。根据块体间各平面的相对关系，建立了确定参考面的判据，参考面法不

19、仅适合对凸形块体的接触判定，而且对凹形块体的接触同样适合。（4）进一步提出可采用形象法来寻找参考线、参考面，并给出了具体的寻找方法和步骤。利用形象法亦可迅速完成参考面的识别，该方法形象、直观，这是一种高效率的搜寻方法。（5）建立了块体静动力学约束条件块体间无侵入发生及不承受拉伸力的控制标准和位置锁定原则。（6）实际工程中非连续面的大小是有限平面，它们在空间分布上具有很大的随机性。针对这一情况，作者分析了目前处理节理网络技术方面的不足，提出利用拓扑学单纯同调原理解决块体的形成问题。根据空间各有限平面，可以为有限元法、离散元法 块体理论以及其他要求具有一般形状的数值方法提供模型数据。该方法不仅可以

20、生成规则的块体，而且还可以生成形状十分复杂的多边形、多面体，甚至是带孔洞的块体。（7）根据总刚矩阵 的特点，提出了一种高效三角分解方法图解法。这种借助于图解理论的非零元素存贮方法使得计算过程直观而高效。通过图形和连接矩阵的配合使用，使整个过程直观明了，具有存贮量小、计算量少、舍入误差小等优点。利用图解技术，事先对块体优化编号，确定存贮量，并使三角分解过程中产生的新元素个数最少，从而达到减少数据存贮和计算，提高解题效率的目的。K（8）根据非连续变形分析理论，编制和开发了相应的程序。并结合三峡船闸高边坡稳定性分析实例进行了研究和分析，其计算结果与实际观测结果大体接近。（9）从理论上解决了非连续变形

21、分析理论与有限元、边界元及传统极限平衡方法的耦合，并给出耦合原理及耦合的公式。（10）针对非连续变形分析法存在的问题，提出解决方案。时步的大小事关非连续变形分析理论计算结果的正确性，通过众多的计算表明仅当时步小于0.005s时，计算结果是稳定可靠的。（11）针对非连续变形分析理论自身的特点及高边坡滑坡位移复杂性，从理论上深刻分析、研究了边坡滑坡时产生混沌的机理及过程，并提出防止边坡进入混沌状态的实际措施。（12）将不连续变形分析与AutoCAD动态技术结合，使块体系统非连续变形分析法的计算实施过程形象、直观，完成了块体系统的变形或破坏全过程的模拟。2、不连续变形及其特点不连续变形及其特点 非连

22、续变形分析以各个块体的位移作为求知量，通过块体间的接触和几何约束形成一个块体系统，块体单元受非连续面的控制，在块体运动的过程中单元之间可以接触，也可以分离。在块体运动过程中，满足块体间不侵入和不承受拉伸力的条件，总体平衡方程是由系统的最小势能原理求得的。求解方程组即可得到块体当前时步的位移场、应力场、应变场及块体间的作用力，反复形成和求解总体平衡方程式，即可得到多外时步后块体的变形，应力及位移情况，也可求得块体系统最终达到平衡时应力场、位移场等情况以及运动过程中各块体的相对位置及接触关系。当然，非连续变形分析方法也和其它任何一种数值方法一样，有其一定的适用范围和局限性，这主要表现在：（1）岩体

23、的种类繁多，其性质极为复杂。在通常情况下，非连续岩体的力学表现为非均质、各向异性、非线性等复杂的特性，随着应用工程环境的不同，使岩体数值计算有许多非确定因素。这对非连续变形分析在对岩体参数（ ， ）取值时带来麻烦，岩体参数的取值直接影响到计算结果的正确性。目前，非连续变形分析的数值计算中，一般假定岩体是弹性体，对塑性、粘性体等尚不适用。c（2）计算时步的大小肯定会影响非连续变形分析的最终结果，时步过大将导致计算结果与实际测量值相差较大，从更深理论上研究时步的合理取值对非连续变形分析的计算具有较大的指导意义。计算时步的大小与实际问题联系起来，合理地考虑时间因素也是今后需要研究的问题。（3）非连续

24、变形分析方法比较适合模拟硬岩，对于软岩或软硬相间的情况在判定块体接触时会遇到困难。也就是说非连续变形分析比较适合于块体自身小变形的情况。（4）在非连续变形分析实施的过程中，有大量的数据需要处理，需要大量的计算机时间来存贮、转换和调用，耗用了大量的计算机CPU时间，应优化和改进计算方法，缩短计算时间。3、DDA今后的发展方向今后的发展方向 DDA是近10年来迅速发展起来的一种较新的岩土数值计算方法，已于1996年在美国召开过一次国际性研讨大会。之后，在日本、英国相继举行了研讨会。若能将非连续变形分析和已发展成熟的有限元、边界元以及极限平衡等方法相结合，发挥各自的优势，则会使这一方法更趋于完善。对

25、非连续变形分析方法的研究，应重点改进前面提及的不足之处。除上述工作外，正在进行和将要进行的工作还应包括如下的内容：（1）继续完善三维非连续变形分析理论以及计算程序的实现仍然是今后一个时期内需要解决的问题，优化计算机程序和设计方法是提高解题效率和推广非连续变形分析法应用的前提。（2）从理论和计算机程序设计上解决颗粒形状单元在不同颗粒级配条件下的非连续变形分析。对常见的几种概率分布形式下的颗粒级配沿途进行非连续变形分析法的模拟计算。（3）渗流问题是影响岩体稳定的重要因素之一，因此有必要在非连续变形分析理论和相应的程序中加入解决渗流问题的功能，加强渗流问题与非连续变形分析法的结合具有重要的意义。（4

26、）加强流行、分形、混沌等新兴学科在非连续变形分析中的应用与推广，进一步研究图形理论与非连续变形分析法的结合。DDA方法作为一种新型的岩土数值计算方法，有着广阔的应用前景。邬爱箐将DDA应用于三峡高边坡，梁国平提出的LDDA也是当前研究的热门课题，LDDA（Lagrange Discontinuous Deformation Analysis）方法即具有拉格朗日乘子的不连续变形方法，是模拟非线性界面的有效方法，其在界面上设置接触单元来模拟张裂和滑动。拉格朗日乘子即是接触面上的正压力和切向力。由于区域分解方法的使用，使求解的效率大大提高；由于使用了DDA方法中接触判断准则，因此不必定义滑移面。 无

27、单元法是一种新的数值计算方法，其特点是采用滑动最小二乘法所产生的光滑函数来近似场函数，计算形函数，从而只需计算域的几何边界及计算点，摆脱了单元限制，大大简化了前处理工作。由于提供了场函数的连续可导近似解，在材料分析中，使得位移及应力、应变计算结果均连续，不需进行后处理修匀。无单元法的结点生成非常容易，根本不存在其它方法所遇到的网格重构问题，因此在开裂计算中将有很好的应用前景。总之，无单元法保留了有限元的一些特点，克服了有限元的不足，适于进行岩土工程数值模拟，尤其便于跟踪裂纹扩展，提供了岩土工程数值模拟的新途径。无单元法的最早思想由Nayroles等于1992年提出，并被称之为虚拟单元法(Dif

28、fuse Element Method)。但Nayroles的近似比较简单，使用了低阶积分，边界条件的引入也不准确，而且在拟合函数的求导中忽略了一项，计算比较粗糙。Belytschko等对Nayroles的方法进行了改进，提出了无单元伽辽金法(Element-free Galerkin Method)，采用了高阶高斯积分，引入了被Nayroles忽略的项，并用拉格朗日乘子引入本征边界条件。Lu等又将无单元伽辽金法作了进一步改进，在滑动最小二乘法中采用了正交基，大大加快了计算效率和精度。周维垣等4用罚函数法引入了本征边界条件，并对一些关键问题作了探讨和改进，使得无单元法得到了与有限元法类似的方程

29、组形式，并提高了计算精度。 由于无单元法具有前后处理简单、计算精度高、在计算过程当中便于增减结点等优点，已被广泛应用于板弯曲计算，弹塑性问题分析以及线弹性开裂分析中。1、基本原理基本原理 滑动最小二乘法由Lancaster提出，但直到近年才由Nayroles等将其应用于求解边值问题。它是用加权最小二乘法来近似场函数的一种方法，其详细原理参考文献。为了克服滑动最小二乘法在系数矩阵形态较差时计算精度低的缺点，在实际计算中采用了将基函数正交化的方法。由变分原理即：平衡位移使系统总势能取驻值，并采用罚函数法引入位移边界条件，在位移边界上作用刚性很大的弹簧，并将已知位移作为弹簧的初始拉伸位移。无单元法计

30、算流程如图N.7。开始生成节点并计算节点影响半径存储各网格所包含的节点号确定高斯点位置及其积分权检索高斯点的影响节点计算高斯点对整体平衡方程的贡献并将其集入整体平衡方程解平衡方程求应力、应变等所需变量并进行结果分析结束高斯点上逐个计算图N.7 无单元法的计算流程2、方程组解法方程组解法 由无单元法平衡方程的集成过程可以看出，其刚度矩阵具有对称正定性和稀疏性。如果节点号编制合理，无单元法还可以做到刚度矩阵非零元素的带形分布。但有两个问题：（1）带宽较大因为无单元法中形函数是由最小二乘法得到，在求每个高斯点上的函数值时，为了使最小二乘法的正规方程系数矩阵非奇异，通常应使每个高斯点有较多的影响节点，

31、而不象有限元法中只有位于同一单元中的节点才相互作用。这使得无单元法即使节点号编制合理，其带宽一般也是较大的。（2）不能发挥无单元法的优越性 无单元法最大的特点是消除了单元，而只用节点进行计算，这使得它在某些方面的应用有很大的优越性。首先，它可以任意加密节点，而不需象有限元那样重新划分网格；其次，在追踪断裂扩展时很方便，可在裂尖不断布点，而不需重新划分网格。但在这种情况下，要做到带宽优化就很困难了，除非每加一个节点或每移动一个点都重新进行一次节点编号的优化，但如果实际操作起来，其计算量将是非常大的。为了克服这些缺陷，无单元法方程组的求解采用了一种迭代解法共轭梯度法，共轭梯度法是用来求解系数矩阵为

32、对称正定稀疏矩阵的大型线性方程组的一种有效方法，其基本理论及算法参见文献。 界面元又称之为刚体弹簧元、刚体元、界面应力元等，它是基于1980年日本东京大学Kaiwai教授和周维垣、杨若琼提出的刚体弹簧元模型而建立有限元计算的。它综合了目前岩土工程数值分析中连续介质变形体模型和不连续刚体模型的特点，其建模思路可概括为三个方面：（1）假定单元的变形累积在界面层，单元本身就可视为刚体，因而可采用分片刚体位移模式构造结构的应力场，使其逼近真实的位移场。（2）基于相邻块体单元的相对变形和不同的材料特性（非均质、各向异性、弹塑性、非线性以及粘性），导出界面应力公式并以此表征结构的应力场。（3）根据微分方程的弱解形式加权残量法或虚功原理或者放松界面位移连续条件的广义变分原理，建立相应的支配方程。1、界