数值计算方法

数值计算方法的共同思想

节理岩体数值计算方法及其应用(一)：方法与讨论 摘要 叙述了工程岩体的典型破坏类型、发生条件和判断方法，简单介绍了针对节理岩体的数值计算方法发展过 程、现状和特点，重点叙述了利用3DEC 基本功能生成三维随机节理网络和节理力学参数的随机赋值方法。 其中 三维节理网络技术可以非常逼真地模拟现实条件下各种复杂形式的节理分布，甚至可以直接利用节理的地质编录 资料在数值模型中生成确定位置上的确定性节理，在某种程度上一体化地实现了节理岩体的计算机模拟和数值计 算，促进了不连续岩体、岩石计算方法在应用技术方面的发展。 最后讨论了不连续力学数值计算技术的工程应用， 指出了合理运用数值计算方法解决工程...全部

  节理岩体数值计算方法及其应用(一)：方法与讨论 摘要 叙述了工程岩体的典型破坏类型、发生条件和判断方法，简单介绍了针对节理岩体的数值计算方法发展过 程、现状和特点，重点叙述了利用3DEC 基本功能生成三维随机节理网络和节理力学参数的随机赋值方法。
  其中 三维节理网络技术可以非常逼真地模拟现实条件下各种复杂形式的节理分布，甚至可以直接利用节理的地质编录 资料在数值模型中生成确定位置上的确定性节理，在某种程度上一体化地实现了节理岩体的计算机模拟和数值计 算，促进了不连续岩体、岩石计算方法在应用技术方面的发展。
  最后讨论了不连续力学数值计算技术的工程应用， 指出了合理运用数值计算方法解决工程问题的一般原则。 关键词 岩石力学，数值模拟，不连续体，离散元，随机，节理岩体 1 概述 数值计算已经被普遍应用于工程设计中解决各 种岩石力学问题。
  随着数值计算技术的发展，它在 很大程度上已经取代了一些传统的研究手段，使得 岩土工程中的许多岩石力学问题有可能较快地得到 解决。目前已经开发出来的岩土工程数值计算软件 及相关领域的工程软件，已为岩土工程设计和施工 中相关问题的研究提供了充分的选择余地。
  在目前 条件下如何正确选择和有效地使用这些软件，则成 为岩石力学数值计算实际应用中的一个现实问题。 对岩体几何特征的研究开始于 20世纪70 年代， 主要是采用概率统计方法研究节理的产状、间距和 迹长等几何参数的统计分布[1]以及岩体节理网络计 算机模拟等，其最终目的还是为研究节理岩体的力 学特性服务。
  然而，迄今为止，能正确地描述岩体 中的节理分布，并把它们与岩石力学数值计算方法 相结合的应用实例还很鲜见。 本文的研究，从岩体破坏类型及其发生条件入 手，首先界定结构面控制型破坏的发生条件，从而 规定了需要采用不连续力学方法进行岩石力学分析 方法的一般范围，利用3DEC 程序的基本功能直接 在3DEC 计算模型中模拟节理网络方法和节理力学 指标值的随机赋值方法，针对中国水电工程特点， 讨论了进行数值计算时需要考虑的有关岩体非连续 性问题。
   2 岩体破坏类型与发生条件 地表和近地表岩体工程如人工边坡、深度不大 的矿坑工程涉及的岩体初始应力水平一般都不高， 岩体的潜在破坏方式主要取决于岩体本身特性。对 于坚硬和中等坚硬岩石地区的这类工程，一般而言， 结构面对岩体的潜在破坏方式起决定性作用。
   地下工程围岩破坏可以概括地分为两大类：一 类是应力控制型破坏；另一类是块体控制型破坏(也 称结构面控制型破坏或重力型破坏)。这两类破坏的 岩体中各有其发生条件，前者主要发生在初始地应 力水平相对较高(相对于岩体承载能力)的岩体中， 后者则多出现在地表和地下浅埋岩体工程中。
   修建在以硬岩为主的岩体地下工程，其围岩的 初始应力水平在很大程度上决定了地下工程围岩二 次应力的量级，而岩体的承载能力取决于岩石强度 和结构面发育程度。这3 者之间的关系直接影响着 围岩的状态和潜在问题的性质[2]。
  图1 描述了这3 个因素与围岩破坏特征之间的一般关系。其中的左 侧第1 行描述了初始地应力水平，它是用初始地应 力中的最大主应力与岩石(硬岩)的单轴抗压强度之 比来表示的(即图中的s1/sc，反映了岩体中应力与承 载力之间的关系)，其比值小于0。
  15 时属于低地应 力水平，0。15～0。4 为中地应力水平，高于0。4 则是 高地应力水平。岩体的完整性则由岩体RMR 分类 的分值表示，分别以50 和75 为界限划分出3 种完 整性类型的岩体。
  在低地应力水平条件下，不管岩 体的完整程度如何，都不会表现出应力控制型破坏， 视岩体完整程度，围岩分别表现出弹性响应、块体 滑移和散体型破坏。相反地，在高地应力水平的条 件下，应力控制型破坏则为显著特点。
  岩体完整程 度的不同，总体上会影响具体的破坏范围和方式， 但不会从根本上改变破坏的性质。 图 1 所示关于地下工程围岩破坏性质与岩体初 始条件(地应力水平、岩石强度、岩体完整性)之间 的关系，可以很好地帮助判断工程中的岩体潜在稳 定问题的性质和表现形式，从而可在工作开始阶段 制定正确的计算、设计方法。
  这种判断仅是初步的； 地下工程布置、开挖方式等都可以显著影响岩体的 稳定性和破坏方式，这是具体的问题。随着工程开 发进度的推进，我们可以获得更详细的信息，对问 题可以作进一步的判断。
   在工程建设中，特别是大型土木工程建设中， 与矿山工程不同的是，工程场址是经过严格考查和 论证以后选定的，一般不会选在RMR 值小于50 的 破碎岩体中。自然界中RMR 值大于75 的岩体是比 较少见的，在大多数情况下，这类工程场址区岩体 以RMR 值在50～75 者占主导地位，且岩石以中～ 上等强度为主。
  当这些工程在地表以下300 m以内 时，如果不考虑个别因素的影响(如边坡坡脚应力集 中)，最大初始主应力一般不会超过15 MPa，而中 等坚硬岩石的单轴抗压强度可以接近100 MPa，这 种条件属于低地应力的上限状态。
  也就是说，大体 上，对地表边坡工程和浅埋地下工程而言，沿结构 面的块体滑动将是起主导作用的潜在破坏形式，相 关的地质调查、岩体力学试验、岩体稳定性研究和 岩体加固设计都应该围绕这个主题进行。
  显然，一 些建立在连续介质力学理论基础上的分析方法并不 适合于进行这类工程问题的相关研究。 随工程埋深增大，初始应力水平不断增高。工 程开挖后二次应力状态对围岩稳定性的作用增加， 结构面之间的咬合作用也可能随结构面法向应力 (围压)的增高而不断增强。
  在这种条件下，对于RMR 值为50～75 的岩体而言，岩石本身的特性对岩体变 形和稳定的影响可能会将逐渐显示出来，在这就要 求在强调结构面作用的同时，岩石特性不能忽略。 当埋深进一步增大，进入高地应力状态时，应 力起主导作用。
  这种条件下进行一般性的岩石力学 研究时，次一级的结构面如节理的作用可以通过折 减岩体力学参数实现，使得连续介质力学方法能比 较接近地描述岩体的力学行为。但在对特定问题如 岩爆、微震的研究时，结构面的不连续特性仍需特 ? 3446 ? 岩石力学与工程学报 2004 年 图1 地下工程围岩初始条件与围岩稳定状态、破坏特征的一般关系 Fig。
  1 The relationship of initial condition，stability and failure characteristics of rock mass around underground projects 别注意。
  深埋地下工程中的岩爆和大型微震现象的 发生部位多与断层或岩性接触面有关的事实，说明 了结构面在高地应力条件下有时仍可以起重要作 用。 3 非连续介质力学方法与使用条件 前面的分析表明，对地表和浅埋地下工程进行 岩体不连续力学方法的分析研究是非常必要的，但 基于连续介质的力学理论的一些分析方法和计算程 序，在基本理论上，尚不能正确有效地描述岩体的 工程行为。
   目前针对岩体不连续特性的力学方法大致有 5 种，各种方法的特点分述如下[3]： (1) 极限平衡理论方法：主要考查块体由静态 转化为运动状态的力的矢量平衡条件，它不考虑块 体系统内一旦块体发生运动后荷载的重新分布，忽 略块体和结构面的变形。
  典型的例子如石根华等的 块体理论方法和Warburton 的矢量稳定分析方法。 (2) 动量转换方法：该方法认为块体和结构面 均为不可变形的刚体，块体破坏是沿结构面的滑动， 在破坏产生的瞬间，动量在两个接触块体之间发生 转换。
   (3) 不连续变形分析(DDA)方法：该方法认为结 构面是刚体，块体可以是刚体、也可以是变形体。 (4) 模态方法：当把块体处理为刚体时，与下 面要讨论的离散元方法相似。对变形体，模态方法 采用了几何重叠。
  代表性计算机程序为CICE。 (5) 离散元方法：离散元方法的标志性特点包 括两个方面，其一是允许有限位移和离散体的转动 完整岩体块状岩体 破碎岩体 RMR＞75 50＜RMR＜75 RMR＜50 线弹性响应 块体滑动破坏 块体在开挖面的散 体型破坏 临近开挖面的脆性破坏 咬合部位的局部脆性破坏 和沿结构面的散体型破坏 绕开挖面的脆性破坏 完整岩石的脆性破坏和 块体滑移 岩石的挤压与鼓胀，呈弹 性连续体特性 (c) 高地应力水平(s1/sc＞0。
  4) (b) 中地应力水平(0。15＜s1/sc＜0。4) (a) 低地应力水平(s1/sc＜0。15) 咬合部位的局部脆性 破坏和块体滑动 第23 卷第 20 期朱焕春等。
   节理岩体数值计算方法及其应用(一)：方法与讨论? 3447 ? 及脱离；其二是计算过程中可以自动判别块体之间 可能出现的新的接触关系。代表性软件包括早些时 候的TRUBAL，DIBS 和后来出现的3DSHEAR， UDEC 以及近期出现的3DEC 和PFC 系列。
  其中的 PFC 系列，即粒子流方法程序，主要是针对岩体细 观非连续性开发的。 离散元方法认为结构面是可以张开和滑动的， 对块体而言，视问题的需要，既可以处理成可变形 体、也可以处理成不可变形的刚体，结构面本身可 以产生法向和切向变形。
  离散元方法可以对结构面 进行粘合处理，局部或全部地把结构面切割的不连 续岩体蜕化成连续体。因此，离散元的解法是多样 的，对地表和浅部工程，当岩石变形相对结构面变 形可以忽略不计时，可以把块体处理为刚体；当岩 石相对软弱或埋深较大时，可以把块体处理成变形 体；对完整性很好或很差时，特别关心的主要是深 埋条件下的应力水平和状态时，可以蜕化成类似于 连续介质。
   上述 5 种方法所涉及的理论假设以及由这些假 设所规定的基本功能可以归纳如表1。一些高级方 法兼容了连续介质力学的基本特点和功能。 由此可见，在进行工程岩体岩石力学数值分析 之前，正确认识问题的基本性质和潜在破坏方式的 主要类型及其控制因素，是合理选择计算程序和确 定计算方法的重要前提。
  下面就这个问题进行进一 步讨论。 比值(3s1－s3)/sc 也可以用来帮助判断地下工程 开挖中潜在问题的性质，其中的s1，s3 和sc 分别为 垂直于开挖前地下工程纵轴线的横断面上的最大、 最小主应力和岩石单轴抗压强度。
  一般地，该数值 大于0。8 时，岩体开挖后可能出现应力控制型破 坏[4]。 在应力问题不起决定性作用时，岩体不连续面 发育程度可能会起重要作用。一般地，如果地下工 程开挖跨度与岩体中优势结构面间距的比值达到 10∶1 时，就应该考虑结构面可能引起的围岩破坏 问题。
  但对具体工程而言，其破坏形式是否一定会 是结构面控制型为主，还取决于洞室高度和洞室方 位与结构面优势方位、地应力场方位之间的关系等 因素相关。 对于多组结构面的情形，如果已知每一组优势 结构面的平均产状和平均间距，则在任何一个水平 方向上结构面的平均间距可以由下式估计： ?? ??n i i i i d d 1 cos sin 1 q b (1) 式中：di 为第i 组结构面的平均间距， i b 为平均倾 角， i q 为与研究对象考察方向在水平面上的交角。
   现实条件下典型的应力控制型和结构面控制型 破坏可以在现场很好区分。常见的结构面控制型破 坏沿某一个(或几个)特定的结构面发生，破坏后块 体形状完整，不会出现片理化现象。而典型的应力 控制型破坏后，块体因高应力的突然释放而呈碎片 状。
  图2 给出了加拿大一个深埋(埋深达1 000 余 米)矿山工程中发生在岩浆岩中的应力控制型破坏 典型实例，可看到破坏面形态和剥落体形状，显 著不同于一般的块体滑落破坏。 图2 地下工程围岩应力控制型破坏示例 Fig。
  2 Failure example of stress-control-type of rock mass around underground project 强烈的片状化 强烈破碎的碎块 方 法 ? 3448 ? 岩石力学与工程学报 2004 年 4 节理的3DEC 处理 传统的节理网络模拟技术目前在工程上多用来 展示岩体中的节理分布，并未有效地用来作为数值 计算中地质模型的一部分。
  传统的节理网络模拟技 术遵循的是概率统计理论方法，认为节理参数(如倾 向、倾角、间距)服从某种随机分布形式，在根据样 本获得其分布形式的数学表达式以后，再设法(如利 用求解反函数的方式)在计算机上生成服从这种分 布形式的随机数，以模拟节理网络。
  这种方法在工 程中的实用性往往受到很大限制，原因在于： (1) 一些节理参数的分布形式可能非常复杂， 很难寻找到满足精度要求的数学表达式来描述其分 布特征； (2) 对于复杂的分布形式，不能保证利用反函 数法求出的节理参数随机数具有足够的精度，即统 计分布特征可能与样本存在严重误差； (3) 生成三维节理网络依赖每一组节理的倾 向、倾角、间距、和迹长等4 个参数，现实中进行 多组节理的节理网络模拟时，实际上很难用传统方 法获得全部这些参数的统计分布形式并进行反函数 求解。
   为了避免传统方法的上述不足，另一种建立在 离散数学理论基础上的节理网络模拟方法——直接 法可以用来模拟几乎是任何情况下的节理网络。用 直接法生成节理网络的精度，从理论上讲，主要取 决于取样的精度。
  利用直接法生成某一条节理随机 数的思路是[5]： (1) 先对每一个节理参数(如倾向)的样本(地质 编录资料)在其变化区间内进行离散化处理，即分成 若干小区间，求每个小区间内该参数样本的频率(占 样本总体的百分比)。
   (2) 根据需要模拟区域的大小和节理间距的统 计结果计算出在该区域范围内、在每一组节理的法 向方向上需要生成节理的数目。显然，如果这些待 生成的节理随机数在各离散区间的分布频率与样本 的分布频率充分接近，则节理随机数代表了样本的 分布。
  由此可以计算出每个小区间内需要生成的节 理随机数数目，作为目标数。 (3) 利用Mount-Carlo 法生成一个在[0，1]范围 内服从均匀分布的随机数，再把它换算到需要的区 间范围内(与该节理参数统计样本的分布范围一 致)。
  现在的问题是需要验证这个随机数是不是需要 的，视其能否代表样本的分布取舍。 (4) 设目前考察的是第i 条待生成节理的某个 参数的一个随机数，如果把它引入到它所属的小区 间内，计算此时该区间内随机数数目。
  如果此时的 数目不超过目标数，保留该随机数，反之，生成并 考察下一个随机数。 (5) 不断重复上述过程，直到生成全部随机数 为止。 获得了节理随机数只是解决了问题的第一步， 如何把这种方法直接应用到某个计算程序中，在数 值模型中形成可供数值计算用的节理网络模型，则 是问题的另外一个方面。
  3DEC 程序对节理的处理 是按条进行的，它需要指定每一条节理的倾向、倾 角、节理面任一点的坐标(这3 项必选)和节理的连 通率(选择项)。其中的连通率是相对于节理切割块 体数目的一个相对概念，有别于现实中关于节理连 通率的定义，该选择项有助于模拟节理迹长。
  3DEC 程序同时提供了在[0，1]区间内服从正态分布和均 匀分布的两个随机“籽”。显然，直接法和3DEC 模 拟节理的技术思路之间有着很好的结合点。直接法 可以直接为3DEC 提供它需要的节理倾向和倾角随 机数，而遵照测线法(scanline)进行节理采样的原理 [1]，在模型内规定一条(或多条)不与节理组平行的测 线和获得在测线上的一个随机起始点以后，就可以 在模型中生成第1 条节理；类似地，利用直接法可 以生成下一条节理的倾向、倾角及其与这条节理的 间距，并计算出下一条节理面上一个点的坐标，从 而生成下一条节理，直至生成节理网络。
   按以上方法在 3DDEC 模型中生成的节理没有 考虑节理的长度，因此它可能贯穿整个模型而过度 估计了节理连通性。这个问题可以通过下述方法得 到修正：首先，根据现场获得的各节理组中长大贯 通型节理的长度和间距等，利用3DEC 的隐藏(HIDE) 功能，按区域生成长大节理；其次，根据已经形成 的模型块体尺寸和一般节理的长度，选择合适的连 通率参数值，实现对迹长的模拟。
   3DEC 在模型中生成节理需要的一些关键参数 如倾向、倾角和位置，都对应于节理地质编录过程 中的一个采样参数。这就使得有可能把节理编录资 料直接移植到模型之中，此时的节理模型，从性质 上就不再是一个“不确定性的随机模型”，而是一个 “确定性模型”，其理论意义和工程意义都是突破 性的，在理论认识上突破了对节理网络只能进行统 计处理的限制。
  体现在工程实践中的不同是：确定 第23 卷第 20 期朱焕春等。 节理岩体数值计算方法及其应用(一)：方法与讨论? 3449 ? 性模型除了能够回答问题的性质和统计量特征等以 外，还可以明确指出问题出现在什么位置。
   由于 3DEC 模型中的每一条节理都有一个特定 的序号，结合3DEC 的其他功能，使得对每一条节 理赋以不同的力学参数值成为可能。类似地，这些 力学参数值可以服从任意给定的数学分布形式，也 可以遵循离散分布特征，甚至可以为每一条节理明 确每一个参数值。
   5 讨论 具备良好工程价值的数值计算应该遵循如下几 条原则： 首先，是作必要的工程调研和查勘，理解岩体 实际条件和存在的各种现象，对问题的潜在性质有 一个正确的估计，并以此为依据制定研究方案、选 择正确的数值计算方法和程序。
  切忌只满足于向研 究委托人索取几个数据，用自己手上仅有的程序获 得计算结果。同时，从工程建设的严肃性出发，一 般不应该将有过重大改动的程序直接应用到工程问 题中。 数值计算只是一种手段，虽然数值计算技术的 发展很大程度上提高了我们认识复杂问题的可能 性，但就像任何科学研究或医疗诊断中的先进手段 一样，数值计算结果给我们揭示的是给定条件下的 某些特定指标值(如应力、位移等)。
  这些指标值到 底说明了什么样的工程问题、计算给定条件和现实 相差会有多大及其可能造成多大的影响，如医生依 据检验结果进行诊断一样，需要全面的知识和必须 的经验才能正确回答。在获得确诊结果之前，当然 需要一个初步的判断，再考虑进行什么样的检验。
   其次，岩体工程建设是一项系统性的实践活动， 强调相关研究人员的知识面和工程经验是必要的。 数值计算人员也是如此，需要具备数值计算全过程 涉及的知识和经验，概括来说，包括数值计算的原 理、岩体工程地质、岩体力学和工程设计意图等。
   对于正在实施的工程，根据已经揭示出来的工程现 象对计算条件进行考证和修正，是用数值计算结果 预测未来问题的一个非常有效的手段；对于尚未实 施的工程，则更需要丰富的经验。 最后，需要根据工程要求和所关心的问题，在 数值计算中采用必要的补充手段，明确回答问题。
   数值计算结果一般只给出基本参数值(如应力和位 移等)。很多情况下，这些基本参数值的绝对大小并 不能明确回答工程中关心的问题。比如，多高的应 力水平才会诱发影响工程运转的问题。
  要回答这样 一个问题，只依赖理论依据是不行的，经验积累和 到工程中去寻找现象、通过建立工程现象与数值计 算结果之间的关系，是回答这样一个问题的有效途 径。收起