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

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