1.2 单裂隙中水流及溶质运移问题的研究进展

1.2.1 单裂隙中水流问题研究进展

1.2.1.1 单裂隙水流理论模型的发展

张开度远远小于其长宽尺寸的断裂称为裂隙。关于裂隙中水流问题的研究已经开展了几十年［21^，22］，这些研究始终集中在四个基本方面［23］：①建立裂隙概念模型；②建立模型解析解或数值解的方法；③对裂隙水力特征进行描述；④运用随机方法描述裂隙隙宽及水文地质参数。早期研究单裂隙时，通常把其简化成平直、光滑且隙宽处处相等，通常忽略裂隙基质本身的渗透性，因此在整个裂隙面上的渗透性是各向同性。等温条件下不可压缩牛顿流体Navier-Stokes方程如下

式中 ρ——流体密度；

μ——液体动力黏滞系数；

u_i——流体在i方向上的速度分量；

P——流体压力；

g_i——方向i上的重力加速度分量。

对于稳定流，式（1.1）中左边项为0，于是得到

式(1.2)也被称为Stokes方程，在平行板模型中，该方程的解如下

式中 b——裂隙宽度。

从式(1.3)中可以看出速度剖面呈抛物线形态，如图1.2所示。通常被称为泊肃叶流(Poiseuille's flow)。在宽度上对其积分可得

对其整理即可得到著名的裂隙水流立方定律

式中 q——单宽渗流量；

J——水力坡降；

b——裂隙宽度；

g——重力加速度；

v——水流的运动黏滞系数。

运用立方定律求得的裂隙内部流场为层流且水流速度处处相等，而在紊流阶段或微裂隙中运用立方定律求得的结果往往不尽如人意［24］。

理想裂隙在自然界是不存在的，天然裂隙面均为粗糙裂隙，其隙宽是沿程变化的。通常上下裂隙面会有一定程度的接触，而且其有效隙宽往往取决于隙面所受法向应力的大小。天然裂隙中水流在流动过程中遇到许多阻碍，即使对于裂隙中的纯粹泊肃叶流，其与裂隙面接触的占裂隙体积10%的区域只传递约5%的流量，Djik和Berkowitz［25］运用核磁共振技术，获得饱和粗糙裂隙中水流图像，其中水流速度剖面呈类似抛物线，但不完全对称，这种速度分布对其中水流的影响十分明显，因此人们开始怀疑立方定律是否一直有效。Hakami和Larsson［26］对瑞典Äspö地区采集的天然花岗岩裂隙进行了水力实验和隙宽测量，结果表明测量得到的平均隙宽是通过立方定律得到的水力隙宽的约1.4倍。随着研究的深入，人们逐渐认识到把裂隙简单地简化为平直光滑的平板模型还具有很多其他方面的局限性。例如，裂隙中局部水流密度和温度的改变都将对流场产生影响，而立方定律只考虑了裂隙的隙宽、尺寸和水的黏滞系数(动力黏滞系数或运动黏滞系数)，虽然流体的黏滞系数和温度、密度密切相关，但是由温度/密度的不均导致的热对流/密度流却被忽略。除了温度、密度之外，裂隙压力也对其中的流场产生一定的影响。Witherspoon和Wang通过室内实验证明，当平板裂隙压力超过10MPa时，用立方定律得出的结果和实际情况相差甚远［27］。这说明只有在一定的压力范围内才能运用立方定律进行简化运算。

考虑到立方定律的局限性，众多学者对其进行了改进，该研究主要是为了揭示和定量分析传统立方定律在实际运用中的偏离情况。而改进的关键问题则在于隙宽的确定。一般来说，隙宽主要有三种定义［28］，分别是平均隙宽<b>、机械隙宽b_m和等效水力隙宽b_h。平均隙宽<b>是指裂隙隙宽函数b(x,y)的平均值；机械隙宽b_m为裂隙面发生的最大闭合变形量；等效水力隙宽b_h是为了应用立方定律于实际裂隙而提出的概念，即是将实验所得裂隙渗流量代入立方定律反求得到的裂隙宽度。对于光滑平行板裂隙，这三种隙宽值是相等的;而对于实际粗糙裂隙，它们通常是不等的。

Lomize通过一系列实验研究了粗糙裂隙对水流的影响，在此基础上引进了裂隙粗糙度的概念［1］。Snow将此粗糙度概念运用于天然裂隙的水流模拟中，并对传统立方定律进行了相应的评价［29］。Lomize和Louis在大量实验基础上提出的立方定律修正公式，见表1.1和表1.2。

表1.1和表1.2中为紊流时裂隙渗透系数，V_f为裂隙水流平均速度，∆为裂隙绝对糙率。由于天然裂隙中粗糙颗粒分布不均且凸起高度差异较大，因此在实际情况中仍无法确定∆的值。紊流时，水头损失与流速呈非线性关系，可用V^m=-表示，式中m为紊流时的非线性指数，其变化范围为1～2[2]。光滑裂隙中紊流公式为q^m==，由此可得

通过室内实验可以得到lgJ和lgq之间的关系直线，其斜率即为m。

Barton等首次提出将岩石节理粗糙系数JRC运用于裂隙水力隙宽的求解，提出水力隙宽、机械隙宽和JRC之间的关系式［31］

于是立方定律修正为

上式中的难点在于JRC的确定，虽然很多学者提出一系列JRC的测量方法［32^-35］，然而由于其尺度效应［36］及粗糙度的各向异性等条件限制，使得测量准确度仍有待提高。

Iwai通过大量单裂隙水流实验发现裂隙的面积接触率ω(裂隙面接触面积与总面积之比)与水流规律存在一定的联系［37］。据此得出的理想裂隙模型中，沿水流法向方向，裂隙开度连续变化；而在沿水流方向，裂隙开度不变。进而Walsh［38］和周创兵等［39］分别推导出如下公式

式中 κ——经验系数，Walsh建议κ=1，周创兵等建议κ=0。

Amadei等则提出如下的修正公式［40］

式中 σ_b——隙宽均方差。

Nolte对等石英二长岩芯试样进行渗流实验,3个试样的成果整理发现，隙宽指数n远远大于3，随隙宽增大分别为7.6、8.3和9.8［41］。张有天等采用计算机生成人工裂隙和有限元数值计算方法，经分析认为，单宽流量q与b_m也不是3 次方关系，n均大于3，并随裂隙粗糙程度的增加而增大［42］。耿克勤根据人工、天然光滑和粗糙裂隙的实验结果分析得到，对于小开度裂隙层流而言，1.7≤n≤3.0，裂隙面几何形态越光滑n值越大;对于中开度过渡状态，0.8≤n≤1.4;对于大开度裂隙，0.3≤n≤0.48［43］。在此基础上，许光祥等进行了归纳总结，提出了次立方定律和超立方定律的概念，并通过试验进行了验证［44］。

Neuzil和Tracy模拟了由一系列不同隙宽的平行板模型组成的裂隙，这些平行板模型的隙宽分布符合对数正态分布。模拟结果显示，此裂隙中的水流主要通过隙宽较大的路径，裂隙中水流存在优势流［45］。

Tsang等利用电阻代替水流发现裂隙的曲折程度对其中的水流有着较为明显的影响，而且隙宽越小，曲折度的影响越明显。和普通的光滑平行板模型相比，裂隙曲折度和隙面粗糙度对水流产生的阻滞影响可以将水流速度降低3个数量级以上［46］。Brown详细研究了裂隙粗糙度对水流的影响［47］，他利用分形方法生成的粗糙裂隙，模拟了其中的水流状态。他认为，在隙宽较小的情况下，裂隙的曲折度对其中的水流有着十分明显的影响，这也进一步证明了Tsang的结论。

Rasmason和Neretnieks历时3年在瑞典一个废弃的Stripa矿里进行了现场实验［48］，通常称这个实验为Stripa-3D实验，它揭示了裂隙岩体中水流和运移的主要特征。Heath［48］和Bourke［49］也在康沃尔做了现场实验，他们的研究发现，水流很大程度上在裂隙岩体中被隔离的管道流里面流动，并没有大面积的水流发生，水流仅仅在5%~20%的裂隙面内流动。在加拿大安大略的白垩纪石灰岩和美国伊利诺伊州东北部的志留纪白云岩中，Novakowski［50］、Raven［51］和Shapiro［52］分别进行了现场示踪剂实验，发现通过裂隙隙宽估计的水力传导系数和溶质的迁移规律与现场观测的不一致，如果裂隙隙宽的变化很小，则与实际观测吻合得较好。因此，上述的现场实验提供给我们这样一个事实：裂隙面的平行板模型假定并不适用，水流和溶质运移只在裂隙形成的管道或沟槽中发生。

Pyrak［53］等在不同压力条件下，将熔化的金属代替水流注入裂隙中，冷却后打开，发现液体在其中的流动路径是曲折的，沿一定的路径流动，具有优势流和沟槽流的特征。为了有效地模拟裂隙中的沟槽流，Y.W.Tsang和C.F.Tsang［46^，54］提出了一种新的概念模型。他们运用对数正态分布函数并在平均、差异和空间相关长度基础上生成隙宽统计分布，进行了水流和溶质运移数值模拟，预测了其中的穿透曲线，预测的结果和Moreno等［55］提供的数据能很好地吻合。

1.2.1.2 单裂隙物理模型实验

由于现场水文地质条件以及岩体的复杂性，开展现场实验工作需要大量的人力和物力，测试方法和测试装置的不同，所得到的结果也不一定完全一致［56］。室内实验相对来说花费较少而且有较精确的测试方法，因此大部分实验通过室内的小尺度模型进行。

研究者最初是利用两块光滑平板，通过控制两块平行板之间的距离来模拟不同隙宽的裂隙，这类装置较为简单，使用的平行板通常为有机玻璃板或玻璃板，国内外众多学者都做过类似的实验［1^{，24，30］}，图1.3为速宝玉等制作的光滑裂隙水流实验装置。在此基础上进一步发展，人们开始在平行板上粘贴粗糙颗粒来模仿粗糙裂隙，大多数此类实验都是利用标准粒径的砂粒来充当粗糙颗粒［57^-59］，并通过粘贴不同粒径的砂粒和控制平行板之间的距离来获得不同的粗糙裂隙。也有研究者利用机械加工后的钢板来模拟粗糙裂隙。许光祥等［44］为了验证其理论，加工了包括光滑裂隙在内的5种不同形态的裂隙进行水流实验，如图1.4所示。由于粘贴的砂粒粒径或机械加工出的粗糙颗粒都比较均匀，虽然获得了粗糙的裂隙，但和实际情况仍相距甚远。于是研究者开始利用混凝土来获得仿真裂隙［56^，60-62］并进行水流及溶质运移实验(图1.5)。这在一定程度上让研究者获得更为接近实际的实验数据，但是混凝土和天然岩体在材料特性上存在较大差别，如透水性、表面可湿性等，都对其内部水流及裂隙本身形态产生影响，于是一些学者开始利用岩芯中的天然裂隙［63］，或者利用天然岩样制作人工裂隙，如图1.6所示。也有直接在野外进行较大尺度的单裂隙水流实验［64］。在过去的几十年中，研究人员逐渐开始运用诸如裂隙粗糙面剖面测量、低熔点金属注入和树脂铸造等技术来获得更加接近真实裂隙的模型来进行相关的模拟［46^{，47，53，65］}，这些模型考虑了裂隙隙宽变化及其对水流影响的实际情况，对研究变隙宽裂隙中的水流运动提供了极大的帮助。

1.2.1.3 裂隙水流问题的数值模拟

对于某些复杂问题，如在寻求裂隙的渗透性质与影响溶质运移的控制因素的关系时，现场实验和室内实验的测试手段都不完美，这就需要求助于数值模拟，对于这些复杂问题，数值模拟是一种比较方便且精度较高的方法。一些室内实验无法模拟的现象通过数值模拟得到了较好的解决［56］。

自20世纪60年代末期以来，关于裂隙中水流等相关问题的数值模拟已经发展了40多年［66］。据统计，截至1994年，已经出现至少30余种求解此类相关问题的专门软件［67］。Streltsova-Adams［68］在求解混合含水层井流问题时描述了几种求解裂隙水流的解析法。Elsworth［69］提出了求解具有特定几何形状层流或紊流问题的解析方法，由于对问题中水流的形状有特定要求，因此其应用受到了限制。Amadei和Illangasekare［70］运用积分变换得到矩形裂隙中流体势能场的连续表达式。这个方法可以在裂隙隙宽和粗糙度各向异性的情况下使用，同时也可以为进一步研究立方定律提供依据。由于积分变换的运用，不需要再对裂隙进行离散，但是裂隙的几何形态不能过于复杂，运用仅限于裂隙形态较为简单的情况。

在求解裂隙中水流及传质问题时，通常用到差分和积分法来求解平衡方程。对于空间求导问题，积分法往往获得比有限差分法(FDM)更为广泛的应用，原因之一就是积分法能更好地适应不规则几何区域。常用的积分法包括有限单元法(FEM)［70^{，71，72］}和边界元法(BEM)。Elsworth［69^，73］曾运用BEM-FEM混合方法来模拟裂隙水流问题，同时积分有限差分法(IFDM)也获得了广泛的运用［74^-76］。积分型有限差分法通常使用Crank-Nicholson法逼近来处理时间导数。为了在高Péclet数条件下获得稳定解，Sudicky和McLaren［77］引进了Laplace变换。

20世纪80年代，在格子气自动机(Lattice Gas Automata)理论的基础上发展出了一种新型流场模拟方法——格子波茨曼方法(Lattice Boltzmann Method)［78^-84］。与传统的模拟方法相比，该方法具有规则简单、复杂边界易处理、能适应大规模并行计算等优点，并迅速成为模拟复杂流场的新工具［85^-92］。

早在20世纪六七十年代，研究人员就已经开始利用诸如时间序列法、谱分析法和蒙特卡罗法来模拟裂隙隙宽概率分布，这些方法的共同点就是通过对随机参数的控制得到所需要的介质特点，如非均质性和各向异性。20世纪90年代，随着分形技术的充分发展，以及对裂隙面分形特点认识的逐渐成熟，人们开始采用分形技术对裂隙面的粗糙度和裂隙侧面进行分析和模拟［93^-95］。众多的研究表明隙宽满足对数正态分布或者高斯分布［26^，96］，同时，隙宽分布的方差较大，说明了模拟隙宽非均质性的重要性［56］。