岩体裂隙网络渗流特性分析

刘 承 朱珍德 牛子豪 阙相成

(1.河海大学 岩土力学与堤坝工程教育部重点实验室,南京 210098;2.河海大学 江苏省岩土工程技术工程研究中心,南京 210098)

随着西部大开发战略的实施,水利水电工程建设在我国大规模展开.水电工程多建于水文地质条件较为复杂的环境,库区内岩石高边坡及深埋引水隧洞围岩在高地应力-渗透压力耦合作用下,其劣化机理受多类水体运移及裂隙渗流作用的影响[1].因此涉及裂隙岩体渗流相关问题的研究成为当下岩体水力学研究的热点.在自然岩体中,主要渗流通道为多条裂隙形成的裂隙网络[2-6],其结构复杂性为相关工程设计及施工带来了困难,因此准确掌握裂隙网络渗流特性至关重要.

在岩石裂隙渗流研究方面,Lomize[7]进行了平行板试验,并建立了著名的立方定律;速宝玉等[8]考虑天然岩体裂隙粗糙度的影响,结合等效水力开度对立方定律进行了修正,Wong等[9]开展了变质沉积岩应力-渗流耦合试验,得到了不同围压下单裂隙岩体渗透率的变化规律.针对单裂隙的研究,忽略了裂隙几何尺寸的突变对渗流特性的影响[10-14].而自然岩体中,裂隙交叉处水流存在明显的偏流效应[15],进而影响裂隙渗流特性.速宝玉等[16]通过交叉裂隙模型试验,论述了局部水头损失的客观存在性.朱红光[17]在局部立方定律研究中,发现在裂隙隙宽突变处存在局部水头损失,并定义了过量水压降损耗系数.

裂隙在岩体中纵横交错分布,渗流通道多存在于复杂的裂隙网络中.徐维生等[18]考虑粗糙裂隙细部特征及其对整体渗流影响,采用子单元法研究岩体粗糙裂隙网络渗流特性;Huang等[19]利用透明试样对三维透明裂缝网络进行了流体流动试验,将裂隙中流体流动行为可视化和量化,阐明了孔径非均匀性对粗糙裂隙水力特性的影响;Jiang等[20]将达西定律推广到整个裂缝网络区域,建立了针对无侧限裂隙网络渗流问题的偏微分方程,并提出了新的变分不等式公式以简化试函数的选择.以上研究从试验与理论的角度对裂隙网络渗流特性进行了描述,但均忽略了裂隙网络中各裂隙连接节点处的局部水头损失.

熊峰等[21]针对Javadi等[22]所提T 模型,考虑黏性压降和局部压降,提出了新低速下粗糙裂隙非线性渗流模型,并通过饱和渗流试验验证了模型的正确性.姜清辉等[23]利用伽辽金法求解雷诺方程和Forchheimer方程,建立了三维离散裂隙网络中非线性流动的数值计算方法,比较了3种交叉模型的数值模拟结果与实验数据,两者吻合较好.宋羿[24]构建了分叉-交叉裂隙试验室物理模型,定量分析了不同进出口方式的交叉裂隙渗流特性.熊峰等[25]通过对4种裂隙交叉点的流动试验,研究了裂隙交叉点对流体非线性行为的影响,进行了相应的数值模拟.采用CT扫描方法,对含不同类型交叉裂隙的岩石试样进行断层扫描,建立了裂隙网络模型,进行了试验和数值模拟,研究了网络尺度下的非线性流动行为.通过以上研究,裂隙网络交叉节点渗流特性的研究取得一定成果,但仍不够系统,有待进一步研究.

为了进一步研究局部水头损失对裂隙网络渗流特性的影响,本文将裂隙网络电路化,引入渗流阻模型[26],针对裂隙网络中裂隙的基本连接形式,建立等效过量压降损耗阻模型;通过等效过量压降损耗阻描述串并联连接裂隙的过量水压降损耗,并结合多裂隙复合连接物理模型试验进行验证.

1 等效渗流阻模型介绍

1.1 岩石裂隙立方定律

隙宽为e的单位厚度平行板裂隙达西流,其裂隙流量满足立方定律:

式中:q为裂隙流量;μ为液体动力黏滞系数;ΔP为裂隙两端压力差值;ΔL为裂隙长度.

天然岩石裂隙中,裂隙壁面粗糙不平,难以满足平行板模型的条件.因此,需对立方定律进行修正,才能应用于实际裂隙情况.Lomize[7]、Louis等[27]相继对粗糙裂隙的水流特性进行了研究,通过增加粗糙度修正系数C对立方定律进行了修正:

1.2 等效渗流阻模型

由于大多数完整岩块的渗透性极为微弱,可认为流体仅在岩体的裂隙网络中流动.裂隙网络是由有限个单裂隙连接而成,最基本的连接方式包括串联和并联两种,如图1所示.

图1 裂隙的两种基本连接方式

1)串联关系

平行板裂隙串联情况如图2所示,其上下两端裂隙中流量相等,结合立方定律,可得

式中:q为裂隙总流量;q1、q2为上下两段裂隙的流量;ΔP1、ΔP2为上下两段裂隙的压降;C1、C2为上下两端裂隙的粗糙度修正系数;ΔL1、ΔL2为上下两段裂隙的长度;e1、e2为上下两段裂隙的隙宽.

图2 串联平行板裂隙

结合上下两段裂隙的总压降ΔP=ΔP1+ΔP2,得

2)并联关系

如图3所示的平行板并联情况,并联裂隙的总流量等于其内裂隙的流量和,即

图3 并联平行板裂隙

在并联裂隙中,其总压降ΔP=ΔP1=ΔP2,则式(5)可以写成:

将式(4)和式(6)与模拟电路串并联情况下总电压与总电流关系式相类比,可以发现12μCΔL/e3这一部分在渗流问题中与电阻在电路问题中所起到的作用相似,类比结果见表1.因此,引入参数等效渗流阻[26]Rf=12μCΔL/e3表征裂隙网络的阻尼作用.从式中可以看出,渗流阻与裂隙长度、粗糙度呈正相关性,与裂隙宽度呈负相关性.

表1 电路问题与渗流问题相似物理量表

2 等效过量压降损耗阻模型

实际情况下,裂隙流量除受裂隙渗流阻影响,还与相连裂隙宽度比有关.朱红光等[17]在应用局部立方定律计算粗糙裂隙流量时,发现裂隙隙宽突变处存在过量水压降损耗,即裂隙宽度突变处的实际压降值大于计算值.并定义了过量水压降损耗系数ξ:

式中:ΔPc表示通过局部立方定律计算得到的水压降;ΔP表示实测水压降.

结合过量水压降损耗系数,对串联模型进行修正.将式(4)代入式(7),得

显然,结合过量水压降损耗系数修正时,需明确节点两端裂隙的实际水压降.对于单个节点的串联裂隙较容易实现,但对于包含多节点连接的复杂裂隙网络,计算网络内各个节点两端的水压降过程繁琐.因此,本文结合等效渗流阻模型,定义了等效过量压降损耗阻Rx:

式中:Rf1和Rf2为串联裂隙两段裂隙的渗流阻;R总为串联裂隙总渗流阻力.定义过量压降阻率η:

通过等效过量压降损耗阻来修正等效过量水压降的影响,简化了多节点连接裂隙网络渗流特性的分析.下面针对几种常见的裂隙连接模型分析其过量压降阻率.

1)不同隙宽裂隙串联连接

不同隙宽裂隙串联包括扩张型和收缩型两种方式.如图4所示,渗流方向由A 端到B 端为扩张型,由B端到A 端则为收缩型.

图4 串联裂隙

计算所得压降差为:

过量水压降损耗系数为:

将水压降无量纲化:

结合式(10)、(12)～(14)得:

结合试验及模拟分析,发现收缩型裂隙的压降损耗远大于扩张型裂隙的压降损耗,可忽略扩张型裂隙的过量水压降损耗,并得到收缩型裂隙满足下式:

将式(17)代入式(16),可得收缩型裂隙的过量压降阻率η:

根据式(18)不难看出,两个不同隙宽裂隙串联的过量水压降损失与雷诺数Re,窄裂隙隙宽和裂隙长度比值e1/L以及隙宽比e1/e2有关.雷诺数越大,即流速越大时,其过量水压降损耗越大;隙宽比值越小时,其过量水压降损耗越大.

2)串并联模型

裂隙的串并联连接组合方式繁多,本文选取一组基础的连接模型进行研究.如图5所示,隙宽为e1的裂隙分别与隙宽e2和e3的裂隙串联,隙宽e2和e3的裂隙并联.将隙宽为e1、e2和e3的裂隙分别编号为裂隙①、裂隙②和裂隙③.

图5 串并联裂隙

基于立方定律可得,裂隙总流量式为:

串并联裂隙由于存在两个并联的隙宽突变节点,且每个节点处的过量水压降损耗系数不同,较难得到准确的过量水压降损耗.传统方法只考虑主要流动通道的过量水压降损耗[28],即忽略隙宽较小的裂隙②的影响,仅分析裂隙①和裂隙③的串联裂隙.对于e3≫e2的情况下,忽略e2的影响,式(19)即简化为式(11),因此裂隙②可以忽略;但当e2和e3相近时,误差就会较大.

为准确分析e2和e3相近情况下的裂隙渗流特性,结合等效过量压降损耗阻模型进行研究.根据式(18),裂隙②通过影响裂隙③的雷诺数来影响裂隙的过量压降损耗阻,进而影响裂隙渗流特性.

本文采用如下分析方法:分别计算两条支路裂隙雷诺数Rf2、Rf3,根据裂隙支路雷诺数计算两条支路的过量压降损耗阻Rx2、Rx3,结合两部分等效压降损耗阻,综合考虑裂隙渗流特性.

以裂隙①与裂隙②串联为例,裂隙雷诺数为:

式中:Re为裂隙①的雷诺数,Rf2、Rf3分别为裂隙②、裂隙③的等效渗流阻.

结合式(10)、(18)和(21),可得裂隙②支路的等效压降损耗阻:

同理可得,裂隙③支路的等效压降损耗阻如式(23)所示:

结合所得各支路等效压降损耗阻,将串并联裂隙整体电路化,绘制等效电路图,如图6所示.

图6 等效电路图

根据串并联电路特性,总的等效压降损耗阻为:

代入式(22)、(23),且当e2、e3≪e1,e2、e3接近时,式(24)可化简为:

由式(25)可知,串并联电路的总等效损耗阻与支路等效渗流阻和支路隙宽有关,而与主路径(裂隙①)无关.这是因为e2≪e1且e3≪e1,而渗流阻的数量级是e-3,所以可以忽略主路径的影响.

3 等效过量压降损耗阻模型的验证

朱红光等[28]开展了一系列基本形式的多裂隙连接渗流特性试验研究,取得了良好的实验数据.基于上述试验数据,对本文提出的理论公式进行可行性验证.其试验模型概念如图5所示,其中,e1,e2,e3分别为7mm、1mm、2mm;L为40mm;入口为左端,水压为ΔP;出口为右端,水压为零;总流量为q.

采用传统方法,即忽略裂隙②,仅考虑裂隙①与③串联连接,计算结果见表2.

表2 忽略裂隙②计算结果表

基于等效过量压降损耗阻模型,综合考虑裂隙②的影响,计算结果见表3.

表3 基于等效过量压降损耗阻模型计算结果表

将两种方法的流量计算结果汇总并与试验数据进行比较,见表4.

表4 计算结果与试验数据统计表

分析表4可以看出,相较于传统忽略次要裂隙方法的计算结果,基于等效压降损耗阻模型提出的计算方法所得结果与试验数据拟合效果更佳.当雷诺数较小时,次要裂隙水流量极小,对裂隙过量压降损耗的影响可忽略,因此两种方法计算结果相近.随着雷诺数的增大,基于渗流阻模型的计算方法显然与试验数据更加吻合.分析计算所得流量小于试验值的原因,改变式(24)形式,得支路裂隙总渗流阻RZ:

为了便于分析,定义裂隙②与裂隙①雷诺数之比为雷诺数分配系数α,即,α=Re2/Re.结合式(22)、(23)和(26),分别代入Re=200,500,800,1000,绘制在不同雷诺数情况下RZ与雷诺数分配系数α的关系曲线,如图7所示.

图7 RZ 与雷诺数分配系数α 关系曲线

由图7可以看出,不同雷诺数情况下,RZ都随着雷诺数分配系数α的增大而减小.在进行雷诺数分配时,仅按照渗流阻比例进行分配,忽略了裂隙过量压降损耗阻等的影响,导致计算所得的雷诺数分配系数偏小,从而使裂隙总的渗流阻力偏大,计算所得流量偏小.

4 算例

结合上文建立的等效过量压降损耗组模型,对单位厚度六棱柱柱状节理试样进行渗流分析.分析对象如图8所示,由裂隙1～12及节点A～H 构成的联通裂隙网络,其中所有裂隙长度为80mm,裂隙1、2、6、7、11、12的隙宽为6mm,裂隙3、8的隙宽为2mm,裂隙4、9的隙宽为3mm,裂隙5、10的隙宽为4mm.设总压降ΔP=40kPa,水温为20.5℃,则μ=0.998×10-6kPa·s.

图8 单位厚度六棱柱柱状节理试样

分析该裂隙网络可得,仅在节点A,C,F,H 处存在隙宽减小的情况,因此只考虑该四处压降损耗即可.将裂隙网络电路化,得到等效电路图,如图9所示.

图9 裂隙网络等效电路图

忽略等效过量压降损耗阻,结合电路学知识及电路对称性,计算各裂隙流量及雷诺数,计算结果见表5.

表5 忽略局部压降计算结果表

将裂隙雷诺数代入式(18),求得裂隙过量压降损耗阻,见表6.结合所求过量压降损耗阻,计算各裂隙流量,计算结果见表7.

表6 过量压降损耗阻表(单位:10-6s·kPa·mm-3)

表7 考虑局部压降计算结果表

分析计算结果可知,过量压降损耗阻与裂隙突变程度有关,相连裂隙隙宽相差越大,过量压降损耗阻越大.忽略局部压降,会使裂隙网络整体流量偏大,但也有部分裂隙流量偏小,如裂隙5、7、10,这是由于裂隙2和5、裂隙7和10的隙宽相差最小,产生的过量压降损耗阻最小,在考虑局部压降时,其他裂隙增加的阻力更大,导致裂隙5、7、10的流量增大.

5 结论

本文基于渗流阻模型,考虑了不同连结情况下裂隙节点处多余水压降的损失,建立等效过量压降损耗阻模型,并与试验结果进行对比验证,得出以下结论:

1)在一个区域的裂隙网络中,将裂隙网络电路化,考虑等效过量压降损耗阻,对整体裂隙网络进行分析,简化了裂隙网络渗流的计算过程.

2)本文提出的等效过量压降损耗阻模型,分析了串并联裂隙中次要裂隙对整个串并联裂隙渗流特性的影响,通过与试验数据对比,相较于忽略次要裂隙的方法,本模型的计算结果更为准确.

3)定义了雷诺数分配系数,分析了不同雷诺数下,支路总渗流阻RZ雷诺数分配系数α的关系曲线,阐明了理论模型误差存在的原因.

本文提出雷诺数分配系数α分析了理论模型误差存在的原因,因此,在后续研究中,可以针对雷诺数分配系数α进行研究,探究如何对α进行修正,建立更为准确的理论模型.