离散裂隙网络模拟在地下硐室稳定性中的应用——以甘肃北山地下实验室为例

纪景仁， 王驹， 唐振平，2， 李亚伟， 凌辉

(1.南华大学资源环境与安全工程学院， 衡阳 421001； 2. 南华大学稀有金属矿产开发与废物地质处置技术湖南省重点实验室， 衡阳 421001； 3. 核工业北京地质研究院中核高放废物地质处置评价技术重点实验室， 北京 100029； 4. 国家原子能机构高放废物地质处置创新中心， 北京 100029)

随着中国核能事业的飞速发展，高放废物的处理和处置，逐渐成为一个重大的安全和环保问题[1]。目前深部地质处置是国际公认的安全、可行的高放废物处置方式，采用的是“多重屏障系统”设计思路，把高放废物处置在数百米深的地下围岩中。一般把废物体、废物处置容器和缓冲回填材料称为“工程屏障”，把周围的地质体称为“天然屏障”[2]。中国高放废物地质处置遵循“选址-地下实验室-处置库”的三步走战略，现已进入地下实验室研发阶段，确认甘肃北山新场花岗岩场址为中国首座高放废物处置地下实验室建设场址[3]，并完成了地下实验室的初步设计。地下实验室已于2021年6月17日正式开工。

对于北山花岗岩这类坚硬的高放废物处置库围岩，岩体中的节理裂隙是主要的导水通道，核素可能作为溶质随地下水迁移，同时其几何形态和组合特征对地下硐室围岩变形和失稳起着重要的控制作用，因此研究岩体裂隙对地下硐室稳定性的影响是高放废物处置工程成功建设必不可少的一环。岩体节理裂隙作为一种地质构造运动的产物，在一定的构造应力场中生成，所以其空间分布具备一定的统计规律，研究节理裂隙几何特征和分布规律及构建与现场实际围岩具有相似统计特征的三维节理网络模型，是进行节理裂隙岩体工程力学行为计算与分析研究的基础[4]。基于这种理论思想，在20世纪80年代，Long等[5]正式提出离散裂隙网络(discrete fracture network，DFN)模型。DFN是一种依据节理裂隙的方向、位置、密度、大小统计规律，随机生成裂隙圆盘，通过展布于空间中裂隙集团来构建的错综复杂裂隙网络模型，并赋予每个裂隙圆盘相应的力学属性，从而实现各种地应力、地质力学效果的建模和拟合。近年来，随着计算机技术的快速发展， DFN模型在工程地质等领域被广泛应用。张光福等[6]利用3DEC中的DFN技术，实现了煤岩井壁的仿真模拟，并分析其稳定性，应用效果较好；罗攀等[7]使用DFN技术探究了陆相页岩裂缝的尺寸、密度和角度等储层参数对二氧化碳滤失速度的影响；Cacciari等[8]利用地面激光扫描仪扫描隧道围岩表面结构面自动生成DFN模型，并估计隧道围岩中结构面的体密度；Wang等[9]基于表征单元体(representative elementary volume，REV)和合成岩体(synthetic rock mass，SRM)建模技术提出了一种用于节理岩体开挖响应模拟的DFN-DEM多尺度建模方法，并验证了其适用性。

现以甘肃北山新场地下实验室场址地区地表、钻孔节理裂隙为研究对象，运用数理统计与空间几何手段，对节理裂隙产状、大小、密度进行描述，在此基础上采用三维离散单元软件3DEC和Monte Carlo随机分析方法，建立离散裂隙网络(DFN)，进而运用离散裂隙网络-离散元(DFN-DEM)耦合方法[10]，即用DFN模型对离散元块体进行切割，结合甘肃北山即将开展建设地下实验室工程实践，构建能反映工程岩体裂隙分布特征的等效节理裂隙岩体计算模型，并进行地下硐室的稳定性分析，这也是DFN-DEM方法首次运用到中国高放废物地质处置地下实验室硐室稳定性的研究，旨在为中国即将开展建设的高放废物地质处置地下实验室深部硐室稳定性评价提供技术支撑和方法参考。

1 模拟对象与模拟方法

根据甘肃北山地下实验室初步建设方案，地下实验室水平实验主巷道最大埋深可达地下-560 m。以地下实验室地下-560 m深度主试验区的水平硐室为研究对象，开展地下硐室稳定性研究，探索离散裂隙网络技术在北山地下实验室深部地下硐室稳定性分析中的应用效果。假设该水平巷道断面为圆形，直径7 m、长度40 m。

采用数理统计和几何分析方法，首先对甘肃北山地下实验室场址地表裂隙和钻孔裂隙数据进行规律统计，获得研究区裂隙的产状、尺寸和密度数据，建立离散裂隙网络模型；然后构建离散裂隙网络-离散元等效岩体，获得能反映北山地下实验室-560 m深度开挖硐室裂隙围岩特征的模型；最后开展等效岩体模型的数值模拟运算，得到硐室围岩的位移特征和应力状态，实现地下硐室稳定性分析与评价。

2 裂隙的几何描述及统计规律

选取甘肃北山新场地下实验室场址地区，测量3个编号为ZK1、ZK2、ZK7的工程勘察钻孔的裂隙倾角、倾向、深度数据，钻孔位置如图1所示；对新场场址勘察孔周边地表露头测量得到的4 585条裂隙产状、迹长信息进行几何特征分析。对上述裂隙统计规律进行总结，获取相关参数数据，为下一步离散裂隙网络的生成提供数据基础。

由图2可知，甘肃北山新场场址勘察孔周边地表露头裂隙可分为4个优势产状组，采用K-means聚类分析法，利用SPSS软件，分析得到各组优势产状，并划分每条裂隙的优势产状组归属，得到勘察孔周边地表露头结构面优势产状组。然后，根据得到的各组优势结构面，分别对倾向、倾角划分区间进行概率分布研究，发现各优势组倾向、倾角符合正态分布特征。接着，选用地下实验室场址地表裂隙条数大于30条，且分布较为均匀的露头，采用杨春和等[11]提出的圆形窗口法估算裂隙迹长，由频数分布发现各优势组裂隙迹长呈对数正态分布。有研究表明裂隙圆盘直径的分布形式与迹长分布形式相同[12]，因此，圆盘直径的概率分布形式也确定为对数正态分布。而后，采用伍法权法[13]估算裂隙圆盘直径分布函数的均值与标准差。

图1 钻孔位置图Fig.1 Drilling location map

图2 地表裂隙极点等密度图Fig.2 Pole and contour plots diagram of surface fractures

最后，因模拟研究对象位于地下560 m深度的空间，故利用3个编号为ZK1、ZK2、ZK7的勘察孔分别在地下500～598、500～568、500～548 m深度使用钻孔电视装置测得的共计87条孔壁裂隙数据，划分每条钻孔裂隙的优势组归属，并求得各优势组线密度P10，即单位长度内的裂隙数目。离散裂隙网络在达到钻孔定义的目标裂隙线密度P10时停止。

综上，通过对甘肃北山新场岩体裂隙特征参数分析，获得了裂隙空间各参数，如表1所示。

3 数值模型建立

3.1 DFN生成与等效岩体构建

确定建模所需裂隙参数概率模型之后，运用3DEC数值模拟软件中的离散裂隙网络(DFN)功能[14]，基于Monte Carlo随机抽样原则和线性同余法理论，生成能反映研究区节理裂隙分布特征的离散裂隙网络。根据裂隙尺寸大小，删除尺寸较小裂隙，将DFN模型简化到一个拥有许多原始模型特征的简单版本，以提高后续力学模型计算效率。本文为了模型计算需要，假设尺寸大于7 m的大尺寸裂隙对硐室围岩稳定性起控制作用。

简化后的离散裂隙网络模型如图3所示，其空间大小为40 m×40 m×40 m，以正方形形心为坐标轴原点。

图3 DFN模型Fig.3 DFN model

表1 甘肃北山新场离散裂隙网络模拟参数表Table 1 Simulation parameter table of discrete fracture network in Beishan Xinchang， Gansu Province

为了进一步模拟实现甘肃北山地下实验室地下-560 m深度，直径为7 m的圆形隧道及其所处的节理裂隙围岩条件，如图4所示，建立一个边长与DFN模型相同，为40 m×40 m×40 m，以正方形形心为坐标原点的块体，并规定其初始地应力最小水平主应力方向为X轴方向，最大水平主应力方向为Y轴方向，模型高度方向为Z轴方向。同时，在该块体中开挖一条长40 m，轴向沿Y轴，圆心为坐标轴原点，直径为7 m的水平圆形隧道，然后用离散裂隙网络(DFN)模型切割离散元(DEM)块体，构建与实际岩体有相似裂隙分布特征的等效岩体，最终实现对研究区目标岩体的数值模型构建。

图4 DFN-DEM等效岩体Fig.4 DFN-DEM equivalent rock mass

3.2 模型参数赋值

3.2.1 初始地应力赋值

根据对甘肃北山新场地下实验室预选区钻孔地应力的实测数据分析发现，研究区地应力的关系式[15]为

(1)

式(1)中：σH、σh、σv分别为最大主应力、最小主应力和竖直主应力，MPa；H为深度，m。

研究假设DFN-DEM等效岩体模型中的圆形隧道底板位于地下-560 m深度，以模拟地下实验室-560 m深度开挖硐室所处的围岩环境，该模型地应力条件按式(1)赋值。同时，模型的边界条件为：上边界为自由边界，下边界约束竖向(Z向)位移，左右边界为侧向约束。

3.2.2 岩石与裂隙物理力学参数

在DFN中的裂隙为直径大小有限的圆盘状，但是在3DEC中，块体不能被部分切割，因此在3DEC中所有的裂隙必然都比DFN模型中相应的裂隙要大。这个问题可以通过在环形裂隙的内部和外部分别赋予不同的裂隙属性来解决。如图5所示，圆盘内部为真实裂隙，赋予真实裂隙物理属性，圆盘外部为虚拟裂隙，赋予与完整岩块力学属性近似一致的虚拟裂隙物理属性[16]。

在DFN-DEM等效岩体数值模拟计算中，岩石块体采用摩尔库伦本构模型，节理裂隙采用库伦滑移模型，北山新场花岗岩体物理力学参数与裂隙物理参数取值如表2所示[17-18]。

图5 裂隙圆盘赋值属性示意图Fig.5 Schematic diagram of fissure disc assignment attribute

4 模拟计算结果

4.1 位移场分析

经过3DEC迭代运算，甘肃北山地下实验室圆形开挖巷道的整体围岩变形运动情况如图6(a)所示。硐室围岩的整体变形趋势表现为向临空面方向移动，顶拱下沉，底板回弹，且顶拱向下的变形量大于底板向上的回弹量；离硐室表面越远变形量越小，硐室围岩变形几乎关于硐口左右对称。硐室顶拱的最大下沉量为5.0 mm，底板的最大的回弹量为3.7 mm，左右边墙变形为1.0～3.5 mm。

表2 岩体及裂隙物理力学参数Table 2 Physical and mechanical parameters of rock mass and fractures

图6(b)为硐室以轴线为剖面的竖直方向(Z方向)位移云图，变形经300倍放大，可以看出，硐室在与贯穿硐室围岩的裂隙相交处发生了轻微错动，顶拱下沉量与底板回弹量较完整围岩增大，未贯穿硐室围岩的裂隙对硐室围岩的变形量影响较小，且并未扩展至贯穿硐室围岩。

图7为硐室围岩变形特征点竖向(Z方向)位移数据监测曲线，从总体趋势来看，顶拱位移量大于底板位移量，顶拱位移随计算时步逐渐变大，顶拱裂隙处围岩竖向变形可达4.3 mm；底板位移在前期逐渐增加，但是后期由于底板应力重新分布，底板位移变形表现出一定的回弹，底板裂隙处围岩竖向变形值为2.9 mm；同时可以得知，无论顶拱或底板位置，裂隙与硐室相交的薄弱处位移量大于完整围岩处位移量，这也与图6(b)位移云图表征结果相符。

图6 位移云图Fig.6 Displacement cloud map

图7 硐室围岩变形特征监测点竖向位移曲线Fig.7 Vertical displacement curve of the monitoring point for the deformation characteristics of the surrounding rock of the chamber

4.2 应力场分析

硐室开挖引起围岩应力重分布，硐壁围岩因开挖卸荷而产生位移变形，应力状态发生较大改变[19]。

图8(a)为硐室围岩的最大主应力云图，可以看出，硐室周边围岩产生较大的切向压应力集中。最大应力集中的位置出现在硐室的两侧边墙裂隙与硐壁围岩交错处，最大主应力值可达48.2 MPa，顶拱、底板为应力集中区域，但是应力集中程度与范围都小于边墙，其最大主应力应力为13.0～30.0 MPa；硐室附近应力变化显著，但距硐室5 m之外的区域，应力分布较为均匀。

图8(b)为硐室沿轴线剖面最大主应力云图，可以得知，顶拱与底板发生应力集中现象，表现为越靠近开挖面应力集中程度越大，且硐室与裂隙相交处应力发生突变，应力集中现象最为明显最大主应力量值达30.2 MPa。

5 结论

(1)通过分析统计研究区地表裂隙与深部钻孔裂隙分布规律，基于3DEC离散元仿真软件，建立离散裂隙网络，构建能反映深部地下空间岩体裂隙分布特征的等效岩体模型，进行数值模拟运算，为北山地下实验室深部地下硐室的稳定性分析提供一种简单快速的技术手段。

图8 最大主应力云图Fig.8 Maximum principal stress cloud diagram

(2)通过对甘肃北山地下实验室场址地表裂隙进行数理统计与几何特征分析，发现研究区节理裂隙可划分为4个优势组，各组裂隙倾向、倾角均符合对数正态分布，裂隙圆盘直径符合对数正态分布规律；对场址地区地下500～600 m深度3个勘察孔裂隙数据分析可知，地下实验室-560 m深度主试验区附近围岩裂隙各优势组线密度为0.222 2～0.027 7 条/m。

(3)在现有分析条件下，基于等效岩体模拟结果表明北山地下实验室地下硐室开挖后，围岩位移较小，最大变形量仅为5.0 mm，在地下实验室围岩弹性应变范围内，围岩硐室周边围岩应力集中程度不高，最大主应力为48.2 MPa，远小于北山花岗岩的单轴抗压强度141.1 MPa，说明硐室开挖后稳定性较好。