基于等效连续介质模型的单裂隙渗流数值模拟研究

高　瑜，叶　咸，夏　强



基于等效连续介质模型的单裂隙渗流数值模拟研究

高瑜1,3，叶咸2,3，夏强3

(1.云南省地质环境监测院，云南 昆明 650216；2.云南省交通规划设计研究院陆地交通气象灾害防治技术国家工程实验室，云南 昆明 650041；3.成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059)

基岩裂隙介质的空间变异性导致裂隙渗流情况十分复杂。在等效连续介质模型的基础上，建立数值分析模型研究了裂隙宽度和渗透系数对单裂隙渗流的影响。利用单裂隙立方定律对裂隙流量进行解析计算，并在单一变量研究原则下，利用MODFLOW创建了光滑平行板单裂隙介质的数值模型。据等效理论，将单裂隙宽度放大10倍的条件下，解析解计算结果比模拟结果高出6个数量级；不同裂隙宽度条件下，数值模拟结果与解析计算结果都表现出通过单个裂隙的水流量随着裂隙的宽度的增加而增加；不同渗透系数条件下，稳定流与非稳定流条件下得到的流量都随渗透系数的减小而变小。

等效连续介质；单裂隙；裂隙宽度；渗透系数；数值模拟

岩体由于应力作用、风化作用等往往存在裂隙，这使得岩体的渗流性质与一般的多孔介质有很大差别，具有如各向异性、非均匀性等复杂性与特殊性。单裂隙面是岩体复杂裂隙网络构成的基本单元，其几何特征如裂隙的方向、宽度、充填性和粗糙性等等，控制着岩体的渗透机制[1]。因此，只有首要研究单个裂隙面的渗透特征，才能为合理地预测工程岩体中复杂的渗流状态提供可靠依据[2]。

目前，单裂隙渗流问题的研究成果多基于理论解析或者现场试验，存在一定的局限性，因而将理论研究与数值模拟结果有利的结合，开展单个光滑裂隙的水流特征研究显得十分必要。而国内大部分数值模拟研究都是基于等效连续介质开展的，如隧道涌水量的计算、岩溶管道水流量[3]等。模拟结果与实际情况之间存在着怎样的误差，误差有多大，其中有哪些重要的影响因素，关于这些方面的研究则寥寥无几。本文在等效连续介质模型的基础上，利用MODFLOW进行单裂隙渗流的数值分析，对于岩体复杂裂隙网络的渗流特征研究具有显著的意义。

1　单裂隙渗流概念模型及其解析解

岩体裂隙在分布和发育上表现为不均匀性和各向异性。与之相对应的，裂隙介质的地下水渗流显然也存在非均质各向异性的特点。裂隙介质存在气态、固态、液态三相特征，气态主要是空气，并且大多存在于非饱和带中；固态称为固体骨架；液态是地下水、水与其他物质的混合物或是其他流体(如石油等)。地下水有在裂隙岩体中的存在形式较多，一般有薄膜水、吸附水、毛管力和重力水等，本文主要研究的是重力水的运动。由于岩体中的裂隙系统十分复杂，由多条平行或者交叉的单裂隙组成，而裂隙岩体的透水性比岩体基质部分的透水性大很多，故研究将岩石本身的渗透性忽略不计，可以视为流体仅在裂隙中流动[4]。

光滑平行板模型是单裂隙最基本的概念模型(如图1)，假设裂隙由两块一定间距且表面光滑的平行板构成。该概念模型将复杂的单裂隙渗流模型简单化，利用单裂隙的立方定律可以很容易的计算和预测出两平行板之间的稳定渗流量。由于光滑平行板模型构成简单，在预测裂隙中的渗流时，平行板模型和立方定律是极易让人接受的[5]。本文在光滑平行板模型的基础上展开研究。

图1　光滑平行板裂隙模型

然而，光滑平行板模型是一种较为理想的假设情况，实际上岩体裂隙并非平面或平行的空间，且岩体裂隙面的粗糙度和隙宽具有空间变异性。因此立方定律在粗糙裂隙中的渗流的预测中是不合适的，使得我们对于实际工程问题中涉及到裂隙岩体的溶质运移问题和水流运动问题难度增大。大多数情况下，要定量描述单裂隙的渗流问题需利用线性达西定律或基于达西定律的立方定律(简称LCL)。1856年，达西定律出现以后，被广泛应用在空隙介质和裂隙介质中，达西认为裂隙介质中的流体平均流速与水力梯度成线性关系，就服从达西定律。但是，许多研究成果都表明当雷诺数较大时，该线性关系就不再成立。假设裂隙内速度场是稳定、一维、一致的，线性的，构成裂隙的板面是平行且光滑的，同时是忽略裂隙宽度所造成的边界效应，则立方定律就演变成Navier-Stokes(NS)方程的简单解。基于以上假设，立方定律得出了渗流量与裂隙开度的立方成比例的结论，相关实验已经验证了该成果的可靠[6]。

Zimmerman 和Bodvarsson(1996年)[7]认为，当雷诺数小于1 时的流动，惯性力的影响将会变小，它导致渗流场非线性地依赖于水力梯度。粗糙裂隙中渗流的完全描述需要用到三维的NS方程(包含连续性方程)，NS方程仅被用来模拟具有理想几何形状的裂隙中的渗流[8]。

水在裂隙内流动又分为层流和紊流，层流的运动规律和紊流是同样的。当水流损失和与流速呈线性关系时，就称为达西定律[9]。达西定律公式：

(1)

式中：V为水流的平均速度；J为水力坡度；K为渗透系数。

对于等宽度的光滑裂隙中的线性流运动，裂隙内平均流速为：

(2)

式中：g为重力加速度；b为裂隙的宽度；vw为水的运动粘滞系数。并得裂隙的渗透公式：

(3)

则，光滑裂隙的单宽流量计算公式为：

(4)

光滑裂隙的总流量计算公式为：

Q=qL

天际线层次面积比例与评价结果间并不能视为单纯的线性相关。同时，由于天际线样本及受访样本数量的限制，当前无法确定具体的临界点比例值。但为了从统计学的角度，对天际线层次面积比例与评价结果间的相关性进行验证与说明，以下以0.442的比例值为大致临界点，将“天际线层次面积比例”评价结果散点图划分为a、b两段，分别视为两段大致呈现出线性相关的数据集，便于后文相关系数的计算。

(5)

其中，水的运动粘滞系数：

(6)

式中：u为水的动力粘滞系数，ρ为水的密度。

2　MODFLOW数值模拟试验

2.1数值模型的建立

通过查阅相关资料，比对多组数据，选取与模拟条件最接近的一组数据进行模型建立。模型为一个长20 m，宽20 m，高20 m的立方体。本次研究不考虑裂隙的倾斜角度，故将裂隙设置为水平裂隙，裂隙宽度b。基本模型中b=0.2 m。非稳定流模拟时，模拟时间总长设定为10 a，即3 650 d，模型无降雨入渗补给。根据有限差分法的基本思路，在对模型离散化处理时，每个单元体在水平方向上的长度为0.5 m，因此模型被剖分成为40列，40行的网格。垂直方向上，将模型分成7层，1、2层和6、7层的渗透系数kx，ky相同，3、4、5层的渗透系数kx，ky相同，岩体和裂隙各向异性，kz取值均设置为kx或者ky的百分之一，即垂向渗透系数比水平方向的小。其中3、4、5层即为裂隙，1、2、6、7层为岩体。在基础模型中，在稳定流的基础上，假设流经裂隙的水流不变，将左右两边均设置为定水头边界。如图2，将所有层位(第1、2、3、4、5、6、7层)的左侧水头高度设置为20 m，右侧水头高度设置为18 m。

图2　边界条件示意图

2.2基本模型数值解与解析解

在稳定流条件下运行模型，水均衡中in=37.059 m3/day，out=37.059 m3/day，平均值为37.059 m3/day，即为4.29×10-4m3/s。这个值即为通过裂隙的地下水流量。

本文根据单个裂隙的定义，建立一个用于基本研究的几何模型，模型参数如表1所示，根据公式6，算得水的运动粘滞系数为Vm=1×10-6m2/s2。通过光滑裂隙的单宽流量为653.3 m3/s。则总流量Q为20×653.3 =13 066.7 m3/s。

表1　基本模型参数表

由上可知，基本模型的公式求解结果与模拟结果相差较大，考虑到解析算法和模型中的水流为两种物理机制的水流，可以进行几何近似调整，把实际的裂隙宽度放大10倍，再放到等效连续介质模型中。所以，下文中解析算法的结果均为模型裂隙宽度的0.1倍计算所得。

2.3数值试验方案的设定

本研究主要是通过设计一个单裂隙岩体的模型，对不同裂隙宽度和不同岩体和裂隙的渗透系数比值条件下的裂隙岩体中水流情况进行模拟，并初步探讨裂隙宽度、渗透系数比值对单裂隙的渗流是否有影响，如果有影响，那么裂隙宽度和渗透系数比值分别是如何对单裂隙的渗流产生影响的。由立方定律可知，流经平行光滑单裂隙的单宽流量只受到裂隙宽度的影响，与渗透系数无关。所以此处，为了使试验具有可信度，在稳定流条件下，分别对不同裂隙宽度和不同渗透系数比值情况下做了五组数值模拟，另加五组在非稳定条件下不同渗透系数比值的数值模拟。非稳定流条件下，需要增加参数贮水系数，裂隙的贮水率值参考岩溶地区的值，而除裂隙的岩体可以考虑用砂岩的。据《地下水流数值模拟技术要求》，设定岩体的贮水系数为3e-3，裂隙的贮水系数为5e-5[10]。以下所有模拟均建立在上一节单裂隙基本模型的基础上。数值模拟试验方案如表2所示。

表2　数值试验方案表

3　数值模拟结果对比分析

3.1不同隙宽b模拟结果对比分析

通过数值模型模拟计算和解析计算，得到不同隙宽下的模拟结果和解析计算结果，见表3及图3。如表3所示。随着裂隙宽度的增加，流量也在增加。在模拟结果中，裂隙宽度与流量呈线性关系；在解析结果中，裂隙宽度与流量呈指数关系，如图3。模拟结果与解析解计算结果对比分析显示，计算结果与模拟的结果相差很大，解析解计算结果比模拟结果高出6个数量级。

表3　不同裂隙宽度流量对比表

图3　不同裂隙宽度的流量计算结果曲线

3.2不同渗透系数比值模拟结果对比分析

3.2.1稳定流模拟

从立方定律可知，流量的变化与渗透系数无关。但是从表6的模拟结果来看，随着渗透系数的减小，流量变小。表4为隙宽保持在0.2 m时，改变渗透系数，其他参数设置按照数值试验方案表(表2)模拟所得结果与解析计算结果对比表。模拟的裂隙流量随着裂隙渗透系数的减小和岩体渗透系数和裂隙渗透系数比值的增大而减小，且呈抛物线形式如图4所示，与改变隙宽的相比，流量曲线由直线变成抛物线。但是，计算的裂隙流量保持不变。

表4　不同渗透系数比值模拟流量表(稳定流)

3.2.2非稳定流

在利用Visual Modflow模拟非稳定流问题时，模拟计算可以随时停止、暂停，考虑到模型运行问题引起的误差，现增加五组不同渗透系数下的非稳定流数值模拟，模拟3 650 d的时间，将时间划分为十个步长，输出时间分别为140、309、511、754、1 046、1 396、1 816、2 319、2 924、3 650 days。

非稳定流模拟条件下，模拟结果如表6所示，随着渗透系数的减小，流经单裂隙的流量变小。从数值上看，与稳定流相比差别不大。模拟流量与渗透系数的关系如图5所示，呈线性关系。

将模型第三层的速度参数导出，求平均值，得到在改变渗透系数比值的条件下的流速如表5。模拟流速随着裂隙渗透系数的增大而增大，且呈线性关系，但是解析计算流速与渗透系数的大小无关，保持不变。

图4　不同渗透系数流量模拟结果(稳定流)

编号裂隙渗透系数kx,ky(m/s)模拟流速(m/s)计算流速(m/s)11.00E-036.67E-0421.00E-046.67E-0531.00E-056.67E-0641.00E-066.67E-0751.00E-076.67E-083.27E+01

表6　不同渗透系数比值模拟结果表(非稳定流)

图5　不同渗透系数流量模拟结果(非稳定流)

将模型第三层最后时间步长的速度参数导出，求平均值，得到在改变渗透系数的条件下的流速。解析计算的流速不随裂隙渗透系数的改变而改变，模拟流速随着裂隙渗透系数的变小而变小，且变化情况与稳定流的一致，即岩体渗透系数与裂隙渗透系数的比值越小，模拟所得裂隙中的水流速度越小。

4　结语

在单一变量研究原则的指导下，本文利用MODFLOW创建了光滑平行板单裂隙介质的数值模型，分别对宽度和渗透系数做了五组模拟，共15种数值模型。通过模拟分析了不同宽度和渗透系数比值下单裂隙水流渗透的情况，得出以下结论：

(1)分析不同裂隙宽度条件下的流量，数值模拟结果中裂隙宽度与渗流量呈一次函数关系，解析计算结果中裂隙宽度与渗流量呈指数关系，均表现出随着裂隙的宽度的增加，通过单个裂隙的渗流量增加。

(2)分析不同渗透系数比值条件下的流量，无论是稳定流还是非稳定流条件下得到的流量都随渗透系数的减小而变小。稳定流条件下，模拟结果中裂隙渗透系数与流量呈二次函数关系，非稳定流条件下，模拟结果中裂隙渗透系数与流量呈一次函数关系。稳定流和非稳定流条件下的模拟流速都与裂隙系数呈线性关系，且随裂隙渗透系数的增加而增加。

(3)以等效理论，将单裂隙宽度放大10倍的条件下，解析计算结果比模拟结果高出6个数量级。结果相差较大，说明利用等效连续介质模型的模拟方法是存在问题的，现今对于裂隙水流的研究很多模拟都是基于等效连续介质的模拟，这种方法的模拟结果值得怀疑，这将下一步研究工作的重点。

[1]艾万民.岩体卸荷应力场-渗流场耦合试验研究及数值分析[D].重庆:重庆大学.2005.

[2]宋晓晨,徐卫亚.裂隙岩体渗流概念模型研究[J].岩土力学.2004,25(2):226-232.

[3]叶咸.昆明盆地北东侧岩溶水系统探析[D].成都:成都理工大学.2014.

[4]黄勇.多尺度裂隙介质中的水流和溶质运移随机模拟研究[D].江苏:河海大学.2004.

[5]赵剑.公路隧道防排水机理及喷射防水试验装置的研究[D].重庆:重庆交通大学.2007.

[6]宋晓晨.裂隙岩体渗流非连续介质数值模型研究及工程应用[D].江苏:河海大学.2004.

[7]Zimmerman R W,Bodvarsson G S.Hydraulic conductivityof rock fractures[J].Transport in Porous Media,1996,23:1-30.

[8]Mourezenko V V,Thovert J -F,Adler P M.Permeabilityof a single fracture; validity of the Reynolds equation[J].Journal of Physics Ⅱ France,1995,(5):465-482.

[9]周志芳,王锦国.裂隙介质水动力学[M].中国水利水电出版社.2004,1:1-302.

[10]徐军祥,邢立亭.济南泉域岩溶水数值预报与供水保泉对策[J].地质调查与研究.2008,31(3):210-225.

2016-05-03

高瑜(1991-)，女，云南昆明人，助理工程师，主要从事水文地质研究。

P641.73

A

1004-1184(2016)05-0040-04