图 1 裂隙网络管单元
Fig.1 Pipe element of fracture network
表 1 裂隙网络参数统计
Table 1 Statistical parameters of fracture network
WANG Junqi,LI Pengyu,DONG Ye,et al.Study on the seepage tensor of multi-parameter model of pipe elements in fractured rock mass[J].Journal of Zhejiang University of Science and Technology,2017,29(05):358-365.[doi: 10.3969/j.issn.1671-8798.2017.05.007 ]
[PDF全文]
根据现场实测的岩体裂隙产生状况,采用圆盘裂隙假设,用Monte-Carlo方法生成3组裂隙形成裂隙网络。基于水在裂隙中的沟槽流形态将三维裂隙网络简化成空间一维管单元模型,且认为管单元的等效管径与隙宽呈线性关系。设置每组裂隙等效管径与隙宽比值为待求参数,利用遗传算法反演得到,进而求取等效渗透张量,并与整体单参数模型获得的等效渗透张量相比较。与工程实测渗透张量对比,结果表明多参数模型较单参数模型能够更好地表征渗透张量。
Based on the fracture occurrence of rock mass measured on the spot, the Monte-Carlo method is employed to generate three groups of fracture networks by referring to the discfracture hypothesis. In accordance with the hypothesis of trench flow in the fracture, the 3-D fracture network is simplified to a model consisting of one-dimensional pipe elements, and the equivalent diameter of the pipe elements is considered to be linear with the fracture width. The ratio of equivalent diameters to apertures of each fracture being set as the unknown parameter, the equivalent seepage tensor is obtained by inverting the genetic algorithm,which can be compared with the equivalent seepage tensor obtained by the whole single parameter model. In comparison with the field data of seepage tensor, the results show that the multi-parameter model can better characterize the seepage tensor than the single parameter model.
渗流是流体通过多孔介质或裂隙介质的流动,裂隙岩体中流体的渗流分析对人类生产生活非常重要,如石油储存、核废料地下贮存等方面[1],因此裂隙岩体的渗流问题一直受到各国学者的广泛关注。研究发现岩体中结构面具有不连续性且发育极不规则,裂隙网络的几何构成非常复杂。为方便研究裂隙的渗透性,通过计算机技术生成离散裂隙网络模型,对仿真模型渗透特性展开研究的较多[2-4]。在这些研究中,利用渗透张量来表征裂隙渗透特性的方法得到广泛的应用,通过张量的性质可以将非连续裂隙岩体的渗流问题转化为连续介质模型求解渗流场的问题[5],从而简化了计算,使得大规模工程问题的解决成为可能。目前,研究渗透张量主要有实验法和数值法两种方法。在实验法确定渗透张量方面,周志芳等[6]基于振荡式微水试验原理,通过单孔微水振荡试验现场确定了岩体的渗透系数张量。杨金保等[7]对取自丹江口水库库区的辉绿岩剪切裂隙进行了不同围压和裂隙水压下的渗透试验,基于该试验对渗透系数进行考虑三向应力和裂隙水压力变化的修正,得到新的渗透张量的表达式。李克克等[8]通过现场实测获得裂隙结构面参数,采用离散元法进行渗流分析。在数值法确定渗透张量方面,王培涛等[9]应用离散裂隙网络模型,通过平面渗流分析方法研究了节理岩体不同几何分布情况下的渗透张量特征。Baghbanan等[10]基于UDEC计算平台,采用离散裂隙网络模型对渗透张量进行了研究,并对渗透的各异性进行了分析。从文献资料中可以看出,数值法主要是在二维模型的情况下对裂隙渗流进行分析计算,这与实际情况相差较远。在三维模型的研究方面,李新强等[11]使用边界元法来计算三维裂隙网络的渗透张量,通过渗透张量来定量分析裂隙渗流。Lee等[12]根据裂隙的几何分布性质生成裂隙网络,在裂隙圆盘上划分三角形单元网络,用数值方法来模拟裂隙中的复杂流动和二氧化碳的运移问题。Cacas等[13]认为水在裂隙内成沟槽流形态,提出将三维裂隙网络简化成一维管单元,为渗流的分析计算提供了新的思路。Wang[14]等利用压水实验对管单元模型结果进行校正,通过回归分析得到渗透张量。于青春等[15]建立了非连续裂隙网络管单元渗流模型并对模型进行了校正。王俊奇等[16]采用单参数模型,通过反演法确定单元管径,利用有限单元法计算得到裂隙岩体的渗透张量。
由于管单元模型接近水流在裂隙间的流动形态,而且能有效地减少计算量,因此本研究也采用管单元模型进行渗流分析; 不同之处在于,以往的管单元模型中只有一个参数作为等效管径和隙宽的比值,本研究将在三维裂隙网络模型中生成多组裂隙网络,并且对各组裂隙进行标记,设置相应的参数作为管径和隙宽的比值,通过有限单元法研究多参数模型确定的渗透张量。
1 三维裂隙网络模型的建立与渗流计算方法在工程实践中,岩体结构面的分布表面上看起来是杂乱无序的,但是通过将其分布的特征数据进行汇总可以发现,结构面的分布符合一定的概率特征。通过野外的实地测量采样,进行统计分析,建立岩体几何参数的统计模型,利用计算机技术生成统计特性与实地测量相符的裂隙网络模型[17]。
1.1 裂隙网络模型的生成经过多年的发展,目前有多种岩体结构面的几何系统模型用来表征岩体结构系统特性的特定组合。在众多的模型中,Baecher模型因其简便实用而被广泛应用于岩石力学的相关分析中[18]。本研究采用这种模型假设岩体中的裂隙形状为圆形,通过实地调查确定裂隙几何参数的分布形式,利用Monte-Carlo方法生成裂隙网络,并通过优化方法剔除对导水无用的裂隙。在生成的过程中分别对每组裂隙中的每条裂隙进行编号。
1.2 一维管单元模型平面渗流假定认为裂隙为有厚度的圆盘,水在整个裂隙圆盘上流动。进行渗流计算时需对裂隙圆盘进行网格划分,裂隙数目较多时,有限元计算工作量非常地大,给研究工作带来很多不便。本研究根据前人的沟槽流假设,将水在裂隙中的渗流路径简化为管单元模型,即把相交裂隙的圆心与交线的中心作为管单元节点,两个相邻节点组成一个管单元,将裂隙之间水的流动简化为管元通道之间的流动。根据管元所在裂隙平面的组号再对管元进行分类,如图1所示。
图1 裂隙网络管单元
Fig.1 Pipe element of fracture network
1.3 渗流的计算法
假设管单元的水流通道为圆形,那么通过水流通道的流量[18]为:
q=(ΔH)/L(πD4)/(128)(ρg)/μ。(1)
式(1)中:ΔH为两个节点之间的水头差; L为两节点之间的距离; D为管单元管径; ρ为水的密度; g为重力加速度; μ为水动力黏滞系数(20 ℃时,μ=0.01 g/(cm·s))。
根据式(1)可以得到任意的管单元m中流过两个节点i,j的流量为:
{qi=(Hi-Hj)/(Lm)(πD4m)/(128)(ρg)/μ,
qj=(Hj-Hi)/(Lm)(πD4m)/(128)(ρg)/μ,(2)
式(2)中,节点流量q以流入为正,流出为负。用矩阵可表示为:
[1/(Lm)(πD4m)/(128)(ρg)/μ -1/(Lm)(πD4m)/(128)(ρg)/μ
-1/(Lm)(πD4m)/(128)(ρg)/μ 1/(Lm)(πD4m)/(128)(ρg)/μ]{Hi
Hj}={qi
qj}。(3)
由式(3)集成为总体渗透方程:
[K]{H}={Q},
其中,[K]称为总渗透矩阵,它由所有的单元渗透矩阵[Km]叠加得到。
2 管径反演与渗透张量确定2.1 参数的反演由式(1)可知,进行渗流计算须知道每个管单元直径D的大小。根据等效管径与裂隙隙宽的正比关系,设第i组裂隙的等效管径与隙宽的比值为Ni,i=1,2,3,…,n; n为研究域中裂隙组编号。
基于改进后的遗传算法,以渗透张量的模拟值与实测值之间的误差最小为目标函数,通过反演的方法整体上得到每组裂隙的等效管径与隙宽的比值Ni。由于Ni的改变可以影响渗透张量主值的大小,所以可以将目标函数表示为:
minf=(|k1-k'1|)/(min(k1,k'1))+(|k2-k'2|)/(min(k2,k'2))+(|k3-k'3|)/(min(k3,k'3))。(4)
式(4)中:k1、k2、k3为模型数值计算得到的渗透主值; k'1、k'2、k'3为岩体实测的渗透主值。
2.2 渗透张量的确定节理岩体渗透性具有明显的各向异性,即各个方向上的渗透系数变化明显,目前大多采用渗透张量来表示。依据文献[19]中提到的在单位体积立方体的空间中,在两个相对面边界上施加水头边界,对其进行数值仿真试验,可得到各方向的流量qx、qy和qz,但仅靠这3个数值并不能有效确定裂隙的渗透张量,所以在正六面体试件的某一对相对面(边界)施加等水头边界,其余各面根据其相邻的情况设定相应的变水头边界条件,分别改变3组相对面的边界条件,如图2所示的计算模型进行试验。
图2 计算模型的边界条件
Fig.2 Boundary conditions of computational model
模型边界沿ii方向施加等水头边界,另两个方向为变水头边界,可得到6个流量数据; 进行3组仿真可得到18个数据,然后基于迭加原理确定二阶渗透张量。根据相应的达西定律:
Q=kAΔP。(5)
式(5)中:Q为流量向量; k为渗透张量; ΔP为水力坡降; A为过流面积。
由相应的流量分量,通过矩阵变换得到渗透张量矩阵:
[kxx kxy kxz
kyx kyy kyz
kzx kzy kzz]=[qxx qxy qxz
qyx qyy qyz
qzx qzy qzz][1/(ΔPx)0 0
0 1/(ΔPy)0
0 0 1/(ΔPz)]。(6)
3 算 例3.1 工程数据参照文献[11]通过实测得到裂隙几何参数的信息如表1所示,共测得3组裂隙。
表1 裂隙网络参数统计
Table 1 Statistical parameters of fracture network
文献[11]利用边界元法算得表1裂隙网络的渗透张量,并通过现场压水试验进行校核得到渗透主值及主方向如式(7),作为后面管径反演的目标渗透张量。
{k1=3.08×10-7,倾向170.00°,倾角19.8°;
k2=4.33×10-7,倾向45.9°,倾角57.3°;
k3=4.70×10-7,倾向269.6°,倾角24.8°。(7)
3.2 多参数模型计算采用自编裂隙网络有限元程序,将表1中信息输入到程序中,对裂隙网络进行模拟。程序生成100 m×100 m×100 m的立方体区域,在其中截取20 m×20 m×20 m的立方体区域作为研究域。利用Monte-Carlo方法生成三维裂隙网络。在生成域中共生成结构面200 255条,其中第1组结构面129 006条,第2组结构面38 224条,第3组结构面33 025条。在研究域中,剔除对渗流无作用的裂隙后共生成结构面911条,其中第1组结构面为480条,第2组结构面为250条,第3组结构面为181条,将其简化成管单元模型,则得到4 104个单元,节点数为3 033个,生成的裂隙网络如图3所示。
图3 三维裂隙网络生成
Fig.3 Generation of 3-D fracture network
在数值计算过程中,设定水头边界H1=30 m,H2=10 m。由于模型中生成了3组裂隙网络,根据2节内容所述,设置参数N1、N2和N3分别作为3组裂隙的等效管径与隙宽的比值。以式(7)作为拟合的标准,通过优化拟合,得到3个参数:13.17、12.06及3.01。
通过计算得到多参数模型水头分布如图4所示,渗透张量(m/s)如下:
K=[5.35×10-7 8.16×10-7 1.62×10-7
8.16×10-7 3.24×10-7 1.02×10-7
1.62×10-7 1.02×10-7 3.48×10-7]。(8)
图4 三组模拟水头分布图
Fig.4 Three groups of simulated head distribution
通过渗透张量得到的渗透主值和渗透主方向如下:
{k1=3.08×10-7,倾向142.35°,倾角24.62°;
k2=2.16×10-7,倾向20.36°,倾角49.14°;
k3=6.71×10-7,倾向247.88°,倾角30.30°。(9)
与式(7)数据相比,误差为:
{ε1=0.000 5,θ1=25.99°;
ε2=0.09,θ2=17.18°;
ε3=0.42,θ3=19.87°。(10)
式(10)中:ε为主值相对误差; θ为主方向绝对误差。
3.3 单参数模型计算采用文献[16]的单参数模型方法,对三维裂隙网络渗流模型进行数值模拟。在计算的过程中设置一个参数N作为3组裂隙等效管径与隙宽的比值。在表1所示信息的情况下,得到参数N的值为12.34,其渗透张量结果如下:
K=[5.44×10-7 7.99×10-7 1.12×10-7
7.99×10-7 4.12×10-7 1.08×10-7
1.12×10-7 1.08×10-7 2.31×10-7]。(11)
渗透主系数与渗透主方向如下:
{k1=3.84×10-7,倾向146.7°,倾角10.6°;
k2=1.66×10-7,倾向30.6°, 倾角66.8°;
k3=6.36×10-7,倾向240.7°,倾角20.4°。(12)
与式(7)数据相比,其误差为:
{ε1=0.25,θ1=24.3°;
ε2=1.61,θ2=11.8°;
ε3=0.35,θ3=27.0°。(13)
式(13)中:ε为主值相对误差; θ为主方向绝对误差。
与多参数模型相比,单参数模型渗透主值的相对误差较大。
3.4 隙宽为定值时不同模型的计算3.2和3.3两节中隙宽的分布服从均值为0.000 1,标准差为0.000 1的对数正态分布,本节将裂隙的所有隙宽设置为定值0.000 1,其他几何参数保持不变。研究域尺寸为20 m×20 m×20 m的情况下多参数模型的参数值为19.51、16.55和14.05,渗透张量为:
K=[6.34×10-7 8.37×10-8 4.39×10-8
8.37×10-8 3.36×10-7 3.54×10-8
4.39×10-8 3.54×10-8 3.97×10-7]。(14)
{k1=3.08×10-7,倾向169.19°,倾角22.78°;
k2=3.94×10-7,倾向47.68°,倾角62.31°;
k3=6.66×10-7,倾向252.78°,倾角14.90°。(15)
{ε1=0.000 7,θ1=35.59°;
ε2=0.09,θ2=36.01°;
ε3=0.42,θ3=19.87°。(16)
单参数模型的参数值为16.46,渗透张量为:
K=[4.52×10-6 6.66×10-8 6.20×10-8
6.66×10-8 3.38×10-7 3.04×10-8
6.20×10-8 3.04×10-8 2.67×10-7]。(17)
{k1=3.07×10-7,倾向154.68°,倾角1.39°;
k2=2.47×10-7,倾向59.94°,倾角73.61°;
k3=5.03×10-7,倾向245.09°,倾角16.33°。(18)
{ε1=0.000 9,θ1=23.73°;
ε2=0.042 98,θ2=17.21°;
ε3=0.070 0,θ3=24.4°。(19)
改变研究域的尺寸大小得到多参数模型的参数均值与单参数模型的参数,如表2所示。
表2 模型参数值
Table 2 Parameter values of the model
从式(16)和式(19)的计算结果可以看出,两种模型在定隙宽的情况下,渗透主值和渗透主方向的误差在工程允许范围内。
3.5 分 析通过3.2和3.3节的计算结果可以看出,当隙宽服从对数正态分布时,多参数模型主方向的绝对误差总和为63.04,单参数模型的误差和为63.1,两者相当; 在渗透主值误差方面,多参数模型的误差明显小于单参数模型。在实际情况中,不同的裂隙组之间的管宽比存在差异,多参数模型更加与实际情况相符合,所以在计算结果上比单参数模型更为精确。从两种模型的计算结果可以发现,主方向上的差异比较明显,这是因为在优化的过程中,在使目标函数式(4)达到要求的条件下,对应主方向与目标张量的主方向保持基本一致,即两者夹角为锐角即可。
当裂隙网络的隙宽为定值时,单参数模型的计算结果优于多参数模型,这是因为当隙宽为恒定值时,每个单元的管径与隙宽的比值接近于同一个数,所以采用单参数模型能得到更精确的结果。同时,由表2还发现恒定隙宽情况下,多参数模型的参数均值与单参数模型参数值近似相等,且在研究域尺寸不同的情况下参数值在17.00左右。
4 结 论本研究采用三维离散裂隙网络模型,将复杂的裂隙网络简化为管单元模型,在计算方面设置多个参数得到渗透张量,从算例的计算结果可以看出:当隙宽服从对数正态分布时多参数模型误差较小,这种情况下适宜选择多参数模型; 当隙宽均值为恒定大小时,适宜采用单参数模型进行计算分析。同时,从式(10)的计算结果可以发现,多参数模型的计算结果与工程实测结果之间的误差比较小,有着较高的计算精度。这表明在工程应用上,使用多参数模型进行模拟可以提高渗透分析的准确性,因而具有良好的理论与实用价值。
- [1] JAVADI M, SHARIFZADEH M, SHAHRIAR K. A new geometrical model for non-linear fluid flow through rough fractures[J].Journal of Hydrology,2010,389(1):18.
- [2] 王者超,张振杰,李术才,等.基于离散裂隙网络法的地下石油洞库洞室间水封性评价[J].山东大学学报(工学版),2016,46(2):94.
- [3] 宋晓晨,徐卫亚.裂隙岩体渗流概念模型研究[J].岩土力学,2004,25(2):226.
- [4] WITHERSPOON P, LONG J. Some recent development in understanding the hydrology of fractured rocks[C]//The 28th US Symposium on Rock Mechanics.Tucson:CRC Press,1987:422.
- [5] ANSARI M A, DEODHAR A, KUMAR U S, et al. Isotope hydrogeological factors control transport of radon-222 in hard rock fractured aquifer of bangalore, karnataka[M].Berlin:Springer International Publishing,2016:259.
- [6] 周志芳,庄超,戴云峰,等.单孔振荡式微水试验确定裂隙岩体各向异性渗透参数[J].岩石力学与工程学报,2015,34(2):271.
- [7] 杨金保,李蕾,谢妮.考虑应力状态的裂隙岩体渗透系数张量计算[J].人民黄河,2013,35(3):133.
- [8] 李克克,张杨松.基于实测结构面的岩体渗流特性分析[J].勘察科学技术,2016(1):1.
- [9] 王培涛,杨天鸿,于庆磊,等.基于离散裂隙网络模型的节理岩体渗透张量及特性分析[J].岩土力学,2013,34(增刊2):448.
- [10] BAGHBANAN A, JING L. Hydraulic properties of fractured rock masses with correlated fracture length and aperture[J].International Journal of Rock Mechanics and Mining Sciences,2007,44(5):704.
- [11] 李新强,陈祖煜.三维裂隙网络渗流计算的边界元法及程序[J].中国水利水电科学研究院学报,2006,49(2):86.
- [12] LEE I H, NI C F. Fracture-based modeling of complex flow and CO2 migration in three dimensional fractured rocks[J].Computer and Geosciences,2015,81:64.
- [13] CACAS M C, LEDOUX E, MARSILY G D, et al. Modeling fracture flow with a stochastic discrete fracture network: calibration and validation: 1. the flow model[J].Water Resources Research,1990,26(3):479.
- [14] WANG M, KULATILAKE P H S W, UM J, et al. Estimation of REV size and three-dimensional hydraulic conductivity tensor for a fractured rock mass through a single well packer test and discrete fracture fluid flow modeling[J].International Journal of Rock Mechanics and Mining Sciences,2002,39(7):887.
- [15] 于青春,武雄,大西有三.非连续裂隙网络管状渗流模型及其校正[J].岩石力学与工程学报,2006,25(7):1469.
- [16] 王俊奇,岳潇,李闯.三维裂隙网络管单元模型确定岩体渗透张量[J].华北电力大学学报(自然科学版),2015,42(1):108.
- [17] 张发明,汪小刚,贾志欣,等.三维结构面连通率的随机模拟计算[J].岩石力学与工程学报,2004,23(9):1487.
- [18] 汪小刚,贾志欣,张发明,等.岩体结构面网格模拟原理及其工程应用[M].北京:水利水电出版社,2010:83.
- [19] ZHANG X, SANDERSON D J. Numerical modeling and analysis of fluid flow and deformation of fractured rock masses[M].Amsterdam: Elsevier Science Ltd,2002:53.