1、第 卷 第 期 年 月人 民 长 江 ，收稿日期：基金项目：国家自然科学基金项目（，）作者简介：马亚楠，女，硕士研究生，主要从事粗糙裂隙渗流方面的研究。：通信作者：杨志兵，男，教授，博士生导师，博士，主要从事复杂地下介质渗流机理与控制方面的研究。：文章编号：（）引用本文：马亚楠，杨志兵，熊小锋，等 粗糙裂隙各向异性对非达西渗流特性的影响 人民长江，（）：，粗糙裂隙各向异性对非达西渗流特性的影响马 亚 楠，杨 志 兵，熊 小 锋，陈 益 峰，（武汉大学 水资源与水电工程科学国家重点实验室，湖北 武汉；武汉大学 水工岩石力学教育部重点实验室，湖北 武汉）摘要：为研究各向异性粗糙裂隙中的非线性流动规
2、律，生成了符合岩石裂隙表面几何特征和开度分布的各向异性程度粗糙裂隙，用数值方法求解 方程来模拟粗糙裂隙渗流。采用 公式拟合数值模拟得到的压力梯度与流量关系，分析粗糙度、开度和各向异性程度对裂隙非达西渗流参数的影响规律。研究发现表观导水度各向异性随相对粗糙度和开度场各向异性的增大而增大。最后提出了基于一阶均方根粗糙度和等效水力开度的非达西系数表征模型，进一步推导了相应的临界雷诺数公式。所提出的裂隙非达西系数表征模型可以很好地描述各向异性粗糙裂隙非线性渗流特征，并与室内试验结果拟合良好。关 键 词：粗糙裂隙；非线性渗流；各向异性；方程；临界雷诺数中图法分类号：文献标志码：引 言在高坝枢纽、高水头抽
3、水蓄能电站、煤矿开采和深埋隧洞等工程活动中，受渗透水压和裂隙形貌的影响，岩体渗流具有强烈的非线性，发生渗漏、突透水和渗透破坏的风险急剧增大，因此非线性渗流对于工程建设和运行稳定有着重要影响。岩体渗流主要是通过裂隙结构面进行的。受地质沉积环境、温度变化和地质构造作用等因素的影响，天然裂隙面存在各向异性特征，使得裂隙渗流模拟及预测存在不确定性。研究各向异性粗糙裂隙中的非线性渗流规律，对于预测宏观尺度裂隙岩体渗流，预防隧洞涌水、内水外渗、煤矿突水和瓦斯突出等事故具有重要理论意义和参考价值。近年来已有许多学者开展了关于裂隙岩体渗流的研究。基于平行板假设的立方定律被广泛用于描述裂隙介质中的渗流。以往的研
4、究针对天然裂隙粗糙特性的影响，提出了一系列修正的立方定律。偏离线性规律的非达西渗流现象在数值模拟、室内外试验和工程应用中广泛存在。研究发现裂隙渗流的非线性受到法向应力、剪切错动、接触面积等因素的影响，其本质上是裂隙几何形状和开度的变化对裂隙渗流的控制作用。非达西渗流细观机制和宏观表征的研究得到了普遍关注。在细观机制方面，理论分析结果显示当惯性效应占主导作用时渗流出现非达西特征；此外，数值模拟结果表面涡旋区的产生和发展是非达西渗流发生的主要原因。在宏观表征方面，公式和 公式常用于描述非达西渗流的压力梯度与流量关系。其中 公式物理意义明确，参数取值稳定性强，被广泛用于裂隙介质非线性流的表征中。目前
5、主要通过建立裂隙形貌与非达西渗流公式参数的关系来对非线性系数进行量化，如采用峰值起伏度、水文弯曲度和表面分形 人 民 长 江 年幂律关系或提出新的参数来表征裂隙形貌，建立起非线性系数表征模型。但这些研究没有考虑裂隙的各向异性及渗流方向性的影响。开度场的各向异性可显著影响裂隙渗流规律。剪切作用是比较常见的导致粗糙裂隙产生各向异性特征的方式。在各向异性粗糙裂隙中，裂隙渗透特性受各向异性程度、相对粗糙度等因素的共同影响；在剪切粗糙裂隙中，垂直于剪切方向的渗透性往往大于平行于剪切方向的渗透性。剪切裂隙各向异性程度与剪切位移、平均开度、粗糙度、法向应力等因素密切相关；同时，剪切过程中也伴随着开度场的改变
6、和通道流的形成与变化。目前对于各向异性粗糙裂隙渗流规律的研究成果多为定性认识，且对渗流的各向异性与开度、粗糙度等裂隙几何形貌特征的关系研究相对较少，缺乏对各向异性裂隙非线性渗流特性的量化表征，因此亟需深入研究。本文采用地质统计随机模型生成表面符合分形特征的不同各向异性程度的粗糙裂隙，通过直接数值求解 方程进行模拟，分析粗糙裂隙的非线性渗流特征，阐明裂隙的各向异性程度、裂隙开度、粗糙度等因素对非线性渗流的影响规律，提出考虑裂隙各向异性和粗糙度的非线性渗流参数表征模型，并推导出相应的临界雷诺数表达式，为裂隙网络及宏观尺度岩体渗流研究提供理论依据。研究方法 粗糙裂隙生成本文采用 提出的数学模型来生成
7、随机粗糙裂隙。在该模型中，通过分形维数 和标准差控制裂隙上、下表面的粗糙度，通过设置不匹配长度来控制上下面的不匹配程度及开度的空间分布。本文选取的分形维数为 ，不匹配波数为（不匹配长度 裂隙长度 不匹配波数）。开度均值和标准差的选取参照文献 试验中的数值。选择平均开度为 ，；开度标准差为，。采用 方法生成的裂隙为随机裂隙，其中开度为 和 的裂隙是通过开度为 的裂隙平移上表面得到的。根据 模型，为考虑裂隙各向异性，将裂隙表面的功率谱密度（）表达为式（）：（）（）（）（）（）式中：，为，方向上的波数，为裂隙表面的分形维数，参数 和 表征裂隙的各向异性程度。如图 所示，当 时，裂隙开度场为各向同性；
8、当 时，裂隙开度场为各向异性，空间关联性或连通性在 方向上优于 方向，即 方向为开度场各向异性优势方向；当 时，方向为优势方向。本文选取 方向为优势方向，各向异性程度表征参数 取值，用于研究裂隙表面的各向异性对渗流的影响。本文中裂隙面的分辨率为 ，生成裂隙面的大小为 ，方向网格大小为 ，截取其中 用于渗流计算。生成的随机粗糙裂隙的编号和几何参数如表 所列。图 为、的开度分布图。图 不同各向异性程度 的开度分布 表 生成裂隙几何参数 裂隙编号开度标准差 开度均值 起伏度峰值 注：、分别代表开度标准差为，；、分别代表开度均值为，；、分别代表 为，。第 期 马亚楠，等：粗糙裂隙各向异性对非达西渗流特
9、性的影响 控制方程对于单个粗糙裂隙中不可压缩牛顿流体的定常流动，根据动量和质量守恒定律，流动过程可用 方程和连续性方程来描述。()（）（）式中：为速度矢量，；为压力，；为动力黏度，；为 流 体 密 度，；等 式（）中，（）和 分别代表流动过程中的惯性力和黏性力，式（）是代表质量守恒的连续性方程。本次模拟中入口处设定为流量边界条件，出口处设定为 压力边界条件，其余设置为无滑移边界。方向划分网格 层。模型网格数量（）。通过有限元软件求解 方程得到不同裂隙的流场。分别将流动方向设置为沿 方向和 方向，即垂直于优势方向和平行于优势方向，来研究裂隙渗流的各向异性。图 为沿 方向和沿 方向的两种流动情形的
10、边界条件。图 流动边界条件 方程常被用来宏观表征非线性流，具体如下：（）（）（）式中：为流量，；为描述黏性力能量损耗的系数，（）；为描述惯性力能量损耗的系数，；为裂隙垂直于流动方向上的宽度，在本文中为 ；为等效水力开度，；为非达西系数或惯性系数，取决于裂隙几何结构；为无量纲化的系数。本文中设 为仅由几何形状控制的变量，沿用 等和 等的定义，通过线性阶段的系数 反算得到。在裂隙渗流中，常用雷诺数 和 数 来表征流态的变化，定义为如下形式：（）（）可以反映惯性力和黏性力的相对强弱，可以反映流动偏离线性的程度。目前常使用非线性程度因子 来区分达西流和非达西流。（）当 大于临界值时，可以认为水流已经发
11、展为非达西流态，参照 等的取值，临界值取。模拟的雷诺数范围为 。渗流数值模拟参数取值如表 所列。表 渗流数值模拟模型参数 参数单位参数值水的密度 水的黏滞系数 入口流量 ，出口压力 雷诺数，裂隙粗糙度及各向异性的表征目前裂隙粗糙程度表征方法可归纳为三大类：经验方法（）、统计参数（平均值、均方根、方差等）和分形维数表征法。本文采用一阶导数均方根 来描述裂隙的粗糙度。()（）式中：为采样点数，为第 个点的 坐标，为采样点间隔。对于三维曲面，可以在流动方向上取多条线计算值，取其平均值作为该曲面的 值。根据生成裂隙表面的精度，沿 方向取 条剖面线，计算平均值作为垂直优势方向的 值，。沿 方向取剖面线平
12、均值为平行优势方向的 值，。如表 所列，对于随机生成的裂隙，在 不变的情况下，其 值随均方差值增大而增大，但，的值维持稳定，说明 可以较好地表示裂隙的粗糙度及各向异性。随机生成多组不同各向异性程度的粗糙裂隙并计算其，值，发现，与裂隙的 值有很强的相关关系，拟合相关度大于 。等的研究表 人 民 长 江 年明，与分形维数和标准差之间也存在很好的相关性。因此，可以用 来表征裂隙粗糙度及各向异性。表 不同开度均方差随机裂隙的一阶导数均方根 开度均方差，图 ，与各向异性表征指标 的关系 ，各向异性粗糙裂隙非线性系数表征模型 等提出了双参数模型，用峰值粗糙度和等效水力开度对式（）中的无量纲非达西系数进行表
13、征。根据不同粗糙度裂隙在法向应力下的流量 压力梯度实验数据，借助 优化算法，得到模型拟合系数的值：，。()（）（）式中：为裂隙的峰值粗糙度，。对于各向异性裂隙，不同渗流方向的 值相同，因此上述表征模型不能很好地体现裂隙各向异性对非线性渗流的影响。本文参考上述双参数模型，用 和 作为参数来表征非达西系数，具体形式如下：（）式中：，为拟合系数。本文选取的表达式为经验表达式，不是严格的理论公式，选取该形式基于以下原因：形式上与双参数公式相似，用 代替 来表征裂隙形貌对于非线性渗流的影响。非达西系数与 近似呈幂函数关系，在多个研究中已经得到证明。对于光滑的平板模型情况，有，从而得到 。公式可以退化为光
14、滑平板模型。数值模拟结果和分析 压力梯度 流量关系根据数值模拟得到的结果，用 公式对裂隙渗流的压力梯度和流量曲线进行拟合。如图 所示，不同粗糙度、机械开度和各向异性程度下裂隙渗流的 曲线均与 公式拟合良好（），曲线呈现非线性特征，说明在选择的进口流量范围内，裂隙均已发展为非线性流动。图 压力梯度 流量关系 从图（）、（）可以看出，随着粗糙度增加和开度的减小（即相对粗糙度的增加），曲线的非线性程度增强。从图（）可以看出，对于 为 的裂隙，垂直于优势方向（方向）与平行于优势方向（方向）的 曲线基本一致。随 的减小，开度场各向异性增强，不同方向的 曲线差别变大，且可以看出，垂直于优势方向的非线性程度
15、更强。第 期 马亚楠，等：粗糙裂隙各向异性对非达西渗流特性的影响 表观渗透性采用表观水力开度来反映裂隙渗流的表观渗透性，以此研究渗流的非达西特性。（）式中：为表观导水系数，；为表观水力开度，。图 给出了不同各向异性程度的裂隙水力开度 和机械开度 之比随雷诺数的变化。可以看出，在雷诺数很小时（），几乎保持不变，此时渗流处于达西状态，渗透率接近固有渗透率。随着雷诺数的不断增大，不断减小，裂隙渗透性不断降低。从图 可以看出，对于 为 的粗糙裂隙，垂直于优势方向与平行于优势方向的 值基本一致，说明两个方向的渗透性大致相等。时，不同方向的 值差异较大，且平行于优势方向渗透性更强。图 水力开度与机械开度比
16、值随雷诺数变化 为了量化导水系数的各向异性，定义 为导水系数比：，（）式中：，为垂直于优势方向与平行于优势方向的表观导水系数。图 为不同粗糙度、不同开度、不同各向异性程度的裂隙导水系数比随雷诺数的变化。从图中可以看出，在雷诺数较小时，值基本不变。随着雷诺数增大，值不断增大，随着水流惯性效应的增大，流动状态更加复杂，平行于流动方向的导水系数变化较小，因此 值不断增大。从图（）和（）可以看出，随相对粗糙度的增大，裂隙内空间结构更加复杂，导水系数差异更大。从图（）可以看出，当 为 时，值接近于，说明垂直于优势方向与平行于优势方向的导水性在线性阶段和非线性阶段都基本相等，裂隙渗流具有各向同性。当 从 减小到 时，导水性各向异性程度增大。图 导水系数比 随雷诺数变化 临界雷诺数非线性程度因子取 时可以得到流动的临界雷诺数。图 给出了临界雷诺数随着各向异性程度的变化。从图 中可以看出，在 为 的裂隙中，垂直于优势方向与平行于优势方向的临界雷诺数基本一致。随着 的减小，不同方向的临界雷诺数差异变大，且平行于优势方向的临界雷诺数大于垂直优势方向的临界雷诺数，即更不容易发生非线性流动。图 非达西程度因子
