等效岩体三维随机节理网络模型构建及其在两河口水电站中的应用

武 娜，梁正召，宋文成，李万润

大连理工大学土木工程学院，大连 116024

在岩土与地质工程领域，岩体是由大小不同、形态各异的岩块及结构面组成的地质体.其中，节理是一种最常见的结构面.岩体中节理的形态、尺寸、空间分布及组合特征决定了其力学性质的不连续、非均匀性和各向异性，同时也引起了岩体力学参数的尺寸效应[1].由于岩体结构的复杂性，如何准确获取节理几何特征参数，建立反映节理真实空间状态的等效岩体三维随机节理网络模型，进而确定合理的等效力学参数是近年来岩石力学研究的一个热点.

由于节理几何参数的分布具有一定的统计特征和概率分布特性，目前常采用统计学方法对其进行定量化描述和分析.传统的节理采样方法主要分为两类：（1）基于岩体露头统计的测线法和窗口法[2]；（2）基于钻孔的岩芯统计法[3].其中测线法和窗口法应用最为广泛，而岩芯统计法往往用于节理网络模型的校核.但由于节理数量庞大，人工罗盘测量效率低，并且测量精度不能得到保证，传统的方法逐渐被全站仪法[4−5]、三维激光扫描法[6]和摄影测量法[7−8]等先进的节理测量方法所替代，并成为主流手段.其中，非接触岩体数字测量技术ShapeMetriX3D（3GSM）在边坡岩体节理几何信息调查中得到广泛应用[9−10].例如，曹俊等[9]借助ShapeMetrix3D测量系统，对甲玛矿五个分段的巷道节理面进行测量统计，得到节理面几何特征参数的数据，获得节理面的概率分布类型和特征值.

三维随机离散节理网络（Discrete fracture network，DFN）模型是建立等效岩体三维随机节理网络模型的基础.受地质活动的影响，节理通常随机地，无序地分布于岩体中.因此，根据每条节理的几何形态建立真实的三维随机DFN模型很难实现.目前，大多以统计学与概率论为基础，假定节理面为等概率出现的光滑平直薄圆盘，并运用Monte-Carlo随机模拟理论，构建与实际节理分布具有相似统计的等效三维随机DFN模型.近年来，岩体三维随机DFN模型的研究在国内外取得了重要成果，并得到广泛应用[11−13].吴顺川等[14]基于Monte-Carlo方法，采用Matlab软件，并通过OPEN GL三维可视化技术构建了三维DFN模型，模拟结果与实测数据具有良好的相似统计性.

随着三维随机DFN技术和离散元方法的发展，等效岩体技术或被称为合成岩体技术逐渐应用于岩体工程中.但是，目前能建立三维等效岩体并应用于岩体力学性质分析的软件是有限的.常用的方法主要包括三种：（1）基于有限元与离散元相结合的方法，采用ELFEN软件建立[15]；（2）基于颗粒流理论，利用 PFC3D软件建立[16−17]；（3）基于块体离散元理论，利用3DEC软件建立[18−20].例如，杨忠民等[21]采用3DEC离散元软件，研究了三维节理网络模型中节理几何参数的特征值和分布类型对岩体表征单元和单轴抗压强度的影响规律，分析了岩体强度对各参数的敏感性.但实际岩体中常常含有大量的非贯通节理，采用上述模型计算时，常常导致计算速度大幅度降低甚至出现计算机内存不足的情况.

因此，本文采用Baecher模型和Monte-Carlo方法，基于有限元软件RFPA3D实现了三维随机DFN模型的重构，并给出了生成三维随机DFN模型的计算模块，实现了三维节理岩体尺寸效应的模拟研究.然后，以两河口水电站左岸边坡坝址区下游处节理岩体为研究对象，建立与实际节理分布等效的三维随机DFN模型，并验证本文所构建模型的合理性.最后，基于等效岩体三维随机节理网络模型，开展了三维节理岩体力学参数的尺寸效应研究.

1 三维随机离散节理网络模型构建

在三维随机DFN模型中，通常采用Baecher模型描述节理面的几何特性，其基本特点是假定节理面形状为平面圆盘形（或椭圆形），其空间特性由其中心点位置、产状和直径来定义.其中，节理中心点位置在研究区域内服从互相独立且具有均匀概率的泊松分布；节理面产状由节理面倾向和倾角定义；节理直径和体密度可根据节理迹长和线密度分别求解[22].本文采用C语言和Fortran语言，基于Baecher模型和Monte-Carlo方法，利用RFPA3D软件实现了三维随机DFN模型的重构，并给出了生成三维随机DFN模型的计算模块，具有对网络模拟所需基本参数输入、编辑和修改的功能，实现了对输入的基本参数自动转化为图形和数据的功能.其中，三维随机DFN模型构建的计算流程主要参考吴顺川等[14, 23 ]学者的论文.

图1为利用RFPA3D内置的DFN模块，生成研究区域为 13 m × 13 m × 13 m 的三维随机 DFN 模型的透视图.在透视图中，三种不同颜色圆盘代表不同的节理组.相对于其它DFN模型，该DFN模块特有的特征如下：（1）显示任何一组或几组节理面的立体图和透视图，使其空间分布状态以图形方式展示；（2）对图形可进行任意角度的旋转、移动、放大和缩小；（3）显示节理的交切面，并沿着平行于坐标轴方向对模型进行剖切，便于查看模型结果.除此以为，该模型还可以对不同节理组选用不同的颜色表示，使图形显示更为鲜明.

图1 节理的透视化模型.（a）一组节理面；（b）两组节理面；（c）三组节理面Fig.1 Perspective model of fracture: (a) a set of joint planes; (b) two sets of joint planes; (c) three sets of joint planes

2 等效岩体三维随机节理网络模型构建

首先，采用RFPA3D软件建立大型黏结颗粒有限元模型，然后将三维随机DFN模型嵌入到有限元模型中，基于八叉树结构快速赋予节理面单元不同的物理力学参数，便可建立工程尺度的等效岩体模型，如图2所示.为了描述岩石和节理的这种非均匀性，RFPA3D软件中假定节理和岩体细观单元力学性质满足Weibull分布，并采用弹性损伤本构模型，可通过对等效岩体三维随机节理网络模型进行单轴、三轴加载等数值模拟，定量研究不同工程尺度条件下节理岩体尺寸效应与各向异性等力学性质.更重要的是，该方法可以模拟节理岩体中裂纹的起裂、扩展、贯通和失稳过程.由于篇幅有限，其计算原理的详细介绍可参考文献[24]～[25].

图2 等效岩体三维随机节理网络模型的构建.（a）岩石模型；（b）DFN 模型；（c）等效岩体模型Fig.2 Construction of three-dimensional random DFN model of equivalent rock mass: (a) rock model; (b) DFN model; (c) equivalent rock mass model

3 岩体工程实例分析

节理几何特性参数的获取是建立三维DFN的前提和基础.本文以两河口水电站左岸边坡坝址区下游处岩体为背景（如图3（a）所示），采用3GSM软件，分析研究区域节理的几何参数，为边坡岩体力学性质分析提供参考资料.两河口水电站位于四川省甘孜州雅江县境内的雅砻江干流上，坝址区位于雅砻江与鲜水河汇合口（即两河口）下游约 0.6～3.6 km河段上.从区域构造环境看，两河口库区地质构造较简单，断层少且规模小，不存在由于强烈构造作用而形成大范围岩体破碎的岸坡.工程边坡岩体中含有大量的节理、裂隙等交错不连续面，影响岩体的力学行为.

3.1 两河口边坡岩体节理面识别

按要求对研究区域岩体拍照，把获得的左、右视图导3GSM软件，圈定出重点测量区域，3GSM系统自动根据像素点合成对应的三维图形.再通过标杆以及地质罗盘对图形的尺寸和方向进行修正，最终实现节理面几何信息统计，其中节理的上半球赤平投影图如图3（b）所示.

图3 基于 3GSM 软件节理分布图.（a）坝址区下游处边坡岩体；（b）上半球赤平投影图Fig.3 Fracture distribution diagrams based on 3GSM: (a) slope rock mass downstream of dam site; (b) stereogram of the upper hemisphere

3.2 节理优势组划分

岩体中不同节理面的形成时期和原因不尽相同，导致节理面也存在差异.但岩体中节理面的发育又具有一定的规律性和方向性.岩体中通常包含多组节理面，而节理面的统计分析是针对某一分组而言.因此，需根据产状进行节理面优势组划分.本文在合成三维立体图形的基础上，根据节理的倾向和倾角并借助上半球赤平投影图对节理面进行优势组划分，研究区域内节理可分为3组，如图3(b)所示.另外，3GSM软件可直接输出每组节理的倾向、倾角、迹长、间距、线密度的均值以及每组节理中每条节理面信息的几何参数（倾向、倾角、迹长、间距）.

中国水电顾问集团成都勘测设计研究院通过对两河口水电站边坡坝区13287条节理统计分析表明：坝区优势面主要有5组（占统计总数的70%），其中3组节理面的倾向和倾角如下：235～270∠5～30、330～360∠5～35、90～120∠70～85.而借助3GSM软件得到两河口下游处，1#、2#和3#节理面的倾向和倾角分别为 246.15∠25.98、351.34∠28.16和94.73∠83.8.通过对比发现，岩体优势节理面的倾向和倾角与工程地质勘察结果比较一致，例如 246.15∠25.98与 235～270∠5～30、351.34∠28.16与 330～360∠5～35、94.73∠83.8与90～120∠70～85较一致.从而验证了研究区内数据统计的准确性.

3.3 节理概率模型确定

由于岩体内每组节理面具有相似统计性，因此，节理面的几何参数可通过相应的概率分布类型来描述.但其相对的分布类型无法直接从3GSM软件中输出.因此，根据3GSM软件中每条节理面几何参数分布信息，利用origin软件，分别建立频率分布直方图，再对其数据进行非线性回归分析，可得到相应的概率分布类型.以1#节理面为例，通过分析节理面倾向、倾角、迹长和间距的概率分布直方图，可得到节理面概率分布类型，相关系数（R2）及相应的均值（E(X)）和方差（D(X)），如图4示.结果表明，节理面的倾向和倾角均服从正态分布，迹长服从对数正态分布，而间距服从负指数分布.不同概率分布类型的相关系数均大于0.9，可认为结果拟合较好.因此，利用origin软件对数据进行非线性回归分析获得节理几何参数的概率分布和方差是合理的.最后，可分别得到三组不同节理面几何参数概率分布的类型、均值和方差，如表1所示.

表1 节理面几何参数的概率分布及数值Table 1 Probability distribution and value of geometric parameters of fractures

图4 1#节理面几何参数概率分布直方图.（a）节理面倾向；（b）节理面倾角；（c）节理面迹长；（d）节理面间距Fig.4 Probability distribution histogram of geometric parameters of 1# fracture surfaces: (a) joint plane dip; (b) joint plane inclination; (c) trace length of joint surface; (d) joint plane spacing

4 三维随机 DFN 模型校核

基于Monte-Carlo方法所得到的三维随机DFN模型并不是真实的节理网络模型，而是对真实节理图形的一个随机抽样，只有数学意义上的相似统计性.所以需对本文建立的节理网络模型进行校核.校核方法主要有图形对比法和参数对比法[26].图形对比法只能依靠肉眼来判别，无定量指标，可靠性较差.本文采用数据对比法验证三维随机DFN模型的准确性.首先，对三维随机DFN模型中平行于Y轴的多个平面进行剖割，获取相应的剖面，如图5（a）所示.然后，在所获得的剖面上，采用类似测线法的原理获得不同剖面内节理面出露迹线的几何参数，并将结果与现场资料进行对比，如图5（b）所示.表2给出了根据三维随机DFN模型获得模拟值与现场实测值的对比.评价模拟效果精度常采用相对误差来表示，若相对误差小于30%，即可认为生成的三维DFN模型效果良好[26].

表2 等效岩体三维节理网络模拟数据检验Table 2 Verification of three-dimensional DFN model data of equivalent rock mass

图5 等效岩体三维随机 DFN 模型剖面截取.（a）节理面剖面；（b）节理面剖面出露迹线Fig.5 Profile capture of three-dimensional random DFN model of equivalent rock mass: (a) joint surface profile; (b) exposed trace of joint surface profile

通过对比三维网络模型的模拟值和实测值可以发现，三组节理面几何参数的相对误差均小于30%.其中，相对误差的最大值和最小值分别为11.74%和2.84%.可见，本文构建的岩体三维DFN模型可应用于三维节理岩体力学性质和实际工程问题的研究中.

5 三维节理岩体尺寸效应

本文以两河口水电站左岸边坡坝址区下游处边坡岩体为例，基于构建的等效岩体三维随机节理网络模型，进行三维节理岩体尺寸效应的研究.首先，根据表1，生成模型尺寸为 15 m × 15 m × 15 m三维随机DFN模型，获取节理中心点坐标，法向量和半径等信息；而后，以三维随机DFN模型的中心为中心分别建立模型边长为 1，2，···，12 m 的共12种不同尺寸的等效岩体三维随机节理网络模型，模型边边长用L表示，如图6所示；最后，基于RFPA3D软件，开展三维节理岩体尺寸效应的研究.模型中节理与岩石的力学参数见表3.在计算过程中，模型底面（YOZ面）中心点在水平和竖直方向分别固定，其余点竖直方向固定，水平方向自由.在模型的上表面（与X轴垂直的上表面）施加位移荷载直至模拟失稳破坏，施加的位移荷载的大小为 2×10−5倍的模型边长.

表3 研究区域内岩石和节理的参数Table 3 Parameters of rocks and fractures in the study area

图6 三维节理岩体尺寸效应研究.（a）L=1 m；（b）L=2 m；······；（f）L=12 mFig.6 Scale effect of three-dimensional fractured rock mass: (a) model side length is 1 m; (b) model side length is 2 m; ······; (f) model side length is 12 m

图7（a）为三维节理岩体单轴抗压强度和弹性模量随模型尺寸的变化规律图.结果表明，三维节理岩体单轴抗压强度和弹性模量均随着模型尺寸的增加先减小然后趋于稳定，表现出明显的尺寸效应.节理岩体表征单元体的尺寸可以根据变化系数的大小确定[27−28].单轴抗压强度和弹性模量变化系数随模型尺寸的变化规律如图7（b）所示.取变化系数小于等于5%时模型的最小尺寸为节理岩体的REV，那么节理岩体单轴抗压强度和弹性模量的 REV 尺寸分别为 8 m × 8 m × 8 m和5 m ×5 m × 5 m.岩体的 REV 尺寸为不同参数 REV 的最大值，因此，综合考虑岩体的单轴抗压强度和弹性模量参数，岩体的 REV 可评估为 8 m × 8 m ×8 m.当岩体的 REV 取 8 m × 8 m × 8 m 时，岩体的等效单轴抗压强度和弹性模量分别为25.159 MPa和19.443 GPa.关于该模拟方法与实际工程问题的有效性验证，接下来作者会进行深入分析和研究.

图7 数值模拟结果.（a）岩体力学参数尺寸效应；（b）岩体力学参数变化系数随模型变化规律Fig.7 Numerical simulation results: (a) scale effect of rock mass mechanical parameters; (b) variation law of rock mass mechanical parameter variation coefficient with model size

6 结论

为研究节理岩体的力学性质，基于有限元软件分别建立三维随机DFN模型和等效岩体三维随机节理网络模型，并以两河口水电站左岸边坡坝址区下游处边坡岩体为例进行工程应用，其主要结论如下：

（1）基于 Baecher模型和Monte-Carlo方法，利用RFPA3D软件，实现了三维随机DFN的重构，并给出生成三维随机DFN模型的计算模块.该方法具有对网络模拟所需基本参数输入、编辑和修改等功能，现实了对输入基本参数自动转化为图形和数据的功能.

（2）考虑岩石和节理的非均匀性，采用弹性损伤本构模型，结合嵌入的随机DFN模型，利用RFPA3D软件建立了工程尺度等效岩体三维随机节理网络模型，并实现了对等效岩体模型进行不同组合荷载下数值模拟的研究.

（3）利用数字摄影测量系统3GSM软件获取两河口水电站左岸边坡坝址区下游节理岩体的几何特征参数.基于研究区内节理几何特征参数和三维DFN模型，建立了研究区内节理岩体三维随机DFN模型，并验证了三维DFN模型的正确性.

（4）以两河口水电站左岸边坡坝址区下游处边坡岩体为例，构建了等效岩体三维随机节理网络模型，并进行了三维节理岩体尺寸效应的研究.结果表明，研究区域内的岩体单轴抗压强度和弹性模量表现出明显的尺寸效应，并确定其REV尺寸为 8 m × 8 m × 8m，对应的等效单轴抗压强度和弹性模量分别为 25.159 MPa和19.443 GPa.