锦屏一级左岸岩体渗透特性的数值模拟研究

田洪圆，张 亮，李 宏

(1.大连理工大学 海岸和近海工程国家重点实验室， 辽宁 大连 116024；2.中国电建集团 昆明勘测设计研究院有限公司， 云南 昆明 650041)

锦屏一级左岸岩体渗透特性的数值模拟研究

田洪圆1，张 亮2，李 宏1

(1.大连理工大学 海岸和近海工程国家重点实验室， 辽宁 大连 116024；2.中国电建集团 昆明勘测设计研究院有限公司， 云南 昆明 650041)

基于RFPA2D-Flow软件，并结合现场压水试验，从表征单元体积(REV)及渗透张量入手，对裂隙岩体的渗透特性进行数值模拟研究。以锦屏一级水电站左岸岩体为例，基于岩体裂隙几何特征，运用蒙特卡罗方法生成岩体裂隙网络模型；模拟计算不同尺寸不同方向上的等效渗透系数，依据类常量与类张量特性初步确定裂隙岩体渗透表征单元体积(REV)及渗透张量，再结合现场压水试验结果修正渗透张量。随着模型尺寸的增大，岩体的渗透性逐渐趋于稳定，表征单元体尺寸的大小约为迹长的4倍，表明基于RFPA2D-Flow数值模拟并结合压水试验确定裂隙岩体渗透张量和REV的方法合理有效。

裂隙网络模型；表征单元体积；渗透张量；蒙特卡罗方法；RFPA2D-Flow

裂隙岩体中不同产状、迹长、隙宽等几何特征，使岩体渗透特性在空间上存在客观的各向异性。大型岩体工程的渗流分析中等效连续介质模型应用最为广泛[1]，虽然通过等效连续介质模型并不能对裂隙局部渗流进行详细描述，但它能够阐述清楚裂隙岩体宏观渗流特点，且计算量相对较小。

Sagar Budhi等[2]考虑了裂隙几何特征以及裂隙几何参数引起的渗透系数的误差，提出二阶统计分析方法；Oda Masanobu[3-4]基于裂隙空间分布的随机性，通过统计理论对渗透张量进行推导，同时针对荷载条件作用下岩体渗透性进行了研究；Zhang X等[5]采用多孔介质模型和离散裂隙网格模型，考虑裂隙几何特征和应力，数值模拟得到二维裂隙岩体渗透张量；Chen S H等[6]的复合单元法结合蒙特卡罗方法生成裂隙网络，计算不同尺寸岩体试件不同方向上的等效渗透系数过程中考虑了岩块与裂隙间的相互渗透作用；Rong Guan等[7]基于节理网络模拟技术，应用渗流能量叠加原理，推导得到计算节理岩体渗透张量公式，提出明确裂隙岩体渗透表征单元体积之法，并讨论了裂隙几何特征对岩体渗透特性的影响；王培涛等[8]通过离散裂隙网络模型，构建平面渗流分析法，对节里岩体中的渗透张量进行了探讨，利用对渗流定向性系数进行定义，分析岩体渗流定向特征；杨金保等[9]基于裂隙几何特征的基础上，还考虑了应力作用对渗透系数张量的影响；车富强等[10]把裂隙样本法和现场测试法相结合，用裂隙样本法初步确定裂隙岩体渗透张量，采用现场测试法修正裂隙岩体渗透张量。王俊奇等[11]利用离散裂隙网络的管单元模型研究了裂隙岩体渗透张量，通过渗透椭圆拟合确定。

本文以锦屏一级水电站左岸岩体为例，根据岩体地质勘探资料，生成随机裂隙网络模型，将模型数据导入RFPA2D-Flow软件；选取岩体不同尺寸、不同方向的分析域进行渗流数值模拟，得到等效渗透系数；基于类常量与类张量特性[12]来确定裂隙岩体渗流表征单元体积和渗透张量；再结合现场原位压水试验数据，对模拟结果进行修正。

1 裂隙岩体渗透张量及表征单元体积分析的基本原理

确定裂隙岩体渗透REV、渗透张量应用的方法是：按照现场实测区域中内节理几何参数，通过Monte-Carlo法构建出裂隙网络生成域，然后于该域里的中心位置取出一系列尺度不一的正方形作为REV试算域。

一般情况下，REV的尺寸为裂隙平均迹长的3～4倍[13]，而分析域的尺寸必须大于REV尺寸，才能使计算有意义。所以根据裂隙迹长的统计数据大致确定出最大分析域的尺寸，按大于分析域的2倍大小来确定裂隙网络模型生成域的大小，然后于该域里的中心位置取出一系列尺度不一、不同方向的正方形作为REV试算域。

1.1 裂隙岩体的平面渗流模型

二维裂隙网络模型里，为了便于操作，设定节理是短直线相互交错所形成。根据现场试验资料来看，节理倾角以及间距还有迹长和密度等几何特征做了研究分析，获得了某种概率密度分布，然后再通过Monte-Carlo法生成节理网络，可参考图1。

图1 REV分析模型

根据图1，定水头边界是MM′、NN′，而MM′、NN′间渗透系数如式(1)，不透水边界是MN、M′N′。

(1)

在这个公式里：Δq表示研究区域中产生的总流量，m3/s；ΔH为M′M、NN′二者间水头差，m；K则是MN方向等效渗透系数，m/s；MN、M′M是模型边长，m，二者相等长。

1.2 渗透张量的计算

渗透张量与其他的二阶张量一样，也具有坐标转换特性。图2中假定整体坐标系xoy以原点o为圆心，逆时针旋转角度形成局部坐标系x′oy′，x′轴与x,y轴的方向余弦记为lx′x=cosθ，lx′y=sinθ，y′与x,y轴的方向余弦记为ly′x=-sinθ，ly′y=cosθ。

则[K′]各元素Ki,j与[K]中Kij间满足转换关系：

Ki′j′=Kij·li′i·lj′j

(2)

(3)

整体坐标系下x轴和y轴分别是主渗透系数的方向，则局部坐标系下轴方向的渗透系数为：

Kθ=K1·cos2θ+K2·sin2θ

(4)

极坐标系(l,θ)下，椭圆方程可以表示为：

(5)

图2 坐标系相对位置

图3 椭圆关系曲线

1.3 渗透张量的修正

上述1.1节计算得到的渗透系数并不能代表某一工程岩体的平均渗透性，此处取几何平均值作为岩体等效渗透系数：

(6)

式中：Kω为等效渗透系数；Ki为计算得到的渗透主值。取m=K0/Kω，对渗透张量进行修正，其中K0为现场压水试验得到的渗透系数。

1.4 表征单元体积的确定

裂隙岩体渗透REV即可以在宏观方面等效表现出岩体渗透性相对应的体积最小值，如果岩体尺寸比REV更大时，表明有二阶对称正定渗透张量。

把裂隙岩体等效为连续介质，在此引入类张量和类常量特性[11]。当所研究的岩体的渗透性参数在给定的误差范围内同时满足类张量和类常量特性，此时研究区域大小即为相应的渗透REV。

在裂隙网格生成域内以模型中心为基点，以2 m起始至28 m构成14个研究域，每个对应尺寸每隔30°旋转一次，构成12个方向，通过统计下侧边界流量求得裂隙岩体等效平均渗透系数。此处利用裂隙网络对称性，仅对0°～150°上6个方向进行模拟计算，根据对称性得到其余6个方向上等效渗透系数，确定渗透特性表征单元体REV大小。

2 工程实例

国内现在已建成的拱坝中，最高的即锦屏一级水电站，该水电站的水库正常蓄水位1 880 m，库容77.6亿m3。岩层的倾向在285°～295°范围内，与之相应的倾角则在35°～45°范围内。拱坝左边是反向坡，右边是顺向坡，整个区域中岩体结构面明显，同时受地质构造作用的影响强大[14]，坝址区节理裂隙划分成5组：

(1) N15°—80°E，NW∠25°～45°；

(2) N50°—70°E，SE∠50°～80°；

(3) 近SN—N30°E，SE∠60°～80°；

(4) N30°—50°W，NE∠60°～80°；

(5) N60°W—EW，NE(SW)或S(N)∠60°～80°。

各级岩体岩体结构、风化卸荷以及所处岸坡不同部位情况，各级岩体的渗透系数如下：

Ⅱ级岩体：q<3 Lu，新鲜岩体；

Ⅲ1级岩体：3≤q<10 Lu，深部微卸荷岩体；

Ⅲ2级岩体：10≤q<100 Lu，浅表弱卸荷岩体；

Ⅳ1级岩体：q≥100 Lu，岸坡浅表部强卸荷岩体，卸荷裂隙发育且多张开。

Ⅳ2级岩体：q≥100 Lu，左岸深部裂缝发育段岩体，多松弛破碎张开呈空缝。

根据表1中裂隙的分布规律，可以发现迹长均值为2.29 m，结合参考文献[15]，可取分析域尺寸为10 m×10 m，则生成域应大于分析域的2倍，此处取为25 m×25 m，由于几何参数中缺乏对裂隙张开度的统计数据，所以此处参考文献[2]，取裂隙张开度的估计值微0.2 mm。基于蒙特卡罗方法生成裂隙网络模型，得到裂隙的数据文件，然后将裂隙几何信息导入RFPA2D-Flow软件，边界条件如图1所示，通过数值模拟结算得到裂隙岩体的渗透特性变化如图4所示。

表1 岩体结构面几何参数

图4 锦屏一级左岸岩体不同尺寸渗透特性变化

给定误差限15%，根据图4，由类常量特性图4(a)可见，当尺寸达到11 m×11 m以上，Kxx、Kyy、Kxy满足类常量特性；由类张量特性图4(b)可见，当尺寸达到10 m×10 m，变异系数在允许误差范围内；再通过图4(c)可以进一步验证，当尺寸达到12 m×12 m以后，渗透主值以及两比值系数亦趋于稳定，所以取REV尺寸为12 m×12 m，相当于最大裂隙迹长的4倍，与文献[15]的结论相符合。

上文裂隙张开度取为估计值来确定REV大小无可厚非，但以此得到的渗透张量并不能用来表述整个岩体的渗透特性，因为岩体各岩层、各分区的渗透性亦不同，为了更好的反应岩体渗透特性，结合现场试验数据，将各岩层渗透系数相近者归为一类，取其平均渗透系数近似推算出裂隙开度均值，取分析域尺寸为15 m×15 m，通过数值模拟分别计算各研究域渗透特性，得到初始渗透张量如表2所示。

表2 初步确定渗透张量

由于随机裂隙网络具有随机性，并不能完全反应出真实的裂隙分布以及裂隙连通情况，所以再次根据现场压水试验数据(见表3)进行修正，修正后的渗透张量如表4所示。

表3 各岩层水压试验水力传导系数[16]单位：m/s

表4 修正后各分区的渗透张量

3 结 论

(1) 结合数值试验法和现场单孔压水试验，本文以锦屏一级水电站坝区左岸岩体为例，对渗透张量及渗流REV进行了研究。

(2) 在研究域内，根据岩体内节理的分布特征，运用蒙特卡罗方法生成随机裂隙网络，再导入RFPA2D-Flow进行渗流计算，能同时考虑节理的分布及连通情况，此法特别适用于裂隙岩体二维渗透张量的研究。由于节理裂隙分布的随机性，数值试验中采用的裂隙产状及几何参数也不能完全的反映岩体内节理的实际分布，由此计算得到的渗透张量存在理论误差。特别地，通过现场压水试验结果来修正模拟结果，保证数值和工程要求相一致。

(3) 该研究通过渗透张量计算方法，考虑了节理的实际分布及连通情况、不同岩层及岩体等级对渗透张量的影响。该法具有计算方便、易于理解的优势，非常适合用于研究裂隙岩体渗透性。

[1] 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.

[2] Sagar Budhi, Runchal Akshai. Permeability of fractured rock: effect of fracture size and data uncertainties[J]. Water Res. Research, 1982,18(2):266-274.

[3] Oda Masanobu. Permeability tensor for discontinuous rockmass[J]. Geotechnique,1985,35(4):483-495.

[4] Oda Masanobu. An equivalent continuum model for coupled stress and fluid flow analysis in jointed rock masses[J]. Water Resources Research, 1986,22(13):1845-1856.

[5] Zhang X, Sanderson D J, Harkness R M, et al. Evaluation of the 2-D permeability tensor for fractured rock masses[J]. International Journal of Rock Mechanics and Mining Science & Geomechanics Abstracts, 1996,33(1):17-37.

[6] Chen S H, Feng X M, Isam S. Numerical estimation of REV and permeability tensor for fractured rock masses by composite element method[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008,32(12):1459-1477.

[7] Rong Guan, Peng Jun, Wang X J, et al. Permeability tensor and representative elementary volume of fractured rock masses[J]. Hydrogeology Journal, 2013,21(7):1655-1671.

[8] 王培涛,杨天鸿，于庆磊，等.基于离散裂隙网络模型的节理岩体渗透张量及特性分析[J].岩土力学,2013,34(S2):448-455.

[9] 杨金保,李 蕾,谢 妮.考虑应力状态的裂隙岩体渗透系数张量计算[J].人民黄河,2013,35(3):131-134.

[10] 车富强,陈益峰，等.基于等效连续介质坝区裂隙岩体渗控效应分析[J].中南大学学报(自然科学版),2014,45(6):1967-1974.

[11] 王俊奇，董 晔，李鹏宇，等.三维离散裂隙网络管单元模型确定岩体渗透张量的尝试[J].水利与建筑工程学报,2016,14(3):47-54.

[12] 张宜虎.岩体等效水力学参数研究[D].武汉：中国地质大学(武汉)，2006.

[13] 张贵科,徐卫亚.裂隙网络模拟和REV尺度研究[J].岩土力学,2008,29(6):1675-1670.

[14] 李杨鹏,薛新华，杨兴国，等.锦屏一级水电站降雨入渗条件下左岸边坡稳定性分析[J].水利与建筑工程学报,2014,12(1):16-20.

[15] 安玉华,王 清.基于三维裂隙网络的裂隙岩体表征单元体研究[J].岩土工程,2012,33(12):3775-3780.

[16] 荔 强.锦屏一级左岸边坡雾化渗流稳定性问题研究[D].北京：清华大学，2010.

Numerical Simulation of Permeability of Rock Mass on the Left Bank of Jinping-I

TIAN Hongyuan1, ZHANG Liang2, LI Hong1

(1.StateKeyLaboratoryofCoastalandOffshoreEngineering,DalianUniversityofTechnology,Dalian,Liaoning116024,China;2.PowerChinaKunmingEngineeringCorporationLimited,Kunming,Yunnan650041,China)

Based on the numerical simulation and finite element analysis software RFPA2D-Flow and single-hole packer tests, the two-dimensional representative element volume and permeability tensor was analyzed for fractured rock mass. According to the geometric characteristics of the joints, such as the length, the density and the yield of the joints. The random fracture networks of Jinping I hydropower station are generated by using Monte Carlo method, numerical simulation of RFPA2D-Flow is adopted to calculate the equivalent permeability coefficient of different size in different directions of fractured rock masses, in order to obtain the representative elementary volume and permeability tensor. Then they are further modified by the single-hole packer test results. The results show that as the model size increases, the permeability of rock gradually begins to stabilize, the size of REV is about 4 times the length of trace. The method is reasonable and effective to determine the permeability tensor and REV of fractured rock masses combining single-hole packer tests and numerical simulation of RFPA2D-Flow.

fracture network model; representative element volume; permeability tensor; Monte Carlo method; RFPA2D-Flow

10.3969/j.issn.1672-1144.2017.03.013

2017-02-14

2017-03-10

国家重点基础研究“973”计划项目(2011CB013503)

李 宏(1969—)，男，辽宁辽阳人，教授，主要从事储层流动数值模拟、边坡节理岩体渗流研究。E-mail:hong.li@dlut.edu.cn

TV139.14

A

1672—1144(2017)03—0067—05