岩石力学与工程学报 23(20)：3444～3449
2004年10月 Chinese Journal of Rock Mechanics and Engineering Oct . ，2004
第23卷 第20期
节理岩体数值计算方法及其应用(一) ：方法与讨论
朱焕春1，2 Brummer Richard 2 Andrieux Patrick 2
(1武汉大学水利水电学院 武汉 430072) (2Itasca(加拿大) 咨询公司 索德柏立 加拿大)
摘要 叙述了工程岩体的典型破坏类型、发生条件和判断方法，简单介绍了针对节理岩体的数值计算方法发展过程、现状和特点，重点叙述了利用3DEC 基本功能生成三维随机节理网络和节理力学参数的随机赋值方法。其中三维节理网络技术可以非常逼真地模拟现实条件下各种复杂形式的节理分布，甚至可以直接利用节理的地质编录资料在数值模型中生成确定位置上的确定性节理，在某种程度上一体化地实现了节理岩体的计算机模拟和数值计算，促进了不连续岩体、岩石计算方法在应用技术方面的发展。最后讨论了不连续力学数值计算技术的工程应用，指出了合理运用数值计算方法解决工程问题的一般原则。
关键词 岩石力学，数值模拟，不连续体，离散元，随机，节理岩体
分类号 TU 12 文献标识码 A 文章编号 1000-6915(2004)20-3444-06
NUMERICAL METHODS AND APPLICATION FOR JOINTED ROCK
MASS ，PART 1：APPROACHES AND DISCUSSIONS
Zhu Huanchun12，Brummer Richard2，Andrieux Patrick2
，
(1School of Water Resources and Hydropower Engineering，Wuhan University， Wuhan 430072 China )
(2Itasca Consulting Canada，Inc . ，Sudbury ，ON P 3E 1G 1， Canada )
Abstract This paper narrates the typical modes of rock mass failures in engineering activities，the conditions under which such failures occur，and the approaches to predict these failures. A brief description is presented of the evolution ，the current state，and the feature of numerical modelling methods for jointed rock mass. Attention is paid to the technique of generating 3-dimensional joint net in rock mass and randomly simulating the properties of each individual joint on the basis of accessible functions in a 3DEC numerical model. Such technique enables 3DEC to build up a model that contains various joints corresponding to the realistic distributions of joints in rock mass and，furthermore ，achieve a determinative description of joint distributions by transplanting joint mapping data into a 3DEC model. This，in some degree，fulfills the numerical modelling of geomechanical behaviours of jointed rock mass combined with the computerized simulations of its geometry. The essential principles are discussed of applying numerical modelling techniques for discontinuum into engineering projects. Key words rock mechanics，numerical modelling，discontinuum ，discrete element，random ，jointed rock mass
1 概 述
数值计算已经被普遍应用于工程设计中解决各种岩石力学问题。随着数值计算技术的发展，它在很大程度上已经取代了一些传统的研究手段，使得
2003年5月30日收到初稿，2003年7月4日收到修改稿。
岩土工程中的许多岩石力学问题有可能较快地得到解决。目前已经开发出来的岩土工程数值计算软件及相关领域的工程软件，已为岩土工程设计和施工中相关问题的研究提供了充分的选择余地。在目前条件下如何正确选择和有效地使用这些软件，则成为岩石力学数值计算实际应用中的一个现实问题。
作者 朱焕春 简介：男，1963年生，博士，1985年毕业于河海大学水文地质与工程地质专业，现任武汉大学教授、ITASCA 高级岩石力学专家，主要从事岩石力学及相关方向的技术咨询与研究。E-mail ：[email protected]。
第23卷 第20期 朱焕春等. 节理岩体数值计算方法及其应用(一) ：方法与讨论 • 3445 •
对岩体几何特征的研究开始于20世纪70年代，主要是采用概率统计方法研究节理的产状、间距和迹长等几何参数的统计分布[1]以及岩体节理网络计算机模拟等，其最终目的还是为研究节理岩体的力学特性服务。然而，迄今为止，能正确地描述岩体中的节理分布，并把它们与岩石力学数值计算方法相结合的应用实例还很鲜见。
本文的研究，从岩体破坏类型及其发生条件入手，首先界定结构面控制型破坏的发生条件，从而规定了需要采用不连续力学方法进行岩石力学分析方法的一般范围，利用3DEC 程序的基本功能直接在3DEC 计算模型中模拟节理网络方法和节理力学指标值的随机赋值方法，针对中国水电工程特点，讨论了进行数值计算时需要考虑的有关岩体非连续性问题。
视岩体完整程度，围岩分别表现出弹性响应、块体滑移和散体型破坏。相反地，在高地应力水平的条件下，应力控制型破坏则为显著特点。岩体完整程度的不同，总体上会影响具体的破坏范围和方式，但不会从根本上改变破坏的性质。
图1所示关于地下工程围岩破坏性质与岩体初始条件(地应力水平、岩石强度、岩体完整性) 之间的关系，可以很好地帮助判断工程中的岩体潜在稳定问题的性质和表现形式，从而可在工作开始阶段制定正确的计算、设计方法。这种判断仅是初步的；地下工程布置、开挖方式等都可以显著影响岩体的稳定性和破坏方式，这是具体的问题。随着工程开发进度的推进，我们可以获得更详细的信息，对问题可以作进一步的判断。
在工程建设中，特别是大型土木工程建设中，与矿山工程不同的是，工程场址是经过严格考查和论证以后选定的，一般不会选在RMR 值小于50的破碎岩体中。自然界中RMR 值大于75的岩体是比较少见的，在大多数情况下，这类工程场址区岩体以RMR 值在50～75者占主导地位，且岩石以中～上等强度为主。当这些工程在地表以下300 m以内时，如果不考虑个别因素的影响(如边坡坡脚应力集中) ，最大初始主应力一般不会超过15 MPa，而中等坚硬岩石的单轴抗压强度可以接近100 MPa，这种条件属于低地应力的上限状态。也就是说，大体上，对地表边坡工程和浅埋地下工程而言，沿结构面的块体滑动将是起主导作用的潜在破坏形式，相关的地质调查、岩体力学试验、岩体稳定性研究和岩体加固设计都应该围绕这个主题进行。显然，一些建立在连续介质力学理论基础上的分析方法并不适合于进行这类工程问题的相关研究。
随工程埋深增大，初始应力水平不断增高。工程开挖后二次应力状态对围岩稳定性的作用增加，结构面之间的咬合作用也可能随结构面法向应力(围压) 的增高而不断增强。在这种条件下，对于RMR 值为50～75的岩体而言，岩石本身的特性对岩体变形和稳定的影响可能会将逐渐显示出来，在这就要求在强调结构面作用的同时，岩石特性不能忽略。
当埋深进一步增大，进入高地应力状态时，应力起主导作用。这种条件下进行一般性的岩石力学研究时，次一级的结构面如节理的作用可以通过折减岩体力学参数实现，使得连续介质力学方法能比较接近地描述岩体的力学行为。但在对特定问题如岩爆、微震的研究时，结构面的不连续特性仍需特
2 岩体破坏类型与发生条件
地表和近地表岩体工程如人工边坡、深度不大的矿坑工程涉及的岩体初始应力水平一般都不高，岩体的潜在破坏方式主要取决于岩体本身特性。对于坚硬和中等坚硬岩石地区的这类工程，一般而言，结构面对岩体的潜在破坏方式起决定性作用。
地下工程围岩破坏可以概括地分为两大类：一类是应力控制型破坏；另一类是块体控制型破坏(也称结构面控制型破坏或重力型破坏) 。这两类破坏的岩体中各有其发生条件，前者主要发生在初始地应力水平相对较高(相对于岩体承载能力) 的岩体中，后者则多出现在地表和地下浅埋岩体工程中。
修建在以硬岩为主的岩体地下工程，其围岩的初始应力水平在很大程度上决定了地下工程围岩二次应力的量级，而岩体的承载能力取决于岩石强度和结构面发育程度。这3者之间的关系直接影响着围岩的状态和潜在问题的性质[2]。图1描述了这3个因素与围岩破坏特征之间的一般关系。其中的左侧第1行描述了初始地应力水平，它是用初始地应力中的最大主应力与岩石(硬岩) 的单轴抗压强度之比来表示的(即图中的σ1/σc ，反映了岩体中应力与承载力之间的关系) ，其比值小于0.15时属于低地应力水平，0.15～0.4为中地应力水平，高于0.4则是高地应力水平。岩体的完整性则由岩体RMR 分类的分值表示，分别以50和75为界限划分出3种完整性类型的岩体。在低地应力水平条件下，不管岩体的完整程度如何，都不会表现出应力控制型破坏，
• 3446 • 岩石力学与工程学报 2004年
完整岩体
块状岩体
破碎岩体
图1 地下工程围岩初始条件与围岩稳定状态、破坏特征的一般关系
Fig.1 The relationship of initial condition，stability and failure characteristics of rock mass around underground projects
别注意。深埋地下工程中的岩爆和大型微震现象的发生部位多与断层或岩性接触面有关的事实，说明了结构面在高地应力条件下有时仍可以起重要作用。
转化为运动状态的力的矢量平衡条件，它不考虑块体系统内一旦块体发生运动后荷载的重新分布，忽略块体和结构面的变形。典型的例子如石根华等的块体理论方法和Warburton 的矢量稳定分析方法。 (2) 动量转换方法：该方法认为块体和结构面均为不可变形的刚体，块体破坏是沿结构面的滑动，在破坏产生的瞬间，动量在两个接触块体之间发生转换。
(3) 不连续变形分析(DDA)方法：该方法认为结构面是刚体，块体可以是刚体、也可以是变形体。 (4) 模态方法：当把块体处理为刚体时，与下面要讨论的离散元方法相似。对变形体，模态方法采用了几何重叠。代表性计算机程序为CICE 。
(5) 离散元方法：离散元方法的标志性特点包括两个方面，其一是允许有限位移和离散体的转动
3 非连续介质力学方法与使用条件
前面的分析表明，对地表和浅埋地下工程进行岩体不连续力学方法的分析研究是非常必要的，但基于连续介质的力学理论的一些分析方法和计算程序，在基本理论上，尚不能正确有效地描述岩体的工程行为。
目前针对岩体不连续特性的力学方法大致有5种，各种方法的特点分述如下[3]：
(1) 极限平衡理论方法：主要考查块体由静态
第23卷 第20期 朱焕春等. 节理岩体数值计算方法及其应用(一) ：方法与讨论 • 3447 •
及脱离；其二是计算过程中可以自动判别块体之间可能出现的新的接触关系。代表性软件包括早些时候的TRUBAL ，DIBS 和后来出现的3DSHEAR ，UDEC 以及近期出现的3DEC 和PFC 系列。其中的PFC 系列，即粒子流方法程序，主要是针对岩体细观非连续性开发的。
离散元方法认为结构面是可以张开和滑动的，对块体而言，视问题的需要，既可以处理成可变形体、也可以处理成不可变形的刚体，结构面本身可以产生法向和切向变形。离散元方法可以对结构面进行粘合处理，局部或全部地把结构面切割的不连续岩体蜕化成连续体。因此，离散元的解法是多样的，对地表和浅部工程，当岩石变形相对结构面变形可以忽略不计时，可以把块体处理为刚体；当岩石相对软弱或埋深较大时，可以把块体处理成变形体；对完整性很好或很差时，特别关心的主要是深埋条件下的应力水平和状态时，可以蜕化成类似于连续介质。
上述5种方法所涉及的理论假设以及由这些假设所规定的基本功能可以归纳如表1。一些高级方法兼容了连续介质力学的基本特点和功能。
表1 各不连续介质的力学方法的基本特点一览表 Table The basic characteristics of discontinuous
mechanical methods
接触块体块体材料岩体力学力-
变形 应变 破裂
方性质 性质 数量 性质 类型 状态 位移 法力 可可小大小大非无有松紧静动刚刚大线仅和变变变变应应少 线破破散密力力性 性 量 性 力 位形 形 形 形 变 变 性 裂 裂 体 型 学 学
移
1 √ √ √ √ O √ √ √ 2 √ √
√
√ √
√ O √ √
√ √ O √
最小主应力和岩石单轴抗压强度。一般地，该数值大于0.8时，岩体开挖后可能出现应力控制型破坏[4]。
在应力问题不起决定性作用时，岩体不连续面发育程度可能会起重要作用。一般地，如果地下工程开挖跨度与岩体中优势结构面间距的比值达到10∶1时，就应该考虑结构面可能引起的围岩破坏问题。但对具体工程而言，其破坏形式是否一定会是结构面控制型为主，还取决于洞室高度和洞室方位与结构面优势方位、地应力场方位之间的关系等因素相关。
对于多组结构面的情形，如果已知每一组优势结构面的平均产状和平均间距，则在任何一个水平方向上结构面的平均间距可以由下式估计：
1
d =n (1)
cos θi sin βi ∑d i i =1式中：d i 为第i 组结构面的平均间距，βi 为平均倾角，θi 为与研究对象考察方向在水平面上的交角。
现实条件下典型的应力控制型和结构面控制型破坏可以在现场很好区分。常见的结构面控制型破坏沿某一个(或几个) 特定的结构面发生，破坏后块体形状完整，不会出现片理化现象。而典型的应力控制型破坏后，块体因高应力的突然释放而呈碎片状。图2给出了加拿大一个深埋(埋深达1 000 余米) 矿山工程中发生在岩浆岩中的应力控制型破坏典型实例，可看到破坏面形态和剥落体形状，显著不同于一般的块体滑落破坏。
3 √ √ √ √ √ √ √ √ O √ √
4 √ √ √ √ √ √ √ √ √ √ O √ √ √ O √ √ √ 5 √ √ √ √ √ √ √ √ √ √ √ √ O √ √ √ √ √ 注：“√”表示该具备该功能或适用于该条件；“O ”表示具备该项功能但可能不充分有效；空白表示不具备该项功能或没有进行这方面的考虑。
强烈的片状化
由此可见，在进行工程岩体岩石力学数值分析之前，正确认识问题的基本性质和潜在破坏方式的主要类型及其控制因素，是合理选择计算程序和确定计算方法的重要前提。下面就这个问题进行进一步讨论。
比值(3σ1－σ3)/σc 也可以用来帮助判断地下工程开挖中潜在问题的性质，其中的σ1，σ3和σc 分别为
强烈破碎的碎块
图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 的其他功能，使得对每一条节理赋以不同的力学参数值成为可能。类似地，这些力学参数值可以服从任意给定的数学分布形式，也可以遵循离散分布特征，甚至可以为每一条节理明确每一个参数值。
涉及的知识和经验，概括来说，包括数值计算的原理、岩体工程地质、岩体力学和工程设计意图等。对于正在实施的工程，根据已经揭示出来的工程现象对计算条件进行考证和修正，是用数值计算结果预测未来问题的一个非常有效的手段；对于尚未实施的工程，则更需要丰富的经验。
最后，需要根据工程要求和所关心的问题，在数值计算中采用必要的补充手段，明确回答问题。数值计算结果一般只给出基本参数值(如应力和位移等) 。很多情况下，这些基本参数值的绝对大小并不能明确回答工程中关心的问题。比如，多高的应力水平才会诱发影响工程运转的问题。要回答这样一个问题，只依赖理论依据是不行的，经验积累和到工程中去寻找现象、通过建立工程现象与数值计算结果之间的关系，是回答这样一个问题的有效途径。本文的续篇“节理岩体数值计算方法及其应用(二) ：工程应用”将结合实例详细讨论这类问题。
参 考 文 献
1
Pirest S D，Hudson J A. Estimation of discontinuity spacing trace length using scanline[J]. Int. J. Rock Mech. Min. Sci.，1981，19：183～191 2
Martin C D，Kaiser P K，McCreath D R. Hoek-Brown parameters for predicting the depth of brittle failure around tunnels[J]. Canadian Geotechnical J.，1999，36(1)：136～151 3
Kaiser P K，Martin C D，Sharp J，et al. Underground works in hard rock tunnelling and mining[R]. Melbourne，Australia ：Technomic Publishing Co.，2000，1：841～926 4
Cundall P A，Hart R D. Numerical modeling of discontinuum[A]. In：Hudson J A ed. Comprehensive Rock Engineering[C]. Oxford：Pergamon Press Ltd.，1993，231～243 5
朱焕春，陈文理. 生成节理随机数的直接法[J]. 武汉水利电力大学学报，1997，30(2)：10～12
5 讨 论
具备良好工程价值的数值计算应该遵循如下几条原则：
首先，是作必要的工程调研和查勘，理解岩体实际条件和存在的各种现象，对问题的潜在性质有一个正确的估计，并以此为依据制定研究方案、选择正确的数值计算方法和程序。切忌只满足于向研究委托人索取几个数据，用自己手上仅有的程序获得计算结果。同时，从工程建设的严肃性出发，一般不应该将有过重大改动的程序直接应用到工程问题中。
数值计算只是一种手段，虽然数值计算技术的发展很大程度上提高了我们认识复杂问题的可能性，但就像任何科学研究或医疗诊断中的先进手段一样，数值计算结果给我们揭示的是给定条件下的某些特定指标值(如应力、位移等) 。这些指标值到底说明了什么样的工程问题、计算给定条件和现实相差会有多大及其可能造成多大的影响，如医生依据检验结果进行诊断一样，需要全面的知识和必须的经验才能正确回答。在获得确诊结果之前，当然需要一个初步的判断，再考虑进行什么样的检验。
其次，岩体工程建设是一项系统性的实践活动，强调相关研究人员的知识面和工程经验是必要的。数值计算人员也是如此，需要具备数值计算全过程