粗糙裂隙各向异性对非达西渗流特性的影响

2023-02-12马亚楠，杨志兵，熊小锋，陈益峰

人民长江 2023年1期

马亚楠，杨志兵，熊小锋，陈益峰

(1.武汉大学水资源与水电工程科学国家重点实验室，湖北武汉 430072； 2.武汉大学水工岩石力学教育部重点实验室，湖北武汉 430072)

0 引言

在高坝枢纽、高水头抽水蓄能电站、煤矿开采和深埋隧洞等工程活动中，受渗透水压和裂隙形貌的影响，岩体渗流具有强烈的非线性，发生渗漏、突透水和渗透破坏的风险急剧增大[1-4]，因此非线性渗流对于工程建设和运行稳定有着重要影响。岩体渗流主要是通过裂隙结构面进行的。受地质沉积环境、温度变化和地质构造作用等因素的影响，天然裂隙面存在各向异性特征[5]，使得裂隙渗流模拟及预测存在不确定性。研究各向异性粗糙裂隙中的非线性渗流规律，对于预测宏观尺度裂隙岩体渗流，预防隧洞涌水、内水外渗、煤矿突水和瓦斯突出等事故具有重要理论意义和参考价值。

近年来已有许多学者开展了关于裂隙岩体渗流的研究[6-7]。基于平行板假设的立方定律被广泛用于描述裂隙介质中的渗流[8]。以往的研究针对天然裂隙粗糙特性的影响，提出了一系列修正的立方定律。偏离线性规律的非达西渗流现象在数值模拟、室内外试验和工程应用中广泛存在[9-10]。研究发现裂隙渗流的非线性受到法向应力[11]、剪切错动[12]、接触面积[13]等因素的影响，其本质上是裂隙几何形状和开度的变化对裂隙渗流的控制作用。非达西渗流细观机制和宏观表征的研究得到了普遍关注。在细观机制方面，理论分析结果显示当惯性效应占主导作用时渗流出现非达西特征[14]；此外，数值模拟结果表面涡旋区的产生和发展是非达西渗流发生的主要原因[15]。在宏观表征方面，Izbash公式和Forchheimer公式常用于描述非达西渗流的压力梯度与流量关系。其中Forchheimer公式物理意义明确，参数取值稳定性强，被广泛用于裂隙介质非线性流的表征中。目前主要通过建立裂隙形貌与非达西渗流公式参数的关系来对非线性系数进行量化，如采用峰值起伏度[16]、水文弯曲度和表面分形幂律关系[17]或提出新的参数[18]来表征裂隙形貌，建立起非线性系数表征模型。但这些研究没有考虑裂隙的各向异性及渗流方向性的影响。

开度场的各向异性可显著影响裂隙渗流规律。剪切作用是比较常见的导致粗糙裂隙产生各向异性特征的方式。在各向异性粗糙裂隙中，裂隙渗透特性受各向异性程度、相对粗糙度等因素的共同影响[19-20]；在剪切粗糙裂隙中，垂直于剪切方向的渗透性往往大于平行于剪切方向的渗透性[21]。剪切裂隙各向异性程度与剪切位移、平均开度、粗糙度、法向应力等因素密切相关[22-24]；同时，剪切过程中也伴随着开度场的改变和通道流的形成与变化。目前对于各向异性粗糙裂隙渗流规律的研究成果多为定性认识，且对渗流的各向异性与开度、粗糙度等裂隙几何形貌特征的关系研究相对较少，缺乏对各向异性裂隙非线性渗流特性的量化表征，因此亟需深入研究。

本文采用地质统计随机模型生成表面符合分形特征的不同各向异性程度的粗糙裂隙，通过直接数值求解Navier-Stokes方程进行模拟，分析粗糙裂隙的非线性渗流特征，阐明裂隙的各向异性程度、裂隙开度、粗糙度等因素对非线性渗流的影响规律，提出考虑裂隙各向异性和粗糙度的非线性渗流参数表征模型，并推导出相应的临界雷诺数表达式，为裂隙网络及宏观尺度岩体渗流研究提供理论依据。

1 研究方法

1.1 粗糙裂隙生成

本文采用Brown[25]提出的数学模型来生成随机粗糙裂隙。在该模型中，通过分形维数D和标准差控制裂隙上、下表面的粗糙度，通过设置不匹配长度来控制上下面的不匹配程度及开度的空间分布。本文选取的分形维数为2.5，不匹配波数为8(不匹配长度=裂隙长度/不匹配波数)。开度均值和标准差的选取参照文献[17]试验中的数值。选择平均开度为0.5，1.0，1.5 mm；开度标准差为0.1，0.2，0.3。采用Brown方法生成的裂隙为随机裂隙，其中开度为0.5 mm和1.0 mm的裂隙是通过开度为1.5 mm的裂隙平移上表面得到的。

根据Brown模型，为考虑裂隙各向异性，将裂隙表面的功率谱密度G(k)表达为式(1)[26]：

G(k)∝[(kx/a)2+(ky/b)2]-(3.5-D)

(1)

式中：kx，ky为x，y方向上的波数，D为裂隙表面的分形维数，参数a和b表征裂隙的各向异性程度。

如图1所示，当b/a=1时，裂隙开度场为各向同性；当b/a<1时，裂隙开度场为各向异性，空间关联性或连通性在y方向上优于x方向，即y方向为开度场各向异性优势方向；当b/a>1时，x方向为优势方向。本文选取y方向为优势方向，各向异性程度表征参数b/a取值1，1/2，1/4用于研究裂隙表面的各向异性对渗流的影响。

本文中裂隙面的分辨率为512×512，生成裂隙面的大小为128 mm×128 mm，x-y方向网格大小为0.25 mm，截取其中120 mm×120 mm用于渗流计算。生成的随机粗糙裂隙的编号和几何参数如表1所列。图1为F1H3A1、F1H3A2、F1H3A3的开度分布图。

图1 不同各向异性程度b/a的开度分布Fig.1 Aperture field of different anisotropy factors b/a

表1 生成裂隙几何参数

1.2 控制方程

对于单个粗糙裂隙中不可压缩牛顿流体的定常流动，根据动量和质量守恒定律，流动过程可用Navier-Stokes方程和连续性方程来描述[27]。

ρ(u·∇u)-μ∇2u=-∇P

(2)

∇·u=0

(3)

式中：u为速度矢量，m/s；P为压力，Pa；μ为动力黏度，Pa·s；ρ为流体密度，kg/m3；等式(2)中，ρ(u·∇u)和μ∇2u分别代表流动过程中的惯性力和黏性力，式(3)是代表质量守恒的连续性方程。

本次模拟中入口处设定为流量边界条件，出口处设定为0压力边界条件，其余设置为无滑移边界。z方向划分网格7～10层。模型网格数量(2～4)×106。通过有限元软件求解Navier-Stokes方程[15]得到不同裂隙的流场。分别将流动方向设置为沿x方向和y方向，即垂直于优势方向和平行于优势方向，来研究裂隙渗流的各向异性。图2为沿x方向和沿y方向的两种流动情形的边界条件。

图2 流动边界条件Fig.2 Flow boundary conditions

Forchheimer方程常被用来宏观表征非线性流，具体如下：

-∇P=AQ+BQ2

(4)

(5)

(6)

式中：Q为流量，m3/s；A为描述黏性力能量损耗的系数，kg/(m5·s)；B为描述惯性力能量损耗的系数，kg/m8；w为裂隙垂直于流动方向上的宽度，m，在本文中为0.12 m；eh为等效水力开度，m；β为非达西系数或惯性系数，1/m，取决于裂隙几何结构；bD为无量纲化的系数。

本文中设eh为仅由几何形状控制的变量，沿用Fourar等[28]和Nowamooz等[29]的定义，通过线性阶段的系数A反算得到。

在裂隙渗流中，常用雷诺数Re和Forchheimer数Fo来表征流态的变化，定义为如下形式：

(7)

(8)

Re可以反映惯性力和黏性力的相对强弱，Fo可以反映流动偏离线性的程度。目前常使用非线性程度因子α来区分达西流和非达西流。

(9)

当α大于临界值时，可以认为水流已经发展为非达西流态，参照Zeng & Grigg[30]等的取值，临界值取0.1。

模拟的雷诺数范围为0.01～300。渗流数值模拟参数取值如表2所列。

表2 渗流数值模拟模型参数

1.3 裂隙粗糙度及各向异性的表征

目前裂隙粗糙程度表征方法可归纳为三大类：经验方法(JRC)、统计参数(平均值、均方根、方差等)和分形维数表征法[31]。本文采用一阶导数均方根Z2来描述裂隙的粗糙度[32]。

(10)

式中：n为采样点数，zi为第i个点的z坐标，Δd为采样点间隔。

对于三维曲面，可以在流动方向上取多条线计算Z2值，取其平均值作为该曲面的Z2值。根据生成裂隙表面的精度，沿x方向取512条剖面线，计算平均值作为垂直优势方向的Z2值，Z2，⊥。沿y方向取剖面线平均值为平行优势方向的Z2值，Z2，∥。

如表3所列，对于随机生成的裂隙，在b/a不变的情况下，其Z2值随均方差值增大而增大，但Z2，⊥/Z2，∥的值维持稳定，说明Z2可以较好地表示裂隙的粗糙度及各向异性。

随机生成多组不同各向异性程度的粗糙裂隙并计算其Z2，⊥/Z2，∥值，发现Z2，⊥/Z2，∥与裂隙的b/a值有很强的相关关系，拟合相关度大于0.99。Yin等的研究表明，Z2与分形维数和标准差之间也存在很好的相关性[33]。因此，可以用Z2来表征裂隙粗糙度及各向异性。

表3 不同开度均方差随机裂隙的一阶导数均方根

图3 Z2，⊥/Z2，∥与各向异性表征指标b/a的关系Fig.3 Relationship between Z2，⊥/Z2，∥ and anisotropy factor b/a

1.4 各向异性粗糙裂隙非线性系数表征模型

Chen等提出了双参数模型[34]，用峰值粗糙度和等效水力开度对式(6)中的无量纲非达西系数进行表征。根据不同粗糙度裂隙在法向应力下的流量-压力梯度实验数据，借助Levenberg-Marquardt优化算法，得到模型拟合系数的值：m=0.022，n=0.666。

(11)

(12)

式中：ξ为裂隙的峰值粗糙度，m。

对于各向异性裂隙，不同渗流方向的ξ值相同，因此上述表征模型不能很好地体现裂隙各向异性对非线性渗流的影响。本文参考上述双参数模型，用Z2和eh作为参数来表征非达西系数β，具体形式如下：

(13)

式中：κ，ψ，γ为拟合系数。

本文选取的表达式为经验表达式，不是严格的理论公式，选取该形式基于以下原因：① 形式上与双参数公式相似，用Z2代替ξ来表征裂隙形貌对于非线性渗流的影响。② 非达西系数与eh近似呈幂函数关系，在多个研究中已经得到证明[35]。③ 对于光滑的平板模型情况，有Z2=0，从而得到β=0。公式可以退化为光滑平板模型[8]。

2 数值模拟结果和分析

2.1 压力梯度-流量关系

根据数值模拟得到的结果，用Forchheimer公式对裂隙渗流的压力梯度和流量曲线进行拟合。如图4所示，不同粗糙度、机械开度和各向异性程度下裂隙渗流的∇P～Q曲线均与Forchheimer公式拟合良好(R2>0.99)，曲线呈现非线性特征，说明在选择的进口流量范围内，裂隙均已发展为非线性流动。

图4 压力梯度-流量关系Fig.4 Relationship between pressure gradient and flow rate for different fractures

从图4(a)、(b)可以看出，随着粗糙度增加和开度的减小(即相对粗糙度的增加)，∇P～Q曲线的非线性程度增强。从图4(c)可以看出，对于b/a为1的裂隙，垂直于优势方向(x方向)与平行于优势方向(y方向)的∇P～Q曲线基本一致。随b/a的减小，开度场各向异性增强，不同方向的∇P～Q曲线差别变大，且可以看出，垂直于优势方向的非线性程度更强。

2.2 表观渗透性

采用表观水力开度来反映裂隙渗流的表观渗透性，以此研究渗流的非达西特性。

(14)

式中：Ta为表观导水系数，m4；ea为表观水力开度，m。

图5给出了不同各向异性程度的裂隙水力开度ea和机械开度em之比随雷诺数的变化。可以看出，在雷诺数很小时(Re<1)，ea/em几乎保持不变，此时渗流处于达西状态，渗透率接近固有渗透率。随着雷诺数的不断增大，ea/em不断减小，裂隙渗透性不断降低。从图5可以看出，对于b/a为1的粗糙裂隙，垂直于优势方向与平行于优势方向的ea/em值基本一致，说明两个方向的渗透性大致相等。b/a<1时，不同方向的ea/em值差异较大，且平行于优势方向渗透性更强。

图5 水力开度与机械开度比值随雷诺数变化Fig.5 Relationship between ratio of hydraulic aperture to mean aperture and Reynolds number for fractures of different anisotropy factors

为了量化导水系数的各向异性，定义η为导水系数比：

(15)

式中：Ta，x，Ta，y为垂直于优势方向与平行于优势方向的表观导水系数。

图6为不同粗糙度、不同开度、不同各向异性程度的裂隙导水系数比随雷诺数的变化。从图中可以看出，在雷诺数较小时，η值基本不变。随着雷诺数增大，η值不断增大，随着水流惯性效应的增大，流动状态更加复杂，平行于流动方向的导水系数变化较小，因此η值不断增大。从图6(a)和(b)可以看出，随相对粗糙度的增大，裂隙内空间结构更加复杂，导水系数差异更大。从图6(c)可以看出，当b/a为1时，η值接近于1，说明垂直于优势方向与平行于优势方向的导水性在线性阶段和非线性阶段都基本相等，裂隙渗流具有各向同性。当b/a从1/2减小到1/4时，导水性各向异性程度增大。

图6 导水系数比η随雷诺数变化Fig.6 Relationship between ratio of hydraulic conductivity and Reynolds number for different fractures

2.3 临界雷诺数

非线性程度因子取0.1时可以得到流动的临界雷诺数。图7给出了临界雷诺数随着各向异性程度的变化。从图7中可以看出，在b/a为1的裂隙中，垂直于优势方向与平行于优势方向的临界雷诺数基本一致。随着b/a的减小，不同方向的临界雷诺数差异变大，且平行于优势方向的临界雷诺数大于垂直优势方向的临界雷诺数，即更不容易发生非线性流动。

图7 非达西程度因子α随雷诺数变化Fig.7 Relationship between non-Darcian flow factor and Reynolds number for fractures of different anisotropy factors

3 各向异性非线性系数表征模型验证

3.1 拟合系数的确定

根据数值模拟得到不同裂隙的∇P～Q数据，采用Forchheimer公式拟合，得到相应的系数A、B，根据式(5)、(6)反算得到eh和β。借助Levenberg-Marquardt优化算法，根据式(13)对模拟值进行回归分析，得到拟合系数κ=0.006 13，ψ=1.443 58，γ=1.558 9，相关系数R2=0.968，表明该模型可以很好表征各向异性粗糙裂隙渗流非线性系数。

注：图中数据点代表数值模拟值，曲面为模型拟合结果。图8 参数模型回归分析结果Fig.8 Regression analysis to obtain model parameters

将拟合结果代入式(13)、(6)中，可以得到各向异性粗糙裂隙非线性渗流的参数化模型：

(16)

(17)

图9为各向异性模型和双参数模型对粗糙裂隙F3H3A3的拟合结果。对于平行优势方向(y方向)的渗流，双参数模型和各向异性模型均拟合较好。对于垂直优势方向(x方向)的渗流，各向异性模型拟合较好，双参数模型有较大偏差。同时，可以观察到不同方向的双参数模型拟合结果较为接近，这是因为不同方向的峰值粗糙度相同，且等效水力开度相差不大。当裂隙形貌和流动方向确定时，各向异性模型能更好地拟合裂隙岩体中非线性渗流特征。