岩体分形裂隙网络系统中水流动研究进展

第43卷增刊1
2011年3月
哈尔滨工业大学学报
JOURNAL OF HARBIN INSTITUTE OF TECHNOLOGY
Vol.43Sup．1
Mar．，2011
岩体分形裂隙网络系统中水流动研究进展
于
贺，李守巨，满林涛，崔
娟
（大连理工大学工业装备结构分析国家重点实验室，辽宁大连116024）
摘要：以分形几何理论和岩体裂隙系统统计地质资料为基础，建立描述岩体裂隙网络系统几何特性的分
形模型．考虑到裂隙长度、
开度、位置和间距等几何参量的分形特性，建立基于Monte －Carlo 模拟的岩体的裂隙网络系统模型随机生成方法．回顾了岩体裂隙网络系统几何参数的分形分布模型以及离散分形裂隙网络系统建模方法、单个裂隙以及裂隙系统水流动数值模拟方法和裂隙系统等效渗透系数估计方法．对岩体裂隙网络系统的分形几何结构进行研究，使得裂隙网络几何模型真实反映岩体裂隙的结构特征、水流动和溶质传输数值模拟模型更接近岩土工程实际原型、预测结果更加准确，在物理本质上揭示分形裂隙岩体水流动和溶质传输机理．
关键词：裂隙网络；分形模型；水流动；溶质传输；数值模拟文章编号：0367－6234（2011）增刊1－0094－06
Some progress of water flow in fractal fractured rock mass
YU He ，LI Shou-ju ，MAN Lin-tao ，CUI Juan
（Dalian University of Technology ，State Key Laboratory of Structural Analysis for Industrial Equipment ，Dalian 116024，China ）
Abstract ：Based on fractal theory and in －situ statistical data of fracture system in rock mass ，the fractal model of the geometry parameters of rock fracture system is proposed．Taking account of the fractal characteris-tics of length ，width ，spacing and location of fracture system ，the fracture network system randomly generated is reconstructed by using Monte-Carlo simulation algorithm．The distribution models of fracture network sys-tem ，reconstruction of fractal geometry system ，numerical simulation methods of water flow in a single fracture and fracture system ，and estimating equivalent permeability of fracture system are discussed ，respectively．Based on randomly walking method ，the randomly walking characteristics of solute in rough surface fracture are analyzed．In order to describe accurately the influence of fracture surface roughness to the solute sorption prop-erties ，the fractal sorption model is proposed by adding surface fractal dimension to the classical sorption mod-el．
Key words ：fracture network ；fractal model ；water flow ；solute transport ；numerical simulation
基金项目：国家重点基础研究发展规划项目（2007CB714006）；国
家自然科学基金重点项目（No．90815023）
1引言
随着经济的快速发展，人类活动正在对裂隙岩体地下水环境造成前所未有的威胁，诸如垃圾填埋山的开采和选矿、煤炭地下气化、铀矿地浸开采、农田过度施化肥和农药的渗滤等．这些与人类生活密切相关的环境问题迫切要求我们对裂隙岩体地下水的流动机理、溶质传输机制进行研究，以
便对其进行预报和有效控制．地下水的污染，不仅
会影响城市生活供水、工业用水以及农业用水，而且还会造成严重的环境问题．地下水是地球上最
重要的水源之一，
也是人类赖以生存的淡水资源之一．在美国，
地下水约占淡水资源的96%，约有97%的城市人口已地下水作为饮用水．在我国北
方，
80%以上的人口以地下水作为饮用水［1］．岩体中含有许多断层、裂隙、节理和层理等不连续
面，在岩石水力学中一般都统称为裂隙［2］
．裂隙岩体中水的流动和溶质传输过程与土体中的渗流过程有着本质的区别．土体渗流学把土体视为多
孔介质，严格地说视为连续介质，因为它并不研究水是如何从一个孔隙流向另一个孔隙的．裂隙岩体系统不同于土体多孔材料，它不是由细小颗粒或散体颗粒组成的多孔介质，而是由随机分布的裂隙和被裂隙切割的岩块组成的不连续介质，裂隙岩体中地下水的流动完全不同于多孔介质中的渗流．其特点在于，裂隙岩体水的流动具有定向性、高度的非均匀性和各向异性性质［3］．等效连续介质模型忽略了裂隙系统和孔隙系统之间的水交换过程，认为裂隙岩体是一种具有连续介质性质的物质，将裂隙中的水流量等效平均到整个岩体中，将岩体模拟为具有等效渗透张量的各向异性连续体，然后运用经典的连续介质渗流理论进行分析．该模型的突出特点是可以沿用各向异性连续介质理论进行分析，无论在理论上还是在求解方法上均有雄厚的基础和经验，而且不用知道每条裂隙的分布特征和水力特性．但是，该模型在应用时有两个问题需要慎重考虑，一是由流量等效确定出的等效渗透张量并不能保证流速和渗透压力的等效性，有时会掩盖裂隙的局部渗流行为；二是等效连续介质模型的适用性和有效性问题．表征单元体（Representative elementary volume，REV）是判定裂隙岩体是否可以作为等效连续介质的标准．如果裂隙岩体系统存在表征单元体，且其尺度远小于研究区域的尺度，则裂隙岩体系统可以作为等效连续介质来处理．与孔隙介质不同的是，裂隙在岩体中所占的体积比例极小，REV 的尺度通常在几十米到上百米，甚至有时不存在，所以作为非连续离散介质处理裂隙岩体系统的水流动问题是非常重要的［4］．离散裂隙网络模型着重考虑裂隙的导水作用，视岩体为简单的裂隙介质，认为其渗流行为由裂隙决定，因而裂隙岩体的渗流特性与裂隙的开度、方位、连通率等参数有关．该模型以单个裂隙内水流基本公式为基础，利用立方定律和广义的Darcy定律来建立流量平衡方程，然后求解方程组获得各裂隙交叉点的水头值．离散裂隙网络模型抓住了裂隙流的本质特征，能够较好的反映水流在裂隙中的运动规律、由于该模型是建立在单裂隙水流运动规律基础上的，力图反映裂隙网络中各点的真实渗流形态，因而具有拟和性好、精度高等优点．其缺点在于，当裂隙数量较多、分布较复杂时，其计算工作量较大．除了等效连续介质模型和离散网络模型之外，还有裂隙－孔隙双重介质模型和分区混合介质模型．由于地下岩体裂隙结构分布的复杂性和随机性有某些局限性和片面性，很难从整体上去把握裂隙岩体水流动和溶质传输特性．分形几何模型通过分数维来描述自然界中广泛存在的随机、无规则而具有近似相似性或者统计意义上自相似性的系统，并在深层次上揭示自然界所遵循的自相似性规律，即“尺度变换不变性”，为裂隙岩体中裂隙结构拓扑特性的研究提供了一种全新的方法，它将不规则的几何体集合视为分形体，按分形几何的方法定量求解分形体空间分布特征参数 分形维，能够有效的刻画出裂隙岩体的几何体的结构形态．
2岩体裂隙网络系统几何特性
岩体裂隙水流动离散裂隙网络模型的研究主要开始于20世纪60年代，最初的离散裂隙网络模型为确定性的，由3个正交方向上的等距平面组成．之后，具有纯随机参数的模型得到了使用．最初的随机模型只有几个随机参数，其余的参数在空间上是不变的，尤其是假定裂隙是无界的，这显然会带来很大的误差．有界裂隙是随着两个模型的产生而引入的：一个是Baecher圆盘模型，另一个是Veneziano多边形模型．大多数随后的用于裂隙岩体渗流、岩体稳定性和岩体变形分析的裂隙网络模型都是基于这两个模型．Veneziano模型在三维空间中的几何形状是相当复杂的，因此大部分研究者使用相对简单的Baecher模型．之后，许多研究者对裂隙网络模型作了发展，并用其来解决实际工程问题．Niemi（2000）指出基于现场尺度的裂隙系统的建模方法包括以下三种：1）确定性的等效多孔连续介质方法，2）裂隙网络方法，和3）随机连续介质方法．在确定性多孔介质方法中，介质的特性常常假设为固定和已知的．由于计算容量的限制，裂隙网络方法常常适合应用到几百米尺度范围．Niemi基于现场三个不同观测尺度的系统实验数据，研究了低渗透率裂隙岩体的水流动特性和尺度相关特点．初步分析结果表明，用标准的连续模型方法解释一维钻孔观测数据可能大大低估了非均质、非连续裂隙系统的三维联通特性［5］．Dreuzy（2004）通过观察裂隙岩体的露头发现裂隙在所有尺度上（从微观到宏观尺度）是相关联的．这种相关模式常常可以用分形维来表征，而裂隙长度分布符合幂指数律．他采用数值模拟方法研究了分形相关模式和幂指数律长度分布对裂隙网络等效渗透率的影响．研究结果表明，增加裂隙的相关性会使得裂隙系统的连通率和渗透率降低．Dreuzy还认为裂隙在空间的
·
59
·
增刊1于贺，等：岩体分形裂隙网络系统中水流动研究进展
分布并不是完全随机的，而是具有相关性，相关性来自于生成裂隙的力学相互作用．裂隙中心的相关性存在于裂隙系统的所有尺度，并且显示出分形特性．更精确的说，裂隙中心的分布是一个用分形维表征的分形概率场．现场观测结果表明，对于平面模型该分形维的值介于1.5至2.0之间［6］．常用的裂隙长度与裂隙数目的分形几何关系如下［7］
N=Cr－D l（1）式中：N为裂隙的数目；C为常数，它表征了岩体的裂隙密度特性；D l为裂隙长度分形维．同样，裂隙开度与裂隙数目的分形几何关系如下
M=βw－D w（2）式中：M为裂隙的数目；β为常数；D w为裂隙开度分形维．
在工程岩体中，岩体裂隙面的分布极为复杂，但从统计意义上来说，反映岩体裂隙面分布特征的几何参数服从一定的概率分布模型［8］．因此，可以通过现场实际调查和室内统计分析，建立裂隙面几何参数的概率统计模型；进而应用随机的方法，在计算机上求得表征裂隙面分布特征的裂隙网络图形．野外露头观察是最普遍和最常用的工程尺度结构面调查方法．工程地质测绘通过调查工程地质条件在地表的分布推测在地下的分布，是结构面调查的基础，也是最主要、最主要的手段．可以从宏观角度研究工程地质条件及其演变规律以及可据此推断地层和构造在地下的分布，建立三维空间地质结构模型．运用数码相机照相进行开挖面工程尺度地质调查是目前广泛采用的工程地质调查方法之一．特别是对于开挖高度大、坡面陡而难以人工直接测量的露头，用数码照相方法就容易、方便地记录开挖面上结构面的分布．工程地质钻探是工程地质勘查的主要手段之一，也是地质体三维建模中获取地质信息的主要手段，通过岩心可以直观地确定地层岩性、地质构造和岩体风化特性，并可以通过钻孔摄影与钻孔电视对岩层的裂隙方向、裂隙发育程度、断层破碎带及软弱夹层等提供直观影像资料．Cowie（1995）采用数值断裂模型模拟了断层的增长过程，认为断层的尺度具有幂指数律特性和分形特征．断层的长度和断距的相关性使得断层系统具有强烈的多重分形模式．提出了断层系统的容量维、信息维和相关维，信息维和相关维与断层系统断裂的单元群和它们的断距有关．裂隙面几何参数模型与统计参数的确定主要包括极大似然估计法和最小．的统计参数使得似然函数获得极大值，具有明确的物理意义．而最小二乘法得出的统计参量可以使得理论分布函数与经验分布函数之间达到最佳的吻合［9］．按照谷德振（1983）的分级方法，可以将结构面分为五级，从I级几千米长度（宽度几十米）的断层构造，到V级延伸5m之内细小结构面，包括隐节理、劈理和片理等．本课题研究的对象是地下工程、大坝工程、石油工程、水利工程和采矿工程中经常遇到的工程尺度裂隙系统问题．裂隙方位的分布服从Fisher分布［10］
θ=cos－1
ln［e k－F（e k－e－k）］
{}
K
（3）
式中：K为Fisher常数，k是分布模型中的参数．裂隙网络系统中每个裂隙中心的位置采用递归方法随机生成．对于二维问题，每个裂隙中心的位置为
x
i
=x
g1
+R
i
（x
g2
－x
g1
）（4）
y
i
=y
g1
+R
i+1
（y
g2
－y
g1
）（5）其中：
x
g2
＞x
g1
，y
g2
＞y
g1
，R
i+1
=27.0R
i
－int（27.0R
i
）Ri是介于［0 1］之间的随机数，int是取整函数；而初始值R0是采用乘积同元算法得到的．基于随机布朗函数，裂隙的表面形状描述为
z（x，y）=∑ɕ
k=1
C
k
λ－αk sin［λk（xcosB k+ysinB k）+A k］
（6）式中：C k为相互独立的服从标准正态分布的随机变量；A k、B k为相互独立的在［0，2］上服从均匀分布的随机变量，λ为某一常数，1.0＜λ＜1.5，α=3－Df；Df为裂隙表面分形维．裂隙间距的累计频率分布服从幂指数函数，也就是满足分形分布模式［7］
N
s
（S）∝r－D s（7）式中：S是裂隙的间距值，D s是裂隙间距大于S 的累计数目，D s是裂隙间距分形维．
裂隙的产状定义了裂隙空间分布状态，用倾向和倾角两个指标表示．裂隙的产状决定了裂隙水流动的方向性和强烈的非均质特性．裂隙密度反映了岩体完整程度和岩石块体大小的重要指标，也可以用裂隙的间距描述．裂隙长度是指岩体中裂隙的延伸范围，通常指沿裂隙走向或者倾向的延伸长度．研究表明，裂隙长度分布具有幂指数律特性和分形特征．裂隙开度是指裂隙结构面相邻岩壁间的垂直距离，同样具有分形特性．对于粗糙的裂隙，裂隙宽度并不是一个常值，而是随着裂隙面上的坐标位置变化的函数，此时可以采用概率分布密度函数描述裂隙开度的变化情况．常用
·
69
·哈尔滨工业大学学报第43卷
概率分布形式包括对数正态分布和负指数分布．同时，裂隙开度的大小与裂隙的长度具有相关性；裂隙的长度越大，裂隙的平均开度越大．光滑平直的裂隙面为2维，极端粗糙的裂隙面接近3维．而实际粗糙裂隙面的维数介于2 3维之间，是一个分数维，不是一个整数维．因此，可以采用分形维来表征裂隙面的粗糙程度．裂隙的组数是指组成交叉裂隙系统的裂隙组数目．原生裂隙面和某些次生的裂隙面分布的规律性较强，属于系统裂隙面；而大多数构造成因裂隙面呈随机分布，属于随机裂隙面．分形几何为表征和模拟自然界中许多复杂非欧拉（non－Euclidean）形状的几何特性提供了合适的数学框架，并且已经被证明非常适合类似于岩体裂隙等自然界中的非连续体［11］．
3离散裂隙网络系统中水流动模拟
离散裂隙网络模型是把裂隙岩体看成由不同规模、不同方向的裂隙个体在空间相互交叉构成的网络状系统，成为裂隙网络，地下水和溶质沿裂隙网络运动．该模型将裂隙及其交叉点上的水动力关系逐个列出，揭示了裂隙水流和溶质运移的内在联系，为裂隙网络流模型的研究奠定理论基础．岩体裂隙网络渗流的特点是渗流沿裂隙网络定向流动，体现出明显的非均质性和各向异性．Davy（1992）的实验研究表明，断层模式具有自相似特性，并且能够用分形维表征．他通过分析大量的实验结果发现与广义分形维、重心分形维和断层长度分布指数相关的尺度律．为了定量研究断层模式，根据实验获得的数据给出了断层长度分布的幂指数表达式．表达式中的幂指数b就是断层重心的分形维，幂指数a表示了断层的数量与尺度之间的相关关系，裂隙累计长度的分形特性用容量（capacity）分形维表征，建立断层系统的广义维与两个幂指数之间的关系，广义维反映了具有分形结构的断层系统的均质或者非均质特性［12］．周志芳认为等效连续介质模型就是把岩块－裂隙系统等效成连续介质，用连续介质理论描述裂隙中水的流动问题．其优点在于可以直接引用比较成熟的孔隙连续介质渗流理论，操作性好；其缺点在于适用范围常常受到限制，因为并不是所有的裂隙岩体均可以等效为连续介质模型，并且典型单元体的大小和等效水力参数较难确定，更为严重的是把裂隙岩体等效成连续体，不能很好地刻画裂隙的特殊导水作用［3］．张有天提出裂隙网络系统虽然十分复杂，但是采用有限元等数
：裂隙网络由一维单元组成；三维裂隙网络由二维单元组成，甚至可以进一步简化为一维单元．对于三维裂隙网络系统，常规的方法是按照各裂隙圆盘的交线对每个圆盘剖分为三角形单元，然后用有限元方法求解．也可以假设裂隙为多边形，生成裂隙网络，对每个裂隙求与其他裂隙相交的交线．在一个裂隙内，用一个管道连接两条交线的中心点．如此进行，三维裂隙网络就简化为管单元组成的空间网络［2］．Jing（2001）采用DDA方法模拟了裂隙岩体中水流动和变形耦合问题．基于经典的平行、光滑裂隙的立方定律（Cubic law），对于非平行（楔形）、光滑裂隙的等效导水系数进行了推导，给出了解析表达式．基于图论理论，对于二维裂隙系统建立了水流动问题线性方程组及其求解方法，并且数值验证了方法的有效性［13］．Bagh-banan（2007）研究了裂隙岩体裂隙开度与迹长之间的相关性，建立了裂隙迹长的分形表达式．基于经典的适用于各向异性和均质的多孔介质的Dar-cy定律，数值分析了裂隙系统等效渗透率的变化规律，研究了对数正态分布裂隙开度二阶矩对表征单元体（REV）存在性的影响以及裂隙网络系统的等效渗透率张量的可能性［10］．Zazoun （2007）采用二维几何分析方法研究了某区域的裂隙网络系统的几何特性．研究结果表明，表征裂隙中心空间分布的模型符合幂指数分布和具有分形特性．同时，裂隙长度和裂隙间距分布同样服从分形几何律［7］．
对于光滑、平直的单个裂隙，等效导水系数为
c=
g
12μ
b3（8）式中：c为裂隙的等效导水系数，g为重力加速度，b为裂隙的开度，μ为水的运动粘性系数．对于光滑非平行的裂隙，其等效导水系数的计算方法为：将上式中的裂隙开度用等效裂隙开度替换即可．其中，等效裂隙开度为
b
e
=b
m
/F（9）
b
m
=（b
a
+b
b
）/2（10）
F=
（1+r）4
16r2
（11）
b
a
，b
b
分别为楔形裂隙的最小开度和最大开度，
r=b
a
/b
b
．
于青春（2006）在建立渗流模型时，首先产生圆盘形三维裂隙网络，然后令圆盘形裂隙网络退化为三维空间上的管状网络．管元的直径通过模型校正的办法确定，即通过优化方法寻找适当的管元直径，使计算的渗流试验排泄流量值能够满
·
79
·
增刊1于贺，等：岩体分形裂隙网络系统中水流动研究进展
意地拟合实测值．由于不可能校正每个管元直径，所以假设在一组裂隙内管元直径服从对数正态分布，这样渗流模型校正参数为每组裂隙管元直径的平均值和标准差［14］．何杨、柴军瑞（2005）利用Monte－Carlo模拟岩体中裂隙网络分布情况，基于裂隙岩体网络非稳定渗流数学模型，编制相应的程序，并通过算例进行验证，进而对工程实例进行分析．模拟计算结论表明，运用Monte－Carlo 方法模拟岩体中裂隙网络分布状况并进行渗流分析是一种较好的分析方法［15］．因此，根据裂隙网络系统的裂隙开度、间距、位置和长度的分形分布模型，采用Monte Carlo随机模拟方法生产裂隙网络，基于质量守恒定律，在每个节点处流入的流量与流出的流量相等，对于每个节点列出平衡方程，得到线性方程组，然后采用Gauss－Seidel迭代方法求解，模拟分形裂隙系统中水的流动过程．Li （2008）研究了分级的分形裂隙网络系统的等效渗透率参数估计方法［16］
k e =N
1
w2
1
12μL
1－（N
s
λ2）m
1－N
s
λ2
（12）
式中：k e为分形裂隙网络系统的等效渗透率，Ns 为裂隙的分支数目，N1为第一级裂隙的数目，m 为裂隙网络系统总的裂隙数目，λ为裂隙开度的比例因子，w1为第一级裂隙的开度，μ为裂隙的弯曲因子；L为裂隙的长度．
用计算机生成的岩体裂隙网络系统，自然不是岩体结构面的真实网络，它只是与岩体具有相同裂隙分布规律的随即抽样．根据Monte－Carlo 原理，大量的重复这种抽样，那么所获得的裂隙岩体水流动参数的平均值就代表了真实岩体的相应的量．有了这样的模拟裂隙网络系统，就可以定量进行裂隙岩体水流动和溶质传输的研究．Yu （2005）基于分形理论、Monte Carlo模拟方法和修正的Hagen－Poiseulle方程以及Darcy定律，建立了多孔介质渗透率的分形解析表达式，Monte Carlo模拟结果与分形解析解所得到的渗透率进行了对比，两者具有良好的一致性，在多孔材料的热传导和渗流等方面郁伯铭教授取得了许多创新性的研究成果［17］．Ghodoossi（2005）研究了微孔道分形网络中热和气体传输特性，提出了多级分形网络的产生方法，建立了热传导参数和气体流动阻力的分形描述机理［18］．Peratta（2006）提出了一种模拟三维裂隙中流体流动的方法，模型中采用了离散元和边界元方法，以便精确和有效的处理介质三维复杂的几何形状和非均质特性［19］．（）测量方法，分析了其宏观和微观特性，对裂隙的质量分形维、密度分形维、长度和方位分形维进行了讨论和分析，以及这些分形维对裂隙岩体等效渗透率的影响［20］．冯学敏（2006）利用蒙特卡罗方法随机生成符合给定概率分布特征的复杂裂隙网络；然后通过交切和拓扑运算将各裂隙段置于常规有限单元内部，形成内含由多个裂隙段分划而成的子单元的复合单元，根据推导的公式计算渗流场，进而分析岩体的渗透特性．该方法考虑了岩体中每条裂隙的结构和渗透性质，以及岩石和不连通裂隙对整个区域渗流场的影响，提出用子单元映射节点水头的概念来描述裂隙引起的渗流场的不连续性，求解过程与传统有限元法类似，但相对于离散裂隙网络水力学方法前处理较为简化，具有较好的工程应用前景［21］．李守巨（2007）基于分形理论和达西定律，研究了孔隙的弯曲分形维数和表面积分形维数与瓦斯渗透率参数之间的关系；数值模拟了煤层孔隙的弯曲分形维数和表面积分形维数对瓦斯渗透率参数的影响．研究表明，瓦斯渗透率随着孔隙弯曲分形维的增加而降低，随着孔隙表面分形维和最大毛细管直径的增加而增加；应用分形理论可以获得煤层孔隙发育程度和分布特征比较精确的定量信息［22］．宋晓晨（2004）发展了用于模拟裂隙岩体渗流的三维离散裂隙网络数值模型，并编制了裂隙岩体渗流模拟程序FracFlow．该模型可以利用野外露头上采集到的裂隙的观测数据，通过计算机处理最终形成三维裂隙网络的人工几何模型，然后用边界元法求解所生成的裂隙网络中的渗流问题［23］．赵阳升等（2005）经过十余年的研究，先后对20多个矿井的煤层与岩层进行现场观测和显微观测，发现不同尺度下裂隙面数量服从分形分布规律．以岩体裂隙面数量分布具有分形特征为基础，建立了岩体裂隙面数量的三维分形分布仿真模型，并采用VisualC++6.0开发出岩体裂隙面数量的三维分形分布仿真系统．利用该仿真系统可方便地研究岩体内部任意剖面的裂隙迹线分布特征以及任何子块内裂隙面的分布特征［24］．谢和平和周宏伟（1997）系统地评述了分形理论在岩土渗流方面的应用，介绍了数值计算方法、测试技术和研究进展，并且指出孔隙与破裂介质中流体流动及其输运性质的研究是一个极富挑战性的研究课题［25］．
4结语
引入分形模型到岩体裂隙网络系统中，裂隙岩体水流动和溶质传输过程是在分数维分布的裂
·
89
·哈尔滨工业大学学报第43卷。