基于等效连续介质的坝区裂隙岩体渗控效应分析

车富强^{1, 2}，陈益峰^{1, 2}，李毅^{1, 2}，张璇^{1, 2}，郑海圣^{1, 2}

(1. 武汉大学 水资源与水电工程国家重点实验室，湖北 武汉，430072；

2. 武汉大学 水工岩石力学教育部重点实验室，湖北 武汉，430072)

摘要：将裂隙样本法和现场测试法结合，先用裂隙样本法初步确定裂隙岩体的渗透张量，获得渗透主值和主方向，再采用现场测试法来修正裂隙岩体的渗透张量。用上述方法对西南某水电站坝区裂隙岩体的渗透张量进行计算，通过对坝区右岸地下水位的数值计算值与钻孔实际水位的比较，验证该方法的可靠性。最后对坝区进行三维有限元渗流分析，验证坝区防渗排水措施的合理性。

关键词：渗透张量；裂隙岩体；裂隙样本法；现场测试法；渗流分析

中图分类号：TV223          文献标志码：A         文章编号：1672-7207(2014)06-1967-08

Seepage analysis of fractured rock mass in a dam foundation based on equivalent continuum model

CHE Fuqiang^{1, 2}, CHEN Yifeng^{1, 2}, LI Yi^{1, 2}, ZHANG Xuan^{1, 2}, ZHENG Haisheng^{1, 2}

(1. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China;

2. Key Laboratory of Rock Mechanics in Hydraulic Structural Engineering (Ministry of Education),

Wuhan University, Wuhan 430072, China)

Abstract: A method was adopted combining the fissure sample measuring method and field test method. First, the permeability tensor of the fractured rock mass was determined preliminarily using the fissure sample measuring method to obtain the principle values and directions, and then it was modified using the field test method. The permeability tensor of fractured rock mass in dam foundation of some southwest hydropower station was calculated with the above method. Through the comparison of the numerical calculation value and the actual value of the groundwater level in the right bank, the reliability of the permeability tensor was verified and the above method was feasible. Then the result of the three-dimensional finite element seepage analysis of the dam area preliminarily has verified the reasonableness of seepage control measures of the dam area.

Key words: permeability tensor; fractured rock mass; fissure sample measuring method; field test method; seepage analysis

裂隙岩体渗流分析的基本模型有等效连续介质模型、离散裂隙网络模型和双重介质模型。等效连续介质模型尽管无法描述裂隙的局部渗流行为，但对裂隙岩体的宏观渗流特征具有较强的描述能力，且计算量较小，因而在大型工程渗流分析中得到了广泛应用。岩体在成岩和地质构造作用下，其不连续面往往成组发育，渗透特性表现出显著的各向异性特征。渗透张量是研究岩体渗透特性各向异性的有效数学工具，岩体渗透张量的确定是进行裂隙岩体渗流分析的基础，目前，确定裂隙岩体渗透张量的方法主要有反演法、现场测试法和裂隙样本法等。反演法是根据钻孔获得的地下水位观测资料，反推岩体的渗透张量，该方法的准确性依赖于动态数据的可靠性及模型选择的合理性。现场测试法是根据野外压水试验及少量抽水试验，通过公式换算来确定岩体的渗透特性，包括单孔压水试验、三段压水试验^[1]等。单孔压水试验的理论基础是各向同性的，所得渗透系数值是无方向的标量，不能反映各向异性，而三段压水试验耗资较大。裂隙样本法是根据裂隙几何测量值(隙宽、间距、产状)的统计来计算渗透张量，该方法简单实用，理论基础完备。但在实际岩体中，除裂隙的空间方位基本能测定外，裂隙的隙宽和间距是很难准确测量的，因隙宽等测量不准带来的计算误差很大，呈指数级^[2]放大。因此，本文将裂隙样本法和现场测试法综合考虑，既能反映裂隙岩体的渗透特性，又具有一定的可操作性。

1  确定裂隙岩体渗透张量的方法

1.1  用裂隙样本法初步确定渗透张量

裂隙岩体的渗透性取决于裂隙的几何性质和空间发育特征。因此，通过对裂隙的空间展布状况(隙宽、间距、产状等)进行测量，可用统计分析方法初步确定岩体的渗透张量^[3]。

假定岩体发育有n组结构面，每组结构面相互平行，第i组结构面的平均间距和平均开度分别为s_i和b_i，Snow^[4]将岩体的渗透张量表示为

          (1)

式中：K为岩体的渗透张量；δ为Kronecher Delta张量；n_i为第i组结构面的单位法向向量。

假定空间直角坐标系x和y轴方向与地理N和W方向一致，则岩体的渗透张量可表示为

 (2)

式中：α_i和β_i分别为第i组结构面的倾角和倾向(从地理N方向顺时针旋转角度)。

根据裂隙渗流的立方定理，裂隙的渗透系数与其开度的二次方成正比，因此裂隙开度的变化会引起渗透特性显著改变。在一般情况下，通过钻孔岩芯取样来量测隙宽，取样后，应力释放为0 MPa，此时裂隙开度随着应力的释放而增大，在荷载作用下，裂隙法向变形为^[5]

                  (3)

式中：u_n为裂隙法向变形；k_n为裂隙的法向刚度; 为有效法向应力。

那么，卸荷前的裂隙初始开度b₀，可采用下式计算：

                  (4)

式中：b为裂隙实测开度。

通过上述方法计算裂隙岩体的渗透张量，需要获得裂隙的法向刚度值，对于缺乏该参数的部分节理组，可以借鉴类似的工程经验，或者他人的室内试验结果。在本文稍后的工程实例中，裂隙开度多集中在0.1 mm左右，而通过查阅相关文献，我们发现Sharifzadeh等^[6]通过室内试验，量测了某裂隙初始开度为0.1 mm时，裂隙开度与法向应力增大的变化关系，根据他们的实验成果，拟合了相应的曲线，如图1所示。

图1  裂隙初始开度随不同法向应力变化图

Fig. 1  Joint initial aperture under different normal stresses

由于渗透张量为对称的二阶张量，按照二阶张量性质，通过求渗透矩阵的特征值和特征向量可求得渗透张量的主值(K₁，K₂和K₃)及对应的主渗透方向。

为了便于对该渗透张量进行修正，取3个渗透主值的几何平均值作为综合渗透系数^[7]，其表达式为

            (5)

1.2  渗透张量的修正

运用单孔压水试验可求得均质各向同性渗透系数。根据野外压水试验得到的岩体透水率q与单位吸水量ω的关系为：ω=0.01q。当q＜10 Lu(1 Lu=1 L/(min·MPa·m))时，可根据巴布什金公式近似计算岩体的渗透系数K^[8]：

             (6)

式中：ω为单位吸水量(L/(min·m²))；L为试验段长度(m)；a为与试验段位置有关的系数，当试验段底与下部隔水层的距离大于L时取0.66，反之取1.32；r₀为钻孔半径(m)。

该方法得出的渗透系数值具有尺寸效应，不能代表某一工程地质单元岩体的平均透水性。前人研究表明^[9]，这种情况下算术平均值比岩体渗透性实际值高，几何平均值更接近岩体渗透性实际值，故取几何平均值为岩体等效渗透系数^[10-13]：

             (7)

式中：K_ω为等效渗透系数；K_i(i=1，2，…，n)为各试验段的计算结果。为了对渗透张量进行修正，定义修正系数为m=K_ω/K₀。

将m代入由裂隙样本法求得的渗透张量中，得到裂隙岩体的修正渗透张量主值为

        (8)

修正渗透张量K_m的主方向与K的主方向一致。

2  各向异性渗流数值模拟

根据达西定律和质量守恒原理，各向异性稳定渗流控制微分方程为

             (9)

式中：h为水头；k_ij(i，j=1，2，3)为渗透张量的分量。上式应满足如下边界条件。

(1) 水头边界条件：

              (10)

式中：Γ_h为已知水头边界；为已知水头。

(2) 流量边界条件：

          (11)

式中：Γ_q为已知流量边界；q_n为边界流量(渗入为负，溢出为正)；n为边界面单位外法向向量。对于隔水边界，q_n=0。

(3) 自由面边界条件(Γ_f)：

              (12)

式中：Γ_f为自由面边界。

(4) 溢出面边界条件(Γ_s)：

         (13)

式中：Γ_s为溢出面边界。

为了评价坝区防渗排水措施的合理性，需要对防渗排水渗控系统的渗流控制效应进行精细模拟。对于含复杂渗控结构的渗流问题，有限元数值模拟的难度主要体现在渗流场的精细模拟和排水孔幕的精细模拟这2个方面。Chen等^[14-15]提出的排水子结构、Signorini型变分不等式和自适应罚Heaviside函数相结合的方法(简称SVA方法)很好地解决了含密集排水孔等复杂渗控结构的渗流问题。子结构法用来模拟排水孔的渗流行为，简化有限元建模；变分不等式法用于确定自由面与排水孔相交时的岀渗点和自由面的位置；而自适应罚Heaviside函数用于保证算法的鲁棒性和解的稳定性。

3  工程实例

卡拉水电站位于雅砻江干流中游河段，为雅砻江干流中游两河口至卡拉河段水电规划一库六级开发之第6个梯级水电站。卡拉水电枢纽主要由碾压混凝土重力坝、坝身泄洪建筑物及右岸地下引水发电系统等组成。坝顶高程1 995.0 m，正常蓄水位1 987.0 m，建基面最低高程1 869.0 m，最大坝高126.0 m。坝址段河道顺直，河谷狭窄，呈基本对称的“V”型峡谷。坝址区河谷走向为S34°E，地层呈单斜构造，岩层总体倾向左岸，左岸为反向坡，右岸为顺向坡，岩层总体产状为走向N20°~40°W，倾角NE∠45~70°。两岸坝顶高程以下及河床出露的岩体为杂谷脑组(T_3Z²)砂质板岩、变质砂岩互层夹含炭质板岩及大理岩。根据地表地质测绘、实测剖面和平洞资料，左岸共统计到8 363条节理，右岸共统计到10 477条节理。

3.1  渗透张量的计算

下面以钻孔Zk83和Zk82附近区域为例来说明渗透张量计算过程。

钻孔ZK83主要揭露了1组节理，其产状为走向N20°W，倾角SW∠15°。揭露的节理平均间距为3.76 m，埋深范围为5~46 m，这段范围内的法向应力值约为0.67 MPa，由图1可知节理开度约为0.063 mm。

钻孔ZK82主要揭露了2组节理，第1组节理走向为N20°W，倾角为SW∠15°，平均间距为3.08 m；第2组节理走向为N89°E，倾角为SE∠9°，平均间距为3.90 m。埋深15~46 m，这段范围内的法向应力约为0.82 MPa，由图1知节理开度约为0.052 mm。

使用式(2)和(5)初步确定渗透张量的主值和综合渗透系数，计算结果如表1所示。

根据钻孔ZK83和ZK82的压水试验结果，使用式(6)计算各试段的渗透系数，取其几何平均值作为等效渗透系数，结果如表2和表3所示，并求出相应的修正系数。修正后的渗透张量结果如表4所示。

用上述方法计算坝区其他部位渗透张量，根据计算结果，可以将渗透张量相差不大的若干分区合并为一个分区。则坝区最后不同分区的渗透张量如表5所示。

3.2  渗透张量修正值的验证

通过对坝区右岸地下水位进行数值计算分析，与钻孔实际水位进行比较，来验证用上述方法计算渗透张量的可靠性。选取坝址区右岸建立三维整体有限元模型(如图2所示)，共划分单元700 798个，节点299 075个，模型计算范围为垂直坝轴线方向(顺河向)2 085 m，平行轴线方向1 190 m，最大高度1 237 m。河道内的表面结点为已知水头边界，上、下游侧面边界和底部边界取隔水边界，模型上表面除地表水淹没区之外的区域设为潜在溢出边界。计算结果如表6所示。

表1  初步确定渗透张量

Table 1  Preliminary determination of permeability tensor

表2  钻孔ZK83压水试验成果汇总表

Table 2  Results of water-pressure test in hole ZK83

表3  钻孔ZK82压水试验成果汇总表

Table 3  Results of water-pressure test in hole ZK82

表4  修正渗透张量

Table 4  Modified permeability tensor

表5  坝区不同分区渗透张量

Table 5  Permeability tensor of different partitions in a dam foundation               cm/s

从表6可以看出：钻孔实测水头与计算水头(修正前)的最大误差为6.02 m，平均误差为2.68 m；钻孔实测水头与计算水头(修正后)最大误差为5.63 m，平均误差为2.5m，渗透张量修正后，最大误差和平均误差均下降，而且修正后的误差基本满足精度要求，说明通过本文确定渗透张量的方法是可行的。

3.3  渗控效应分析

确定了坝区岩体的渗透张量后，就可以对坝区的防渗排水体系进行精细模拟，初步复核坝基及坝体防渗排水措施布置是否合理。

大坝防渗排水工程由防渗帷幕、排水孔幕和排水廊道等组成。在上游侧，大坝左岸分别在1 880，1 895.5和1 930 m高程设置排水廊道；大坝右岸分别在1 887，1 903和1 932 m高程设置排水廊道。为降低大坝坝基扬压力，在上游帷幕的下游设置一排主排水孔，在坝底中部设置一排纵向辅助排水廊道，并利用下游帷幕灌浆廊道作为纵向排水廊道，分别设置一排辅助排水孔。各层排水廊道通过排水孔幕连接，形成排水系统，孔距3 m。

为满足计算分析的需要，建立了坝区整体三维有限元模型(如图2)，上、下游侧边界与坝轴线的距离分别为616 m和634 m，均超过最大坝高(126 m)的4倍，因此，在三维渗流分析中，上、下游侧边界可视为隔水边界；模型底部边界为隔水边界。上、下游水库淹没区内的表面结点为已知水头边界；坝基排水孔为定水头边界，其水头值取与之相连的排水廊道的底板高程一致。坝体中的排水孔及排水廊道视为潜在出渗边界。模型左、右侧边界位于山体中，取定水头边界。

表6  钻孔水位实测值与计算值

Table 6  Measured value and calculated value of drilling water level

图2  坝区三维有限元模型

Fig. 2  Three-dimensional finite element model of dam foundation

共划分单元1 518 637个，节点778 259个。

在水库正常运行工况时，即上游正常蓄水位1 987.00 m，下游水位为1 913.96 m，溢流坝段剖面等水头线及自由面图如图3所示，坝基扬压力分布如图4所示；左岸挡水坝段剖面等水头线及自由面如图5所示，坝基扬压力分布如图6所示；右岸挡水坝段剖面等水头线及自由面如图7所示，坝基扬压力分布如图8所示。

图3  溢流坝段剖面等水头线图(单位：m)

Fig. 3  Isoline maps of water head in overflow dam section

图4  溢流坝段坝基扬压力图

Fig. 4  Uplift pressure maps at dam base in overflow dam section

图5  左岸挡水坝段等水头线图(单位：m)

Fig. 5  Isoline maps of water head in left bank water retaining dam section

图6  左岸挡水坝段坝基扬压力图

Fig. 6  Uplift pressure maps at the dam base in left bank water retaining dam section

图7  右岸挡水坝段等水头线图(单位：m)

Fig. 7  Isoline maps of water head in right bank water retaining dam section

图8  右岸挡水坝段坝基扬压力图

Fig. 8  Uplift pressure at dam base in right bank water retaining dam section

由图3~8可知：对于溢流坝段来说，在上游侧，经上游排水孔幕的排水降压作用之后，坝体内的孔隙水压显著降低，渗流自由面在坝体内部自上游向纵向辅助排水廊道上方急剧下降，下游坝坡不存在渗流溢出点；在下游侧，因河床水位较高，渗流自由面在坝体内部自下游向纵向辅助排水廊道上方急剧下降，说明下游排水孔幕也起到了良好的排水降压作用。坝体中渗流自由面的最低位置出现在纵向辅助排水廊道上方，并穿过了与之相连的纵向辅助排水孔，表明纵向辅助排水也起到了一定的排水降压作用。水头等值线在坝体上游侧和排水孔附近比较密集，坝体下游侧和坝基中比较稀疏。

对于非溢流坝段，在上游侧，渗流经上游排水孔幕的排水降压作用之后，坝体内的孔隙水压显著降低，渗流自由面在坝体内部自上游向下游排水廊道附近急剧下降，下游坝坡不存在渗流溢出点；在下游侧，下游排水孔幕也起到了良好的排水降压作用，渗流自由面也在坝体内部自下游向下游排水廊道附近急剧下降。水头等值线在坝体上游侧和排水孔附近比较密集，坝体下游侧和坝基中比较稀疏。

溢流坝段坝基扬压力最大值为120.0 m，出现在坝踵处；左岸挡水坝段坝基最大扬压力为95.0 m，出现在坝踵处；右岸挡水坝段坝基最大扬压力为88.0 m，出现在坝踵处；各个坝段扬压力变化的一般规律是：在上游排水体系附近急剧下降，然后在坝体中部比较平缓。下游排水孔下游侧的扬压力分布受下游河床水位控制。

三维计算得到坝体和坝基渗漏总流量为1 021.53 m³/d，折算成单宽流量为4.02 m³/d。流入各廊道的渗漏量如表7所示。从计算结果上看，坝区渗漏量已得到有效控制。

由计算结果可知：坝区防渗排水体系的排水降压效果良好，特别是上游帷幕及排水孔幕的作用十分显著，坝区的渗漏量也得到有效控制，初步说明该坝区防渗排水体系合理。

表7  各廊道渗漏量

Table 7  Seepage flow rate of galleries   m³/d

4  结论

(1) 裂隙岩体渗透张量的确定是进行渗流分析及渗控设计的基础。采用裂隙样本法和现场测试法相结合的方法，既能较好地反映裂隙岩体的各向异性，又能减少因裂隙开度和间距测量不准所带来的误差，具有可操作性。

(2) 对坝区右岸地下水位进行模拟计算，通过钻孔实测水头和计算水头的比较，验证了裂隙样本法和现场测试法相结合方法在工程上的可行性。

(3) 用计算所得的岩体渗透张量，对坝区进行三维有限元渗流分析，验证了坝区防渗排水措施的合理性。

参考文献：

[1] He J, Chen S, Shahrour I. Back analysis of equivalent permeability tensor for fractured rock masses from packer tests[J]. Rock Mechanics and Rock Engineering, 2011, 44(4): 491-496.

[2] 荣冠, 周创兵, 王恩志. 裂隙岩体渗透张量计算及其表征单元体积初步研究[J]. 岩石力学与工程学报, 2007, 26(4): 740-746.

RONG Guan, ZHOU Chuangbing, WANG Enzhi. Preliminary study on permeability tensor calculation of fractured rock mass and its representative elementary volume[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(4): 740-746.

[3] 田开铭, 万力. 各向异性裂隙介质渗透性的研究与评价[M]. 北京: 学苑出版社, 1989: 11-12.

TIAN Kaiming, WAN Li. Research and evaluation of anisotropic fracture permeability[M]. Beijing: National Academy Press, 1989: 11-12.

[4] Snow D T. Anisotropic permeability of fractured media[J]. Water Resources Research, 1969, 5(6): 1273-1289.

[5] 陈益峰, 周创兵, 盛永清. 应变敏感的裂隙及裂隙岩体水力传导特性研究[J]. 岩石力学与工程学报, 2006, 25(12): 2441-2452.

CHEN Yifeng, ZHOU Chuangbing, SHENG Yongqing. Strain-dependent hydraulic conductivity for single rock fracture and fractured rock mass[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(12): 2441-2452.

[6] Sharifzadeh M, Mitani Y, Esaki T. Rock joint surfaces measurement and analysis of aperture distribution under different normal and shear loading using GIS[J]. Rock Mechanics and Rock Engineering, 2008: 299-323.

[7] LI Ping, LU Weng, LONG Yuqiao, et al. Seepage analysis in a fractured rock mass: The upper reservoir of Pushihe pumped-storage power station in China[J]. Engineering Geology, 2008, 97(1): 53-62.

[8] 刘佑荣, 唐辉明. 岩体力学[M]. 北京: 中国地质大学出版社, 1999: 120.

LIU Yourong, TANG Huiming. Rock mass mechanics[M]. Beijing: China University of Geosciences Press, 1999: 120.

[9] 切尔内绍夫 C H. 水在裂隙网络中的运动[M]. 盛志浩, 田开铭译. 北京: 地质出版社, 1987: 89-95.

Chernyshov C H. Water moment in the fracture network[M]. SHEN Zhihao, TIAN Kaiming, transl. Beijing: Geological Press, 1987: 89-95.

[10] Lee C H, Deng B W, Chang J L. A continuum approach for estimating permeability in naturally fractured rocks[J]. Engineering Geology, 1995, 39(1/2): 71-85.

[11] 黄义, 张引科. 张量及其在连续介质力学中的应用[M]. 北京: 冶金工业出版社, 2002: 115-124.

HUANG Yi, ZHANG Yin-ke. Tensor and its application to continuum mechanics[M]. Being: Metallurgical Industry Press, 2002.

[12] 王鹏, 蔡美峰, 周汝弟. 裂隙岩体渗透张量的确定和修正[J]. 金属矿山, 2003(8): 5-7.

WANG Peng, CAI Meifeng, ZHOU Rudi. Estimation and Modification of permeability tensor of fractured rock mass[J]. Metal Mine, 2003(8): 5-7.

[13] 纪伟. 裂隙岩体渗流等效性及数值模拟研究[D]. 南京: 河海大学水利水电工程学院, 2005: 20-23.

JI Wei. Study on equivalent characteristic and numerical simulation of seepage in fractured rock[D]. Nanjing: Hohai University. College of Water Conservancy and Hydropower Engineering, 2005: 20-23.

[14] CHEN Y F, ZHOU C B, ZHENG H. A numerical solution to seepage problems with complex drainage system[J]. Computers and Geotechnics, 2008, 35(3): 383-393.

[15] 陈益峰, 周创兵, 郑宏. 含复杂渗控结构渗流问题数值模拟的SVA方法[J]. 水力发电学报, 2009, 28(2): 89-95.

CHEN Yifeng, ZHOU Chuangbing, ZHENG Hong. A SVA method for seepage problems with complex control measures[J]. Journal of Hydroelectric Engineering, 2009, 28(2): 89-95.

(编辑  何运斌)

收稿日期：2013-06-14；修回日期：2013-08-13

基金项目：国家重点基础研究发展计划(“973”计划)项目(2011CB013503)；国家自然科学基金资助项目(51222903，51079107)

通信作者：陈益峰(1974-)，男，福建漳平人，教授，博士生导师，从事岩土多场耦合、岩土渗流、变形与稳定分析的研究工作；电话：027-68774295；E-mail：csyfchen@whu.edu.cn