基于分形理论的裂隙岩体地下水溶质运移模拟

第23卷 第8期岩石力学与工程学报 23(8)：1358～13622004年4月 Chinese Journal of Rock Mechanics and Engineering April ，20042003年8月4日收到初稿，2003年10月5日收到修改稿。

* 国家自然科学基金(50179010)、国家自然科学重点基金(50239070)、河海大学科技创新基金(2003410943)和河海大学双聘院士基金资助项目。

作者 王锦国 简介：男，29岁，博士，1996年毕业于河海大学水文地质工程地质专业，现任讲师，主要从事地下水科学与工程方面的研究工作。

E-mail ：wang_jinguo@ 。

基于分形理论的裂隙岩体地下水溶质运移模拟*王锦国 周志芳(河海大学地质及岩土工程系 南京 210098)摘要 岩体裂隙面是粗糙不平的，具有分形特征。

裂隙本身可以看成由上下两裂隙面叠合而成的，可以应用分形几何理论来模拟粗糙裂隙面和裂隙张开度的分布情况。

在用分形理论模拟裂隙面及裂隙张开度的基础上，用特征有限元方法对粗糙裂隙中的溶质运移进行了模拟。

模拟结果显示，考虑裂隙面的粗糙度，模拟其中的溶质运移更符合实际情况。

与平均张开度下光滑裂隙中的溶质运移相比，粗糙裂隙中的浓度锋面更落后，而且存在不均匀性和各向异性的特点。

关键词 岩石力学，分形，粗糙度，示踪剂，各向异性分类号 P 641.2 文献标识码 A 文章编号 1000-6915(2004)08-1358-05SIMULATION ON SOLUTE TRANSPORT IN FRACTURED ROCKSBASED ON FRACTAL THEORYWang Jinguo ，Zhou Zhifang(Department of Geology and Geoengineering ，Hohai University ， Nanjing 210098 China )Abstract Rough fracture surface is of characteristics of fractal and can be considered as congruence of the upper and lower fracture surfaces. The rough fracture surface and the distribution of fracture apertures can be simulated by fractal theory. The simulation is made on solute transport in rough fractures by characteristic FEM based on the fracture face and fracture apertures simulated by using fractal theory. The result of simulation shows that it is more actual to simulate solute transport in fractures with its roughness considered. In contrast with the solute transport in smooth fractures with average aperture ，the frontal surface of concentration in a rough fracture is more behindhand and is of characteristics of nonuniformity and anisotropy. Key words rock mechanics ，fractal ，roughness ，tracer ，anisotropy1 前 言在核废料地质贮存、垃圾填埋、咸卤水入侵等资源、环境问题研究领域中，污染物随着地下水在裂隙介质中的运移是近年来研究的热点问题之一。

由于裂隙介质本身的复杂性，使得这项研究变得非常困难。

单裂隙系统中溶质运移是研究裂隙介质中溶质运移机制的基础，有助于加强人们了解溶质在裂隙中的对流-扩散机理。

从20世纪60年代开始，许多学者就对单裂隙系统中溶质、放射性元素的运移规律作了室内和现场试验研究[1～4]；考虑溶质从裂隙向基质的扩散、裂隙表面的吸附作用、水流速度变化、弥散度变化、放射性衰变等影响因素，也有大量的数值模拟研究成果[5，6]。

然而，以往的模拟均假设示踪剂运移在平直光滑的裂隙中进行。

考虑裂隙表面形态的变化，文[7]假定裂隙张开度服从对数正态分布，模拟了粗第23卷 第8期 王锦国等. 基于分形理论的裂隙岩体地下水溶质运移模拟 • 1359 •糙裂隙中的溶质运移，并与室内试验数据作了对比，取得了较好的效果。

裂隙结构面粗糙度是一种自然分形现象。

近年来，文[8～11]运用分形几何研究结构面粗糙度特征，得到了结构面粗糙度与分维数D 的关系式。

本文应用分形理论对粗糙裂隙面进行了模拟，采用特征有限元法模拟了分形裂隙中示踪剂溶质运移的浓度分布情况。

2 裂隙几何形态的模拟2.1 裂隙面的模拟由于裂隙面是粗糙不平的，具有一定的分形特征，而裂隙本身可以看成由上下两裂隙面叠合而成的，因此，可以根据分形几何理论来模拟裂隙面及裂隙张开度的分布情况。

分形几何法就是根据裂隙面分形维数、裂隙面位错值和裂隙面上的面积接触率，用逐次累加法或随机布朗函数法来获得裂隙张开度分布情况。

用一般的函数难以描述天然裂隙面复杂的形态。

分形理论可用来描述极不规则的几何图形，而且许多研究都表明，用分形几何来模拟粗糙裂隙面是合理的。

目前，自仿射分形模型被认为是模拟裂隙面的最好的分形模型[9]。

下面应用随机布朗函数法(Weierstrass-Mandelbrot 函数)来模拟分形裂隙面。

由分形几何理论可知，如果裂隙面的分形维数为D ，则指数为H 的随机布朗函数可按下式产生具有分形特性的裂隙面：=)(y x z ，∑∞=−++1])sin cos (sin[k k k k kHkkA B y B x C λλ(1)式中：k C 为相互独立的服从标准正态分布的随机数；A k 和B k 为相互独立的服从[0，2π]上均匀分布的随机数；H = 3－D ，D 为裂隙面的分维数。

对于某一D ，给定一个λ(一般5.10.1＜＜λ)就可求出任意一点处的曲面高度z 。

遍布(x ，y )就可获得具分形特征的裂隙面。

图1为按随机布朗函数法模拟生成的裂隙面，其中，裂隙面的分维数D = 2.3，参数25.1=λ。

2.2 裂隙张开度分布的模拟由于裂隙面的不同部位其张开度是变化的，因此，采用位错法模拟裂隙张开度的分布更为合理。

根据布朗曲面的性质，若上、下2裂隙面的位错值为)(y x ΔΔ，，则)()(y x z y y x x z z ，，−Δ+Δ+=Δ(2)图1 分形裂隙面的生成(D = 2.3)Fig.1 Generation of fractal fracture face (D = 2.3)服从以零为均值，以H y x )(22Δ+Δ为方差的标准正态分布。

按上式产生的z Δ中，有正、负2个序列。

设负序列中绝对值最大的负数的绝对值为*z Δ，于是裂隙张开度)(y x b ，可按下式估算：*)()(z r y x z y x b Δ+Δ=，， (3)调整r 的值，使)(y x b ，满足裂隙面上的面积接触率(张开度小于等于0认为接触)的要求，即可生成裂隙张开度分布的情况。

当r =1时，裂隙面2壁之间只有一个接触点，即接触面积为0。

事实上，天然裂隙的面积接触率一般都大于零。

因此，应当有10＜＜r 。

将式(3)等号左端的)(y x b ，作为张开度的统计平均值，则给定x Δ，y Δ，就可优化确定出r 。

取位错值=Δx 0.15 mm ，=Δy 0.15 mm ，面积接触率12.24%，此时r = 0.34，平均张开度0.142 mm ，采用上述方法模拟得到的裂隙张开度分布如图2。

图2 分形几何法产生的裂隙张开度分布图 Fig.2 Fracture aperture distribution generated by fractalmethod• 1360 • 岩石力学与工程学报 2004年3 示踪剂溶质运移模拟3.1 示踪济运移模型将裂隙面离散成许多相等大小的正方形网格，认为每一个正方形网格上的张开度为常数。

由于每个单元很小，其渗透系数可由立方体定律给出，即)(122i i i y x b g k ，μ= (4)式中：)(i i y x b ，为第i 个单元形心处的张开度；g 为重力加速度；μ为流体的运动粘滞系数，当溶质为示踪剂时，在15℃下，/s cm 4 011.02=μ，当溶质浓度较高时，μ为一随浓度变化的变量。

裂隙中的水流运动可以看作是2个隔水板之间的承压水稳定流运动，则水流控制方程可以描述为0=⎟⎟⎠⎞⎜⎜⎝⎛∂∂∂∂+⎟⎠⎞⎜⎝⎛∂∂∂∂y H b k y x H b k x i i i i (5) 地下水溶质运移控制方程为−⎟⎟⎠⎞⎜⎜⎝⎛∂∂+∂∂∂∂+⎟⎟⎠⎞⎜⎜⎝⎛∂∂+∂∂∂∂y C D x C D x y C D x C D x yy yx xy xx tCy C v x C v y x ∂∂=∂∂−∂∂)()( (6) 地下水运动方程为⎪⎪⎭⎪⎪⎬⎫∂∂−=∂∂−=y H k v x H k v iy ix (7) 式(6)，(7)中：C 为溶质浓度；x v ，y v 为地下水流速度；xx D ，yy D ，xy D 和yx D 为水动力弥散系数张量的分量，可表示为⎪⎪⎪⎭⎪⎪⎪⎬⎫++=+−==++=*2L 2T *T L *2T 2L )(1)(1)(1D v v v D D v v v D D D v v vD yx yy y x yx xy y x xx αααααα (8)式中：L α，T α为纵向和横向弥散度；*D 为分子扩 散系数；v 为单元平均水流速度，且有y x v v v =。

3.2 数值模拟及结果分析在裂隙面形态及张开度分析的基础上，采用有限元和特征有限元分别对裂隙中的水流和溶质运移进行模拟。

假设裂隙面的分维数为2.3，2裂隙面的面积接触率为12.24%，将50 cm ×50 cm 的区域剖分成400个四边形单元。

裂隙面形态及开度分布如图1，2所示。

计算中，纵向和横向弥散度分别取5.0L =αcm ，2.0T =αcm为了说明溶质运移的各向异性现象，计算考虑了2种边界条件。