不同开度和接触情况下结构面剪切错动节理渗流分析

刘晓光

(盘锦河海土木工程咨询有限公司，辽宁 盘锦 124010)

天然岩体结构面的存在使得岩体的水力特性复杂多变，其中节理面的粗糙起伏对结构面的力学、渗透特性具有重要影响[1]。近年来，随着大型水利工程建设数量的不断增加，结构面强度和渗流问题成为工程设计和建设必须要面对和解决的课题。例如，法国在20世纪中期修建Malpasset拱坝过程中，由于没有充分考虑结构面渗透的影响，造成大坝设计标准过低，进而诱发了重大溃决事故[2]。由此可见，结构面力学和渗透特性研究对保证相关工程的建设和运行安全具有十分重要的意义和价值。

水利工程中的岩体渗流问题是一个十分复杂的问题，其中结构面形成的孔隙较天然岩石的其他孔隙小一个数量级，但是其造成的渗流却要大几个数量级[3- 5]。因此，结构面渗透特性研究是天然岩体渗透特性研究的关键与核心。在传统的岩土工程研究领域，均以线性达西定律对结构面渗透性进行评估。但是，近年来的诸多研究成果显示，结构面渗流过程中的流速与水力梯度之间并不是简单服从线性达西定律，而是呈现比较复杂的非线性关系[6]，基于传统的达西定律进行结构面渗流分析容易导致工程安全隐患[7]。天然岩体中的节理等结构面的剪切位移与错动变形十分常见，特别是受到施工扰动时，这种变形几乎是不可避免的。这种变形不仅会对工程的稳定性造成影响，同时与渗透性强烈相关[8]。在节理的剪切错动变形中，其开度、曲折率以及壁面接触均会发生变化，但是这些变化难以通过物理试验观测和界定，因此剪切错动变形如何影响节理渗流方面一直缺乏深入研究。在节理渗流的早期研究中，一般采取立方定律，但是受岩石中天然裂隙边界曲折的影响，该定律的研究结果往往存在较大误差。陈益峰等基于结构面的可压缩性，同时考虑节理裂隙面的粗糙度，提出了复杂应力条件下的广义立方定律，成为目前节理变形条件下渗流特征研究的主要理论[9]。

根据流体力学理论，节理流体运动过程可以通过N-S方程描述[10- 11]，但由于N-S方程是一组非线性偏微分方程，采用传统的流体计算动力学方法对其求解存在诸多问题。因此，本文通过COMSOL软件求解N-S方程，对节理剪切错动过程中的开度、曲折度以及接触面变化对非线性渗流的影响进行研究，以期进一步探究结构面剪切错动节理非线性渗流机理。

1 计算模型与边界条件

1.1 几何模型

研究采用的物理模型为耦合的上下两个节理面构成，其示意图如图1所示。其中，上节理面的长度为50mm，下节理面的长度为100mm，宽度均为25mm。其中，上节理面受到剪切荷载的作用，沿着下节理面做剪切滑动，节理面的粗糙系数在4～6之间。研究中选取节理面的上下面所对应的部分构建裂隙渗流几何模型。采取自由三角形网格对节理面进行网格剖分，最终获得1417002个计算单元。其中，最小单元尺寸为0.0256mm，最大单元尺寸为0.320mm。

图1 物理模型示意图

1.2 流体控制方程

根据李睿劬等学者的研究成果，黏性不可压流体在复杂裂隙中的稳定运动是流体质点运动的平衡过程。这种力的平衡关系可以通过数学形式反映在控制方程N-S方程[12]，其数学表达式如下：

(1)

式中，ρ—流体的密度；μ—流体的动力黏滞系数；U—流体质点的流速矢量；P—外加压力梯度；F—体力矢量。

不可压缩流体连续方程的表达式为：

U=0

(2)

由于本次研究对象为水体渗流，因此按照25℃条件下水的密度和动力黏滞系数取值，其中ρ为0.997g/cm3；μ为0.89×l0-3Pa·s。

2 节理错动方案与开度变化特征

本次模拟研究主要考虑上节理面剪切位移过程中产生的剪胀行为，根据几何模型的尺寸，设计了1、2、3、5mm四种不同的剪切步长，同时保证各个剪切步长条件下上下节理面的接触面积不大于1%。

节理的张开度是节理渗流的重要影响因素[13]，研究中对不同剪切状态下的上下节理面的张开度进行统计分析，结果显示节理的最大开度和最小开度分别为0.53、0mm，平均开度为0.213mm，标准差为0.056。进一步对节理开度进行频率统计分析，获得如图2所示的不同开度频率分布直方图。由图2可知，节理绝大部分位置的开度介于0.15～0.30mm之间，0.10mm以下的小开度和0.40mm以上的大开度所占的比例极小。通过分析开度频率分布特征，发现Gaussian曲线可以对上述特征进行较好的拟合，其表达式为：

(3)

式中，f(x)—节理面开度统计频率；u—平均开度；σ—标准差。

图2 节理开度分布直方图

按照上述方法对四种不同剪切步长下的开度进行统计计算，获得平均开度与标准差计算结果，见表1。从表1中的计算结果可知，节理的开度会随着剪切位移的增加而增大，但是接触面积的变化并不明显，呈现出比较稳定的特征。从开度的标准差计算结果来看，随着节理面剪切位移的不断增大，标准差由0.056逐渐增加到0.476，说明开度分布离散程度不断增加。

表1 不同剪切步长节理开度特征

3 模型计算结果与分析

3.1 主要渗流特征量分析

为了对节理渗流的非线性特征进行模拟，研究中选择雷诺数变化范围为5～30。利用数值模拟的方法对不同剪切步长和雷诺系数条件下的节理顺流向流速进行计算。结果显示，当剪切步长为0mm时，节理面与接触面上的顺流方向流速均为0，但是在开度较大的局部流速分量较大，而局部开度变化比较明显的部位，顺流向流速分量较小，究其原因主要是这些部位的渗流通道较为曲折。此外，随着模型计算中入口流量的增大，裂隙内顺流方向的流速分量由正变负，并不断减小，说明渗流通道内存在局部回流。究其原因，主要是节理裂隙内存在局部漩涡，并且随着雷诺数的增加，裂隙内的漩涡也逐渐增大，因而造成逐渐增加的惯性力损失。由此可见，节理裂隙内的渗流并不服从线性达西定律。

图3 雷诺数为10时不同位移步长节理流线分布图

在不同剪切步长条件下，对节理裂隙中的渗流路径进行模拟研究，获得如图3所示的雷诺数为10时不同位移步长节理流线分布图。由图3可知，节理裂隙中的流速不均匀分布特征十分明显，流线也呈现出强烈的曲折性特点，究其原因主要是开度分布不均。另一方面，随着剪切步长的增加，节理面的接触率随之增大，并直接导致了沟槽优势流的形成，沟槽流效应对节理面渗流有増强作用，该模拟结果与相关学者的研究结果一致[14- 15]。

图4 雷诺数为10时不同位移步长节理渗流压力分布图

对不同位移步长条件下的节理渗流压力分布进行模拟，如图4所示。结果显示，在位移步长为0mm条件下，压力分布的主要特征是模型进口到出口逐渐减小的特征，存在比较明显的压力损失。此外，在压力的空间分布上看，模型进口部位损失较快而出口部位损失较慢，呈现出明显的空间分布不均匀性。其他位移步长条件下的压力分布具有相似的特点，但是压力分布的不均匀性随着位移步长的增加逐渐增大。究其原因，主要是节理面的接触使渗流路径受到阻挡，进而产生愈加明显的绕流效应，这也成为接触面附近压力降低的重要原因。

表2 拟合线性和非线性参数以及临界雷诺数的计算结果

3.2 非线性特征分析

为了进一步研究节理剪切错动条件下的非线性渗流特点，研究中利用模型模拟的方法对5种剪切步长下的渗流模型进行25次数值模拟，根据模拟结果对模型的出口压力梯度与流量之间的关系进行统计。试验结果显示Forchheimer方程能够有效地描述这一渗流过程，拟合结果如图5和图6所示。由图5和图6可以看出结构面剪切错动节理渗流具有显著的非线性特征。

图5 位移步长0mm条件下流量水力梯度关系曲线

图6 其余位移步长条件下流量水力梯度关系曲线

对每一个剪切步长，采用Forchheimer方程拟合实验数据，得到Forchheimer方程线性参数a与非线性参数b，结果见表2。从表2中的结果可知，随着位移步长的增加，线性参数a与非线性参数b均呈现出减小的趋势，并且减小幅度均在两个数量级左右，同时两个系数的减小呈现出先快后慢的特点。这说明随着位移步长的增加，节理张开度的变化对渗透性的影响逐渐减小。根据Zeng和Grigg的研究成果[13]，在每一个剪切步长计算节理中线性与非线性渗流的临界雷诺数Rec的计算方法见公式(4)，计算结果见表2。

(4)

由表2中的结果可知，线性参数a的变化范围为5.10×108～5.37×1010kg·s-1·m-5；非线性参数b的变化范围为1.80×1014～1.23×1016kg·m-8；临界雷诺数的变化范围为13.5～19.2。由此可见，在将来的工作中，建立一个能够考虑节理形貌的临界雷诺数模型对于界定节理非线性流有着重要的工程意义。

4 非线性渗流机理与讨论

根据上文的模拟计算结果，在节理剪切错动稳定渗流状态下，节理开度较小时体力可以忽略不计，流速不会随着时间的变化而变化。因此，N-S方程可以进一步简化，简化后的表达式如下：

P=μ2U+ρU·U

(5)

由此可以看出，水力梯度损失主要包括线性项和非线性项两部分。其中，线性损失主要为黏性力损失；非线性损失主要为惯性力损失。根据流体力学理论，渗流过程中的惯性力损失大小主要取决于流体质点的运动速度。因此，在流速较小的情况下，此时节理面中流体的流动变为蠕动流，相应地，惯性力损失影响变得十分微小，渗流的线性特征就比较显著。另一方面，非线性惯性力的微分表达式可以进一步表达为如下形式：

U·U=U·(ux，uy，uz)

(6)

显然，惯性力的大小主要受流速U及其空间位置变化率U两个因素的影响。研究中对雷诺数为30，不同剪切位移步长下的节理内流速进行模拟计算，并对计算结果进行平均。结果显示，当位移步长为0、5mm时，平均流速分别为0.119m/s和0.022m/s。由此可见，当位移步长增加时平均流速会大幅减小，说明U不变时，模型节理面剪切位移量的增加有助于减小渗流流体的惯性力损失。

另一方面，针对雷诺数为30，不同剪切位移步长下的节理内流线分布进行模拟计算分析，如图7所示。结果显示，随着位移步长的增加，节理内流线变得更为曲折，说明惯性力在流速不变的情况下，位移量的增加可以导致惯性力空间位置变化率的增大，对增加惯性力的损失较为有利。

综合上述，天然岩石中的节理开度均匀，那么节理渗流的非线性特征主要是由流速增大惯性力损失造成的；而节理内部的流速较低时，节理渗流的非线性特征主要是由渗流路径曲折造成的。

图7 雷诺数为30时不同位移步长的节理渗流流线分布图

对上下节理面接触部位的纵截面沿x方向的流速分量进行计算。结果显示：该方向上的流速有正有负，并分别代表正向渗流和反向回流。从具体计算结果来看，位移步长为0、5mm情况下接触部位的最小流速分别为-9.14×10-3、-3.34×10-4m/s。说明节理内部的涡旋与流速之间存在比较密切的关系，流速越大，越容易引发回流涡旋，从而造成惯性力损失，进而导致节理渗流的非线性特征出现。

5 结论

(1)节理裂隙中的流速不均匀分布特征十分明显，流线也呈现出强烈的曲折性特点，导致节理渗流存在着比较明显的沟槽流和优势流效应。

(2)随着雷诺数的增加，节理渗流逐渐表现出显著的非线性特征，流量和水力梯度之间的关系可以用Forchheimer方程拟合。在剪切位移量从0mm增加到5mm的过程中，线性参数a与非线性参数b均减小了两个数量级左右。

(3)节理渗流的非线性特征主要是由流速增大惯性力损失造成的；当节理内部的流速较低时，节理渗流的非线性特征主要是由渗流路径曲折造成的。