摘 要：【目的】岩体的裂隙特征直接影响其渗流特性，随着隧道工程中与岩体裂隙渗流有关的突涌水问题日益突显，为提高隧道建设与维护的安全性及稳定性，建立一套隧道掌子面岩体裂隙渗流量的预测计算程序。
【方法】基于MATLAB平台，利用Monte-Carlo随机模拟方法编制岩体裂隙二维网络模型生成程序及稳定渗流数值计算程序，研究了裂隙岩体表征单元体积(REV)的尺寸效应及渗透张量的确定方法。
将该程序应用于某地铁隧道，在现场调查及裂隙参数统计的基础上，对隧道掌子面裂隙岩体稳定渗流量和渗透系数进行预测。
【结果】结果显示：(1)掌子面岩体裂隙方向角服从正态分布，迹长及间距服从负指数分布；(2)裂隙岩体表征单元体(REV)的最佳尺寸为裂隙迹长均值的14倍；(3)计算了裂隙岩体的渗透系数椭圆、渗透张量及渗透主轴，分析结果表明该隧道掌子面裂隙岩体具有强烈的渗透各向异性特征。
【结论】结果表明：岩体裂隙的网络模型对于隧道工程岩体的水力学问题研究有不可替代的作用。
所编程序可以预测隧道开挖过程中掌子面的渗流量和渗透主方向，为采取经济、合理的地下水控制措施提供依据。
关键词：裂隙岩体;Monte-Carlo方法;MATLAB;裂隙网络模型;岩体渗透张量;作者简介：雷林(1999—),男，硕士研究生，主要从事岩石力学方向的研究。
左双英(1977—),女，教授，博士研究生导师，博士，主要从事地质工程方面的教学和科研工作。
基金：国家自然科学基金(42167025);贵州省科学技术基金重点项目([2020]1Z052);贵州省科技重大专项计划(黔科合重大专项字[2018]3011);引用：雷林, 左双英, 田娇, 等. 岩体裂隙网络渗流模型及隧道掌子面渗流预测研究[J]. 水利水电技术(中英文), 2023, 54(8): 104- 114. LEI Lin, ZUO Shuangying, TIAN Jiao, et al. Study on fracture network model and permeability characteristics of rock mass in tunnel face [J]. Water Resources and Hydropower Engineering, 2023, 54(8): 104- 114.0 引 言近三十年来，地下空间的开发利用在快速交通、节水节能以及城市地下综合体等方面的应用取得了巨大成功,但在建设与运营过程中仍面临着一些技术问题，如涌水突泥、地面沉降、管道裂缝等。
岩溶区碳酸盐岩层理及节理裂隙发育，使得岩体内部存在着显著的物理力学不均衡性和各向异性，对隧道围岩地下水渗透方向、水力效应及失稳方式起到了控制作用，隧道施工后的围岩局部大变形、渗漏水等不稳定状况也常出现于这些薄弱部位，因此在进行隧道掌子面及围岩稳定分析时考虑裂隙引起的渗透各向异性已经在学术界和工程界达成了共识。
近年来，随着裂隙性油气田的研究以及高放射核废料深埋处置工作的开展，对大型工程中的岩体裂隙渗流理论和模型的探讨受到广泛关注，并形成了热点问题。
目前，基于达西定律的孔隙均质岩体渗流、基于立方定律的单裂隙岩体渗流研究相对比较成熟,而对于岩体裂隙网络建模、连通性判别及非线性渗流计算等方面的研究仍处于探索阶段。
在裂隙模型生成方面，李葳等基于UDEC软件模拟生成岩体裂隙二维模型，进行了渗透性REV值及等效渗透系数的研究；叶懿尉使用3DEC软件简化了节理的建模并探究了节理岩体的力学性质随着节理简化的变化规律；ZHAO等基于幂律模型对裂隙网络进行建模，分析了隧道开挖和排水引起的扰动机制；吴顺川等、武娜等基于Monte-Carlo方法建立了岩体裂隙网络模型，但缺乏对模型水力学特性的研究。
在裂隙网络连通性判别及处理方面，王恩志等以“去除法”与“分离法”相结合处理了非连通裂隙网络；张有天根据逾渗理论将裂隙网络中的非连通裂隙及网络簇删除；张俊提出了“地质熵”的概念来表征二维裂隙网络的连通性。
针对岩体裂隙网络的渗流算法研究，王者超等提出了平行流与辐射流两种流动形态的模型粗糙度计算方法，并借助COMSOL软件研究了两者的渗透特性；薛振晓基于面流模型，建立三维裂隙网络，将裂隙面划分成三角形单元进行渗流计算；申林方等采用格子Boltzmann方法研究了粗糙岩体裂隙几何特征及对渗流特性的影响，模拟了不同形式的裂隙渗流。
以往的研究大多根据岩体裂隙几何参数进行平面展布后形成裂隙网络模型，没有考虑裂隙方位角的变化，计算结果一般等效为均质渗透系数，忽视了裂隙方向性导致的岩体渗透各向异性。
在隧道涌水量预测方面，一般采用水均衡法、水文地质比拟法、地下径流模数法、水动力法等,这些方法预测结果的精度也取决于水文地质单元边界的准确划分及水力参数的精准取值。
随着地下工程在特殊领域应用需求的增加及安全级别的提升，对岩体裂隙网络模拟的精细化要求越来越高，本研究对施工期地铁隧道掌子面进行岩体裂隙的现场采集及参数统计分析，基于MATLAB计算平台，利用Monte-Carlo随机理论编写岩体裂隙二维网络模型及稳定渗流模拟程序，研究裂隙连通性判别方法、表征单元体(REV)尺寸效应及反映岩体渗透各向异性的渗透张量确定方法，预测掌子面渗流量，为隧道建设采取经济、合理的地下水控制措施提供依据。
1 裂隙网络模型的生成1.1 基于Monte-Carlo法生成随机数根据统计得到的确定性样本规律，用Monte-Carlo方法随机生成符合特定规律的模型。
标准均匀分布的随机数是生成符合其他概率分布函数随机变量的基础，本研究以混合线性同余法生成均匀分布的随机数式中，a为乘子；xi(i=0)为初值；modM表示对模数M取余；c为常数增量，c与M互为质数，并且M和a-1能同时被同一个质数或4整除；ui为区间[0,1]均匀分布的随机数。
1.2 裂隙二维网络模型的生成以Monte-Carlo方法，基于MATLAB平台编写代码生成隧道掌子面岩体裂隙二维网络模型，相关步骤如下：(1)统计建模所需参数，即裂隙的密度、倾向、迹长、隙宽及其概率分布函数；(2)确定二维模型的生成域和分析域；(3)在MATLAB平台输入每组裂隙的相关参数，并以矩阵形式储存；(4)裂隙中心点坐标(xi,yi)按均匀分布、方向角及迹长按测量数据统计出来的分布形式随机生成；(5)在规定的区域内模拟生成所有裂隙，即可构成生成域的二维裂隙模型。
1.3 裂隙网络连通性优化假定岩块不透水，渗流只发生在连通的裂隙网络中，裂隙网络由服从一定函数分布的随机裂隙组成，因此计算渗透系数时需要对生成的二维网络模型进行连通性判别及优化处理。
优化该模型的关键在于剔除没有形成水流连续通道的裂隙或裂隙网络簇。
引入复杂裂隙网络的衔接矩阵A={aij}M×N、回路矩阵C={ckj}L×M,将裂隙图形用数学矩阵的方法表示出来裂隙节点与线单元之间的衔接关系由矩阵A定义，矩阵C描述了裂隙网络系统中线元之间的组合关系，基于式(3)和式(4)算法，以MATLAB软件编程实现对裂隙网络中线单元以及节点的自动编号并赋值于矩阵中，结合COLIN等提出的评估二维随机裂隙网络连通性的参数η如下式中，εmean和ρseg分别为裂隙网络中裂隙平均密度与段密度；A为整个区域面积。
η≥1时，表明裂隙路径连通性较好，且η越大，表明连通性越好，裂隙的渗透系数越大。
由式(5)可得到裂隙优化模型如图1所示。
图1 裂隙模型及其优化示意(左右为水头边界)2 二维裂隙网络渗流基本方程在复杂的裂隙网络中，采用有限元方法分析网络水力学稳定流效果较好，二维的裂隙网络由一维单元组成，且其可以分解为多个2结点单元，在模型分析域中进行渗流分析时，只需考虑有水流通过的单元。
如图2所示，虚线圆为以i节点形成的表征单元域，ckj为某回路，由水均衡原理，进、出表征单元域的水流量相等，根据域内地下水均衡方程式，可得第一类边界上的渗流量为式中，Q3为第一类边界上的渗流量；A1、A2、A3分别为衔接矩阵，描述裂隙网络系统中线元和节点的衔接关系；H1、H2、H3为节点水头矢量与边界水头矢量；T为对角矩阵，T=diag(T1,T2,...,TM),其中，Tj=ρgbλb3jhj/μLj。
图2 裂隙单元均衡域由此得到模型的总渗透量Q3,将模型等效为一个平行板模型，由立方定律结合达西定律，可推导出模型的等效渗透系数K式中，a为等效水力隙宽；Q为通过光滑平行板的流量，其值等效为Q3;ν为水的运动黏滞系数，本文取1.003×10-6m2/s;g为重力加速度，取9.8 m/s2;J为沿平行板方向的水力坡降，可根据模型尺寸和水头条件求得；u¯为裂隙水流等效平均流速。
3 基于二维裂隙网络模型的渗流预测3.1 二维裂隙网络渗流计算3.1.1 边界条件设置确定裂隙网络模型的边界条件，即计算模型(见图3),假定正方形裂隙二维模型(分析域)的左边界为高水头边界，右边界为低水头边界，上、下边界为同等变化的梯度水头边界。
图3 二维裂隙网络计算模型3.1.2 生成域和分析域以裂隙网络系统二维稳定渗流数值解计算矩阵为基础，基于MATLAB平台编写相关程序进行隧道掌子面岩体裂隙二维模型(分析域)的稳定流渗透系数计算。
计算之前，首先需要确定裂隙网络模型的分析域尺寸，由于二维裂隙网络模型的尺寸效应,为了使研究结果具有意义，应设置分析域的尺寸大于REV值。
基于某地实测的岩体裂隙数据(见表1)模拟生成二维模型及其渗透系数计算，分析迹长均值与分析域的关系，最终确定裂隙二维模型的尺寸。
如图3所示，假定计算模型的边界条件为：高水头边界H0=10 m, 低水头边界Ht=0 m, 梯度水头边界Hs从左到右由10 m梯度变化到0 m, 张有天提出若二维裂隙模型为矩形域，实际裂隙网络生成域尺寸应为分析域的2倍，如图4所示，虚线框内为分析域。
基于此，采用控制变量法，改变生成域和分析域的尺寸模拟生成岩体裂隙模型并计算渗透流量和渗透系数。
如图5所示，可明显看出，进行岩体水力学分析时，正方形分析域的尺寸会影响岩体裂隙二维模型的渗透流量和渗透系数，假定边界条件不变的情况下，随着分析域尺寸逐渐变大，流量和渗透系数逐渐趋于稳定。
当正方形分析域尺寸为7 m×7 m(生成域14 m×14 m)时，即正方形分析域的边长约为迹长均值的14倍时(最大迹长均值为0.49 m),岩体裂隙二维模型的渗透流量和渗透系数达到稳定状态。
图4 分析域与生成域的尺寸关系图5 裂隙模型分析域尺寸与流量、渗透系数的关系4 工程应用4.1 研究区概况研究区属高原溶盆区，地势空间大致格局是北、西、南三面高，东面低，北、西、南三面略向东倾斜呈簸箕状。
黄望隧道区间属河流槽谷及谷坡山脊地貌，槽谷宽缓低平，地面高程1 070～1 090 m, 山脊坡陡，主要构造为贵阳向斜。
右线YDK39+148隧道段采用矿山法开挖，揭露的掌子面主要为破碎状中风化泥质砂岩夹灰岩，现场勘察发现，隧道掌子面存在裂隙渗水的迹象(见图6中反光处)。
图6 掌子面渗水情况4.2 岩体裂隙现场调查参考国际上通用的岩体裂隙测量方法，采用5个指标(产状、间距、迹长、隙宽、充填物)对岩体裂隙进行测量和描述，掌子面岩体裂隙十分发育，采用人工测量的方法精度较高，但测区范围有限且耗时耗力可能会影响施工进度，再考虑到人工测量危险系数高，因此在满足测量精度及特殊地形的勘测要求下，基于岩体裂隙的传统测量方法，制定了两种裂隙采集互补方案。
(1)测网法：人工测量出隧道掌子面岩体的4条典型裂隙(见表2);地质罗盘测量裂隙产状(倾向、倾角);皮尺测量裂隙的迹长；钢尺测量裂隙视间距；隙宽是裂隙网络水力学分析最主要的参数，采用高精度塞尺(0.02～1.0 mm)多次测量某裂隙不同位置处的隙宽(见图7),取平均值作为该裂隙的隙宽值。
图7 塞尺测量隙宽(2)高清照片采集：布置1.5 m高的三脚架、隧道灯位于隧道掌子面正面6 m处，采用佳能80D单反相机拍摄掌子面正面图片，该照片将用于导入CAD描图。
从不同角度、距离拍摄掌子面图片约50张，辅助CAD描图。
4.3 岩体裂隙数据统计分析为获取该隧道掌子面岩体裂隙的分布及数据特征，选取拍摄质量较好的掌子面正面图片导入CAD软件，对裂隙按1∶1进行矢量化。
(1)选取图8中320 mm×180 mm绿色方框区域为ROI(感兴趣区域),以MATLAB统计分析ROI中的裂隙参数如表3所列。
图8 掌子面正面图CAD处理(单位：mm)(2)对表3数据进行拟合分析，分析过程中剔除个别误差较大的数据，结果如图9所示。
拟合分析结果显示，该隧道掌子面岩体裂隙的倾向(方向角)服从正态分布，迹长及间距均服从负指数分布。
图9 裂隙几何参数拟合分析(3)为更明显地反映研究区域各组裂隙的发育程度、优势方位，拟采用表3数据按空间方位间隔10°分组，绘制裂隙等密度图[见图10(a)]和裂隙倾向玫瑰花图[见图10(b)]。
可见，该掌子面岩体发育4组优势裂隙，与人工现场测量的典型裂隙误差较小，可作为建模及计算的基础数据。
图10 裂隙发育等密度图及玫瑰花图4.4 研究区裂隙模型生成裂隙模型生成：(1)统计建模所需裂隙参数如表3所列；(2)确定裂隙模型的生成域及分析域，裂隙最大迹长均值为1.86 m, 根据14倍的REV尺寸效应要求，确定裂隙模型的正方形分析域边长为26 m, 生成域边长为52 m; (3)将表3数据输入基于MATLAB平台编写的程序，自动生成掌子面ROI裂隙二维模型部分区域图，截取尺寸为8 m×8 m如图11所示。
图11 岩体裂隙8 m×8 m生成域网络模型对比CAD图像处理图可知生成的裂隙二维模型与实际ROI范围内裂隙的分布情况相似度较高。
4.5 稳定渗流计算考虑到开挖过程中地下水位的变化，进行隧道掌子面岩体裂隙二维模型(分析域)的渗透流量及渗透系数计算时，设计高水头边界H0=10 m、20 m、30 m三组水头值、低水头边界Ht=0 m, 以及梯度水头边界：从左到右由H0梯度变化到0 m。
基于Monte-Carlo方法随机生成的裂隙二维模型具有细微的差别，为提高计算结果的精确性，对每组水头边界条件进行三次重复计算。
表4结果显示，随着计算边界条件的改变，二维裂隙网络的渗透流量及渗透系数会发生相应的变化，随着水头的增加，水力梯度增大，分别为0.38、0.77、1.15,模型渗透流量Q逐渐增加，由式(8)计算得到等效渗透系数K逐渐减小，分别为2.97×10-4 m/s、2.59×10-4 m/s、1.81×10-4 m/s。
模型每天的渗透总流量Qh分别为188.71 mL、221.72 mL及268.18 mL,对比图6掌子面渗水情况可知，该计算结果与实际情况相符。
4.6 裂隙岩体渗透系数各向异性特征自然条件下，岩体裂隙的方向性会导致同一点在不同方向的岩体渗透性存在差异，称为渗透的各向异性，通常用渗透张量表示，渗透张量椭圆可简单直接地表现出岩体的渗透特征。
为探究研究区岩体是否存在渗透张量及其裂隙网络是否可等效为连续介质模型，拟设计以下试验来拟合渗透椭圆。
保持其他裂隙几何参数不变，固定水头边界为：高水头边界H0=10 m, 低水头边界Ht=0 m, 梯度水头边界Hs从左到右由10 m梯度变化到0 m, 通过MATLAB软件编程实现将所建模型从原始坐标绕x轴逆时针旋转，使得裂隙倾向(方向角)改变，设置每次转动30°,共旋转12次，计算结果如表5所列。
将渗透系数的平方根倒数和角度分别作为极坐标系中的点，拟合成一个椭圆(见图12),渗透椭圆方程如下图12 不同倾向裂隙岩体渗透张量变化规律长轴a=1/√K2,短轴b=1/√K1,由渗透椭圆的基本性质可确定渗透主值K1、K2及其渗透方向，得到渗透张量表达式由图12,主渗透系数K1与x轴的夹角约为15°,综合图8与图11来看，符合实际情况，K1=3.024 6×10-4(m/s),K2=1.384 1×10-4(m/s),则由式(10)、式(11)计算得到渗透张量为验证该渗透张量是否能够描述生成裂隙岩体的渗透性，需要通过均方误差来判断拟合结果是否符合要求，用下列公式计算式中，Ktensor为θ方向上模拟得到的等效渗透系数；N为旋转方向个数，本文取12;Ksim是通过椭圆拟合得到θ方向上的等效渗透系数。
由式(12)计算出RMSnorm=0.070 82;当RMSnorm≤0.2时，拟合效果较好,说明研究区裂隙岩体存在渗透张量。
图12描述了隧道掌子面不同倾向(方向角)裂隙岩体二维模型渗透张量的变化规律，从图中可以看出，K1与K2的差值会随着裂隙倾向绕“x轴”逆时针旋转的角度不断变化，说明裂隙岩体渗透具有各向异性特征，旋转约15°时，K1与K2差值最大，各向异性特征最明显。
通过裂隙网络建模及裂隙岩体等效连续介质渗流的计算，可以预测隧道开挖过程中掌子面的渗流量和渗透主方向，在施工过程中可以采取“防、排、截、堵”等相关措施处理地下水带来的危害，为安全施工提供科学依据。
5 结 论基于某地铁隧道掌子面的岩体裂隙实测数据，采用Monte-Carlo方法，以MATLAB软件编写代码生成了岩体裂隙的二维网络模型，并依据二维裂隙网络渗流的基本方程，编写相关程序实现了该模型渗透流量和渗透系数的自动计算。
得出以下结论：(1)通过改变岩体裂隙网络模型的尺寸，计算渗透系数，当渗透系数达到稳定状态时的尺寸为岩体表征单元体(REV)尺寸，其大小为岩体裂隙迹长均值的14倍。
(2)采用裂隙的两种互补测量方案对工程区隧道掌子面岩体裂隙统计得到4组优势裂隙，其倾向(方向角)服从正态分布，迹长及间距均服从负指数分布。
通过MATLAB编程自动生成的隧道掌子面岩体裂隙二维模型中，裂隙的分布情况与CAD描图及拍摄的掌子面正面图(ROI)符合度较高，总体模拟效果较好。
(3)对裂隙网络模型进行了连通性判别及优化处理，编程计算出了该隧道掌子面岩体裂隙的渗透流量及渗透系数。
(4)为研究裂隙渗透的各向异性，拟合出了隧道掌子面裂隙岩体在不同方向的渗透张量椭圆，由渗透张量椭圆的性质计算出了渗透主轴及渗透张量，证明了裂隙岩体渗透具有各向异性特征。
水利水电技术（中英文）水利部《水利水电技术（中英文）》杂志是中国水利水电行业的综合性技术期刊（月刊），为全国中文核心期刊，面向国内外公开发行。
本刊以介绍我国水资源的开发、利用、治理、配置、节约和保护，以及水利水电工程的勘测、设计、施工、运行管理和科学研究等方面的技术经验为主，同时也报道国外的先进技术。
期刊主要栏目有:水文水资源、水工建筑、工程施工、工程基础、水力学、机电技术、泥沙研究、水环境与水生态、运行管理、试验研究、工程地质、金属结构、水利经济、水利规划、防汛抗旱、建设管理、新能源、城市水利、农村水利、水土保持、水库移民、水利现代化、国际水利等。