裂隙岩质边坡非饱和降雨入渗特征分析

杨秀竹, 康 镜, 刘正夫, 雷金山

(1.中南大学 土木工程学院, 湖南 长沙 410075; 2.江苏省交通科学研究院, 江苏 南京 210004)



裂隙岩质边坡非饱和降雨入渗特征分析

杨秀竹1, 康 镜1, 刘正夫2, 雷金山1

(1.中南大学 土木工程学院, 湖南 长沙 410075; 2.江苏省交通科学研究院, 江苏 南京 210004)

[目的] 研究降雨入渗过程中，裂隙岩质边坡渗流场和暂态饱和区的变化规律，为裂隙岩体降雨入渗规律、工程技术及水土保持方面的研究提供支持。 [方法] 基于非连续裂隙模型，采用Matlab编制裂隙岩体渗透张量的计算程序；基于等效连续介质渗流模型，采用FLAC3D内置的FISH语言编制程序进行非饱和渗流分析。 [结果] 依托实际工程分析得出了裂隙岩体的渗透张量，在降雨过程中随着降雨强度的增加或者降雨时间的延长，边坡表层区域零孔隙水压力面不断向内部发展，暂态饱和区的范围不断扩大，且暂态饱和区的出现发生在降雨强度大于边坡的入渗率的情况。 [结论] 采用非连续裂隙网络模型和等效连续介质渗流模型对裂隙岩质边坡降雨入渗分析能较好反映裂隙岩体各向异性和非饱和渗流的特点。

边坡; 降雨; 渗透张量; 等效连续介质; 非饱和渗流

文献参数： 杨秀竹, 康镜, 刘正夫, 等.裂隙岩质边坡非饱和降雨入渗特征分析[J].水土保持通报，2016,36(4)：143-147.DOI:10.13961/j.cnki.stbctb.2016.04.026

当岩体破碎时，岩质边坡易受风化降雨作用的影响。随着降雨过程的发展，雨水入渗会引起边坡非饱和区的负孔隙水压力发生较大变化，该部分岩体的自重增加，并且使岩体软弱结构面强度下降，引起边坡的失稳。从而导致引起边坡的失稳的可能性增加，为工程及群众生命财产安全造成较大威胁。因此，对降雨时裂隙岩体中的渗流特征以及渗流场的变化规律进行研究，对降低安全风险，防止滑坡灾害发生有较为重要的意义。渗流介质是影响边坡降雨入渗重要因素之一，刘子振[1]研究了降雨对非饱和黏土边坡的入渗规律，但与研究土质边坡降雨入渗不同，裂隙岩体渗流具有的显著各向异性特点是不能忽略的。渗透系数影响分析结果的正确性，张升堂等[2]用年降雨总量与产流降雨总量的统计推求降雨入渗系数，认为产流降雨总量能比较清晰地反映流域治理程度对降雨入渗的影响，高鹏等[3]采用人工模拟降雨的试验方法研究了不同降雨强度对土壤入渗的影响，建立了土壤水分入渗经验模型。但还须从渗流介质及其结构特点的角度进一步研究。分析裂隙岩质边坡降雨入渗的模型主要有：等效连续介质模型、离散裂隙网络模型和双重介质模型[4]。离散裂隙网络模型认为岩块不能导水而忽略岩块的导水性，但当考虑裂隙岩体非饱和渗流时，岩块的导水性是不能忽略的；双重介质模型中裂隙网络和连续介质岩块耦合模型中的水交换量公式直接影响模型计算的准确性，但水交换量难以准确确定；当岩体裂隙节理等结构面充分发育时，等效连续介质模型采用渗透张量矩阵能较好地模拟裂隙岩体渗流各向异性的特点，并且该模型只需知道裂隙的几何信息的统计值，而对其准确位置并不敏感。因此，本研究采用非连续裂隙网络模型确定裂隙岩体渗透张量并通过等效连续介质渗流模型对裂隙岩质边坡非饱和降雨入渗进行分析，其研究成果可为裂隙岩体降雨入渗规律研究，工程技术及水土保持研究提供支持和参考。

1　等效连续介质渗流原理

(1)

需要指出的是对于岩体而言，将裂隙等效为连续介质模型必须对裂隙岩体进行分析，判断其是否满足等效条件。Long[5]指出了等效裂隙岩体满足：表征单元体REV存在；存在等效渗透张量，且裂隙岩体方向渗透率在极坐标下能形成椭圆，即该等效渗透张量是对称，才可以运用等效连续介质模型。

2　渗透张量的确定

渗透张量是等效连续介质渗流计算模型中的一个关键参数，它直接反映了裂隙岩体非均质各向异性的特点。采用Matlab编制基于裂隙网格模型的渗透张量计算程序包括：运用蒙特卡洛法生成裂隙，并对裂隙进行筛选，最终得到有效裂隙；计算不同大小分析区域中的渗透张量，并给出最小REV和对应的渗透张量[6]。

2.1裂隙生成

生成裂隙采用的蒙特卡洛方法是一种依据大数定律的统计试验方法，它采用随机数进行统计试验，产生符合某种随机分布函数的随机数。在裂隙岩体的渗流分析中，可以根据裂隙几何参数的概率模型[7]，在满足给定的空间分布规律的情况下，随机生成裂隙的长度、方向角、中心点位置。采用蒙特卡洛法生成裂隙的具体流程如图1所示。

图1　采用蒙特卡洛法生成裂隙流程图

2.2渗透张量计算

在生成有效裂隙网格的基础上，渗透张量的计算流程如图2所示。

图2　渗透张量计算流程图

2.3渗透张量计算验证

某岩体的生成域为30 m×30 m，分析域为15 m×15 m，共有2组裂隙，裂隙的中心点服从均匀分布，方向角和迹长均服从正态分布，隙宽为常数；具体裂隙参数详见表1。

表1　试验岩体裂隙参数

对裂隙进行筛选之后的有效裂隙表征单元体尺寸为 ，在该尺寸内可认为裂隙是贯通的，该裂隙渗透张量系数为：

当采用裂隙测量法时，计算出渗透系数张量：

对比两种算法计算出来的渗透系数张量，发现两者很接近，说明本研究对渗透张量的计算模块是合理可靠的。

3　FLAC3D饱和-非饱和降雨入渗分析功能的二次开发

FLAC3D是目前国内外岩土工程界应用广泛的有限差分软件。但由于其在非饱和渗流计算中会将负孔隙水压区的饱和度强制设为1[8]，这与实际情况不符，部分学者改进FLAC3D饱和非饱和渗流分析模块[9]，取得了一些进展。因此，为了利用FLAC3D软件平台进行裂隙岩体非饱和渗流分析，利用内置的FISH语言对FLAC3D渗流模块的二次开发，在迭代过程中根据饱和度不断地修正非饱和区的渗透张量，求出比传统渗流分析更准确的渗流场。对FLAC3D的二次开发可分为两个部分，即饱和—非饱和的渗流分析功能，以及降雨入渗及出渗边界功能的FISH函数编制[10]。

进行饱和—非饱和渗流分析时，在FLAC3D的计算过程中，自动根据上一计算增量时间步的结果，调整非饱和区的渗透张量，从而实现FLAC3D的非饱和渗流计算。本研究采用Genuchten[11]提出单裂隙毛管压力—饱和度关系曲线和相对渗透率—饱和度关系曲线：

(7)

(8)

式中：Pc——毛管压力； P0——进气值； Se——水相的有效饱和度； m,n——经验参数； kr——相对渗透率。

进行降雨入渗和出渗边界分析时，降雨入渗率取饱和渗透系数和降雨强度中的较小值。降雨入渗为动态的边界条件，在降雨过程的分析中，需要在每一个时间步对边坡表层节点的孔隙水压进行判断，若孔隙水压小于0，节点应采用流量边界；而当孔隙水压大于或等于0时，则需将该节点由流量边界修改为已知水头边界。

4　算例分析

4.1工程概况

长沙大冰雪世界矿坑酒店项目位于长沙市岳麓区坪塘镇，边坡主要为岩质边坡，岩石为古生代泥盆系中统棋梓桥组浅灰色—深灰色薄层—厚层灰岩，厚度大于210m。边坡大部分基岩裸露，边坡岩体表面裂隙性溶蚀风化现象严重，沿断层、裂隙及层面等结构面溶蚀风化现象较普遍，风化裂隙发育。据1960—2013年长沙市气象站资料统计：平均降雨量1 394.6mm，最大年降雨量1 751.2mm(1998年)，最小年降雨量708.8mm(1953年)，最大月降雨量515.3mm，最小月降雨量1.2mm，最大日降雨量192.5mm。

4.2渗透张量求解

研究区域边坡岩体主要发育4组节理裂隙，将真实的裂隙节理进行转换，得到其在特征截面上的迹线长度和方向角，运用蒙特卡洛法随机生成二维裂隙网络所需的裂隙几何参数统计详见表2。

将以上数据输入Matlab程序中，便可在 的裂隙生成域内生成二维裂隙网络模型(裂隙网络图略)。

表2　岩体裂隙几何参数统计

为确定裂隙岩体的表征单元体(REV)尺寸，在上述全局域裂隙网络的基础上，对不同分析域尺寸下的渗透张量和对应的均方误差(RMS)进行试算。分析域每次旋转30°，共旋转12次，并赋予相同的边界条件；通过程序计算，现将不同尺寸下的RMS整理如图3所示。

Johan[12]根据计算的RMS判断拟合的椭圆质量，认为只有当 时，各个方向渗透系数拟合的椭圆质量较好，采用等效连续介质模型计算时不会出现大的误差。此时，该区域的尺寸可以视为该裂隙岩体的表征单元体(REV)。当分析域尺寸小于5 m时，相应的RMS大于0.2，其相应拟合的椭圆质量较差；当分析域尺寸在5 m以上时，RMS小于0.2(稳定在0.16附近)，拟合的椭圆质量较好，表征单元体存在。参照

图3，表征单元体尺寸取为6 m×6 m。

通过裂隙网络水力学的方法求出6 m×6 m分析域内裂隙岩体的沿梯度方向的渗透系数详见表3。

图3　不同分析域尺寸下的均方误差值(RMS)变化

角度/(°)渗透系数/(m·s-1)椭圆半径(1-m/s)角度/(°)渗透系数/(m·s-1)椭圆半径(1/m/s)05.51E-06426.191805.51E-06426.19307.07E-06376.152107.07E-06376.15609.44E-06325.472409.44E-06325.47901.36E-05270.752701.36E-05270.751201.08E-05304.683001.08E-05304.681506.55E-06390.753306.55E-06390.75

通过最小二乘法将上表中的数据在极坐标下拟合成渗透椭圆。拟合的椭圆的长短轴计算出渗透张量主值K1=1.33×10-5m/s，K2=5.54×10-6m/s，拟合的渗透椭圆的长轴方向与极坐标的0°方向夹角为93.8°, 小于0.2，故可用这个渗透椭圆求出渗透张量来对裂隙网络模型的渗透性进行描述。考虑到裂隙部分填充、局部非连通以及裂隙测量开度变化等因素影响，计算并修正后的渗透张量(m/s)为：

4.3降雨渗流数值计算

根据实际工程地质条件，建立边坡降雨渗流计算模型，X坐标在坡脚处沿水平方向指向坡内，Z方向在坡脚处沿竖直方向向上为正，计算模型地下水位左侧距模型底部为18 m，右侧为9 m。模型中的人工填土、残积粉质黏土、灰岩Ⅱ均按地勘报告中渗透系数选取，灰岩Ⅵ渗透性由上文分析确定。计算模型网格数为869，节点个数为1 789个。计算模型如图4所示。

4.3.1降雨持续时间对暂态饱和区的影响为分析降雨持续时间对零孔隙水压面的影响，计算边坡在降雨强度为300 mm/d情况下，降雨持续时间为1～4 d零孔隙水压位置的变化图，将数据导入TECPLOT，绘制不同降雨时间下零孔隙水压面。分析得出，在降雨强度为300 mm/d时，随着降雨时间的延长，灰岩Ⅵ表层区域和断层的零孔隙水压面向内部发展的范围不断扩大，暂态饱和区的范围不断扩大，且相比于灰岩Ⅵ表层区域断层处受到降雨持续时间的影响更大。但上覆土层的零孔隙水压面增加区域有限。降雨时间在1～4 d范围内，对边坡内部的浸润面的影响不大。

图4　边坡降雨渗流计算模型

4.3.2降雨强度对暂态饱和区的影响为了分析不同降雨强度对孔隙水压的影响，分别对降雨持续4 d降雨强度为100，200，300 mm/d这3个工况下的边坡渗流进行计算，将计算结果导入TECPLOT，绘制出的不同情况下零孔隙水压面。分析得出，降雨时间为4 d时，随着降雨强度增大，灰岩Ⅵ表层区域的零孔隙水压面向内部发展的范围不断扩大，暂态饱和区的范围不断扩大，但上覆土层的零孔隙水压面增加区域有限。当降雨强度为100 mm/d时，边坡表面刚开始出现零孔压，若降雨强度小于边坡入渗率，雨水将全部入渗。当降雨强度为200 mm/d时，断层处还未出现零孔隙水压面，说明断层处仍处于非饱和状态，但随着降雨强度的增大，当降雨强度为300 mm/d时，断层F1区域内出现零孔隙水压，达到了饱和且影响范围远大于灰岩Ⅵ。在给定降雨时间为4 d时，降雨强度在100～300 mm/d范围内，对边坡内部的浸润面的影响不大。

5　结 论

(1) 利用MATLAB软件，编写了基于蒙特卡洛方法的随机生成裂隙程序和渗透张量计算程序。基于实际工程情况，通过模拟计算，得到了其裂隙分布，并计算出了表征单元体尺寸和渗透张量。

(2) 在模拟裂隙岩体的降雨入渗过程中引入渗透张量，并利用FISH语言编写了非饱和入渗分析的修正模块，完善了FLAC3D的各向异性非饱和渗流分析功能，为利用FLAC3D进行复杂三维裂隙岩质边坡降雨入渗研究提供了基础。

(3) 通过实际工程的研究发现，降雨强度大于边坡的入渗率，可能导致暂态饱和区的出现。随着降雨强度的增加或者降雨时间的延长，边坡表层区域零孔隙水压力面不断向边坡内部发展，暂态饱和区的范围不断扩大，但降雨过程中对边坡内部的浸润面的影响较小。

[1]刘子振.持续降雨入渗非饱和黏土边坡失稳机理及其应用研究[D].甘肃 兰州:兰州大学,2014.

[2]张升堂,拜存有,万三强.黄土高原水土保持措施强化降雨入渗分析及灰色预测[J].水土保持通报,2004,24(2)：29-33.

[3]高鹏,穆兴民,刘普灵,等.降雨强度对黄土区不同土地利用类型入渗影响的试验研究[J].水土保持通报,2006,26(3)：1-5.

[4]张有天.岩石水力学与工程[M].北京：中国水利水电出版社,2005.

[5]Long J C S, Remer J S, Wilson C R, et al. Porous media equivalents for networks of discontinuous fractures[J]. Water Resources Research, 1982,18(3):645-658.

[6]刘海军.基于蒙特卡罗法的岩体裂隙网络模型及渗透张量的研究[D].黑龙江 哈尔滨：哈尔滨工业大学,2011.

[7]刘晓丽,王恩志,王思敬.裂隙岩体表征方法及岩体水力学特性研究[J].岩石力学与工程学报,2008,27(9):1814-1821.

[8]陈育民. FLAC/FLAC3D基础与工程实例[M].北京：中国水利水电出版社,2009.

[9]李毅,伍嘉,李坤.基于FLAC3D的饱和—非饱和渗流分析[J].岩土力学,2012,33(2):617-622.

[10]蒋忠明,熊小虎,曾铃.基于FLAC3D平台的边坡非饱和降雨入渗分析[J].岩土力学,2014,35(3):855-861.

[11]Genuchten M T V. A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils[J]. Soil Science Society of America Journal, 1980,44(5):892-898.

[12]柴军瑞,徐维生.大坝工程渗流非线性问题[M].北京：中国水利水电出版社,2010.

Unsaturated Seepage Analysis of Fractured Rock Slope Under Rainfall Condition

YANG Xiuzhu1, KANG Jing1, LIU Zhengfu2, LEI Jinshan1

(1.SchoolofCivilEngineering，CentralSouthUniversity,Changsha,Hu’nan410075,China; 2.JiangsuTransportationInstituteCo.Ltd.,Nanjing,Jiangsu210017,China)

[Objective] Exploring the variation process of seepage field and transient saturated zone of fractured rock slope under the influence of rainfall to provide support for the study of the law of rainfall infiltration, engineering technology and water conservation. [Methods] Based on discontinuous fracture network model, a Matlab program was designed to calculate the permeability tensor. Based on equivalent continuous medium model, a program using FISH language was designed for unsaturated seepage analysis. [Results] Permeability tensor of fractured rock mass was obtained through the engineering case. With the increase of rainfall intensity and rainfall duration, the surface of zero pore water pressure will extend to the inside of the slope and the scope of transient saturated zone will expand. Moreover, the appearance of transient saturated zone will occur in the case when rainfall intensity is greater than the the infiltration rate of the slope. [Conclusion] The characteristics of anisotropic and unsaturated seepage in fractured rock mass can be better elucidated by using discontinuous fracture network model and equivalent continuous medium model and the models can be used to analyze the rainfall infiltration of fractured rock slope.

slope; rainfall; permeability tensor; equivalent continuous medium; unsaturated seepage

2015-02-22

2016-03-19

湖南省发改委项目“湖南城际轨道交通沿线隐伏型岩溶注浆充填材料应用研究”(湘发改高技[2012]1493号)

杨秀竹(1972—),女(汉族),山东省莱州市人,博士,副教授,主要从事岩石力学、隧道与地下工程的教学科研工作。E-mail:xzyang1661@126.com，134812270@csu.edu.cn。

雷金山(1973—),男(汉族),湖南省湘乡市人,博士,高级工程师,主要从事地质灾害治理的研究工作。E-mail:jslei8522@126.com。

A

1000-288X(2016)04-0143-05

TV223.6