搜索
第九讲岩石力学的数值模拟分析一、概述二、解析方法三、岩石力学中的大规模计算——并行算法四、不连续变形分析五、无单元法六、界面元方法一、概述岩石力学在上个世纪的突出进展在于数值计算方法，即计算机的计算方法在岩土工程中的应用获得了巨大的进展。近数十年来，我国已开展了大量的土木岩土工程。有些工程已属于世界第一流的水平，如水电工程中已建立了大量的混凝土高坝和堆石坝，如世界第三高拱坝二滩拱坝、龙羊峡高坝。现今可用于对岩体工程结构进行力学分析的数值方法多种多样，每一种方法有其针对性和特点，对一个具体的问题用数值模拟方法进行分析时，应选择一种最适合该问题的方法进行研究。数值模拟方法的选择，取决于研究对象即岩体工程结构的岩石力学性质和数值模拟的目的。严格地讲，岩石除具有弹性性质外，还具有塑性性质和粘性性质，只不过在特定情况下，某种性质占主导地位而已。在岩体工程实践中，硬岩及应力水平不甚高的中硬岩，其力学性质主要呈现为弹性或弹塑性；高应力环境下的软岩，其力学性质主要呈现为塑性或粘塑性；对于服务时间较长的地下工程，岩石极软或软且应力水平很高，则在计算分析中不能忽视岩石的流变性质。为了达到了解整个岩体工程系统的应力及变形规律的目的，各种数值方法均可采用，但以弹、塑性有限单元法或拉格朗日元法最为适宜。这两种方法的单元划分灵活，计算所需参数较少且易获得，软件也易于得到，成本较低。局部工程结构的应力及变形分析，若岩石中硬以上，则各种方法均可采用；若岩石软弱，则宜采用能进行大变形分析的拉格朗日元法；若岩体可能发生非连续破坏，则宜采用离散单元法、非连续变形法或流形元法或大变形分析的拉格朗日元法。要实现了解岩体结构的破坏特征及动态破坏过程的目的，则只能采用离散单元法、非连续变形法或流形元法，因为这些方法就是针对岩体介质的非连续性而提出的。此外，对于一个具体问题，是进行平面分析或是进行三维分析，也需做出恰当选择。严格地讲，所有的问题都是三维的，但如果采用平面分析既能达到目的，计算结果误差也不大，为了降低费用和快速方便起见，则以采用平面分析为宜。反之，则应采用三维程序代码进行计算分析。表1给出了几种数值模拟方法所依据的基本原理、求解方式、离散化方法及其适用条件，可供选择模拟方法时参考。数值模拟方法基本原理求解方法离散方式适用条件有限单元法最小势能原理解方程组全区域划分单元岩石中硬以上，小变形，岩体不会发生非连续性破坏如滑动、转动、分离等边界单元法Betti互等定理解方程组边界上划分单元同上离散单元法牛顿运动定律显式差分按结构弱面分布特征划分单元岩石中硬以上，低应力水平，大变形，岩体沿弱面发生非连续性破坏非连续变形法最小势能原理解方程组按主要结构弱面实际情况划分单元大变形，岩体发生非连续性破坏数值流形法最小势能原理解方程组全区域划分单元中硬以上岩体的连续或非连续变形拉各朗日元法牛顿运动定律显式差分全区域划分单元岩石软弱，大变形，岩体的破坏以变形为主表1几种数值模拟方法的对比一、概述二、岩石力学中的大规模计算——并行算法三、不连续变形分析四、无单元法五、界面元方法二、岩石力学中的大规模计算——并行计算近30年计算机发展很快，也推动了有限元的发展，并已应用于实际工程问题。但在遇到实际工程问题时，就有了困难，主要是实际工程量的DOF（自由度）太大。第一，工程对象几何形状复杂，如地下厂房的构造，尾水管内墩子，拱坝和地基往往也涉及很大的范围。第二，在求解一些角度点，如镇墩、排水孔结构等，自由度往往都很大。第三，在考虑岩体的开裂、岩体中的裂隙等，则要求更为庞大的自由度。这些都对计算能力提出了很高的要求，虽然现在的计算机的计算能力有了很大的提高，但是要完成更加复杂的计算任务，为了提高计算的规模和减少计算时间，就需要进行并行计算。几乎所有的大型计算机都是并行机。并行计算对于超大型的数值模型运算已经是一个必不可少的工具。并行计算在国外得到了很大的发展，自从美国国家航天局（NASA）的A.K.Noor在1975年发表第一篇关于有限元并行性计算的文章以来，有限元并行处理技术几乎与并行计算机同步发展，并且由算法研究发展到了算法、软件和硬件相结合的研究。在国内，并行计算还处于发展阶段，而且硬件相对比较落后。在硬件上主要基于向量机、分布式并行机和共享存储式并行机；在内容上，涉及范围比较广，但在系统性和深度上有待进一步发展，软件的开发距离实际应用和商品化还很远；专门针对有限元并行计算的硬件研究很少。（一）并行计算的基本概念并行计算机并行计算是指将计算任务分为n份，并将其分配给n个计算节点（计算节点数对网络机群而言指工作站台数；对大型并行机而言指CPU个数）同时进行计算。并行实际上包含了两个概念：同时性（Simultaneity）和并发性（Concurrency）。同时性是指两个或者多个事件在同一时刻发生；并发性是指两个或者多个事件在同一时间间隔内发生。对于并行计算机，按照实际应用，可以分为如下的4类[1]：大规模并行机（MassivelyParallelProcessors，MPP）对称多处理器系统（SymmetricMultiProcessors，SMP）分布存储系统（DistributedSharedMemorysystem，DSM）网络工作站机群（ClusterOfWorkstations，COW）在上面4种并行机中，COW和MPP具有类似的地方，实际上，当前MPP和COW之间的界限正在逐渐模糊。例如IBMSP2被视为MPP，但它却具有机群结构。由于COW的性能价格比远优于MPP，系统扩充方便，网络连接形式多种多样，所以COW是发展可扩放并行计算的主流趋势。并行算法及并行程序开发并行算法是区别于串行算法的另外一大类型的算法，它是适合于各种并行计算机上求解问题和处理数据的算法。并行程序开发是对给定算法构造并行程序的活动，它要求算法设计者和计算机系统的体系结构的设计者进行广泛的交互。并行程序的设计不仅编程困难，而且调试和分析更加困难。目前，并行程序设计在各方面都处于一个初级阶段，它不支持一个成熟、稳定和通用的并行程序开发环境。并行开发环境指的是软件环境，它对并行计算的影响比一般串行机要大得多。软件对计算性能的影响的差别可达几个数量级。并行开发环境包括两部分内容：操作系统、通信平台、编译和调试工具以及性能测试软件等。其中最主要的是操作系统和通信平台。可以用于并行计算的操作系统主要有Windows（NT或者2000）和Unix，其中与Unix内核基本相同的自由软件Linux在并行计算中占优绝对的统治地位，而比较流行的通信平台有P4、PVM、Express、PARMACS和MPI等。（二）有限元法的并行计算一般来说，有限元计算包括如下的步骤：（1）数据准备阶段，也就是前处理阶段，包括结构的离散化、单元信息和节点坐标的生成。（2）单元分析阶段，包括计算单元刚度矩阵、单元荷载向量等。对于动力分析，还需要计算单元质量矩阵。（3）生成总体刚度矩阵和总荷载向量。（4）约束处理。（5）方程组求解，得到全部节点的位移。（6）后处理。在上面的计算中，（1）、（4）和（6）阶段的计算花费的时间很短，而（2）、（3）和（5）阶段则是有限元分析的关键所在。尤其是（5），将占整个有限元分析的70％的计算量。因此，并行计算主要是针对这3个步骤进行的。用于多机并行系统的有限元并行分析策略对于三维有限元的分析，如果在自带存储器的多机并行系统上，可以采用如下的两种策略：SBS（SubstructBySubstruct）策略对要求解的结构采用某种剖分策略，比如，类似于子结构的方法，使得各个机器所承担的部分刚度/质量矩阵可以进行“孤立”自由度（内部自由度）和“公共”自由度（界面自由度）分块。“孤立”自由度将直接在“本地”进行“消去”，凝聚矩阵提交给“中心机”进行装配并完成求解。然后由各个“节点机”完成消去未知量的换算。所采用的“剖分策略”，一般以孤立自由度数目的总和最大为优化目标，兼顾及“同时完成任务”。采用这种策略，并行化计算是建立在子结构一级的水平上。EBE（ElementByElement）策略EBE策略的基本思想就是将一个总体矩阵的向量积转化到一组单个矩阵的向量积的计算。在有限元计算中，总体矩阵就是由多个单元矩阵迭加而成，因此，非常适合使用EBE策略来进行求解。EBE策略是一种既能节省存储量、又能适宜于并行计算而且容易实现的方法。静力分析在有限元的静力分析中，主要是采用两大类算法：直接并行算法和迭代并行算法。对于直接并行算法，一种是采用子结构直接并行算法，其基本思想就是SBS策略的基本思想。另外一种直接解法是仅对占很大计算量的有限元方程组进行并行计算，其他计算步骤仍然采用串行算法。例如，并行求解方程组可以采用并行Gauss消去法、并行Cholesky分解等方法。周树荃，邓绍忠对直接并行算法进行了很深入的研究，并在YH－I向量机上进行了大量的数值试验。对于迭代并行算法，一般都采用同步运算或者基本完全异步控制，因此，迭代并行算法易于发挥并行机的向量处理功能。目前常用的算法有基于EBE策略的并行预处理共轭梯度法（PreconditioningConjugateGradient，PCG）、并行GMRES方法以及基于多色理论的超松弛迭代（SuccessiveOver-RelaxationMethod，SOR）法等。动力分析TLDL动力响应分析在这方面内容中，研究最多的是直接积分法，并取得了很大的成果。例如，高昀采用分解方法，在Transputer及网络机群环境下，实现了基于Newmark方法的的动力响应分析。结构的动力特性研究在这一方面，相对于动力的响应分析来说，研究的相对较少。T.Hwang和I.D.Parsonsp网格法在MIMD并行计算机上求解特征值问题。国内在这一方面也做了大量的工作，主要是采用Lanczos方法对其进行求解。另外，徐甲同采用QIF分解方法对多处理机上矩阵特征值的并行计算进行了研究。发展方向采用可以达到完全并行的EBE策略，可以完全突破大规模计算对串行计算机的内存和计算速度的限制。理论上而言，只要参与计算的节点足够多，就可以对任意规模的数值计算问题进行求解。虽然采用EBE策略的预处理共轭梯度法是解决有限元并行计算很有效的方法，也是有限元并行计算发展的一个重要方向。尤其是在非线性问题中，EBE方法的效率是很高的。但是目前的各种预处理方法，不仅仅是EBE策略的预处理，就是对整体矩阵进行的预处理，都很难从理论上说明那一种更有效。如何给出一个有效的预处理方法仍然是一个需要解决的问题。另外，EBE方法在非线性问题中的收敛性也是需要进行研究的问题。一、概述二、岩石力学中的大规模计算——并行算法三、不连续变形分析四、无单元法五、界面元方法三、不连续变形分析通常状态岩体是被节理和断层等构造切割成的不连续介质，这些不连续介质影响着岩体结构的滑移和破坏，不连续介质数值模型即是为了解决这种地质问题而产生的。非连续变形分析,即DDA(DiscontinuousDeformationAnalysis)，是由石根华博士提出的分析块体系统运动和变形的一种数值方法。它以自然存在的节理面或断层切割岩体形成不同的块体单元，以各块体的位移作为未知量，通过块体间的接触和几何约束形成一个块体系统。在块体运动的过程中，严格满足块体间不侵入。将边界条件和接触条件等一同施加到总体平衡方程。总体平衡方程是由系统的最小势能原理求得。求解方程组即可得到块体当前时步的位移场、应力场、应变场及块体间的作用力。反复形成和求解总体平衡方程式，即可得到多个时步后以至最终平衡时的位移场及应力场等情况以及运动过程中各块体的相对位置及接触关系。因此，DDA方法可以模拟出岩石块体的移动、转动、张开、闭合等全部过程并据此判断岩体的破坏程度、破坏范围，从而对岩体的整体和局部稳定性作出正确的评价。DDA方法在三峡船闸高边坡稳定分析以及锦屏的稳定分析中得到应用。1、近期的发展块体系统非连续变形分析是近年来发展起来的一种崭新的数值计算理论。在满足弹性理论的基本方程条件下能够反