水库下伏采空区渗流场耦合模拟分析

杜子建 高永涛 邓代强 韩浩亮

(北京科技大学土木与环境工程学院,北京 100083)

目前,国内外学者对天然岩体的水力学特性研究较多,而对于扰动后的岩体渗流特性却研究较少,例如,对于目前比较常见的地下采空区,矿层开采之后破坏了围岩的应力状态,造成应力的重新分布,导致围岩变形、破坏,改变了天然岩体的裂隙分布,岩体的渗透性也会随之改变,采空区范围的渗流场与天然岩体必然有较大区别[1,2].特别当采空区位于水体渗流区域,如本文将提到的水库下伏采空区,水库水体会沿着新产生的采动裂隙流动,并在采空区裂隙带形成水力通道,水体对裂隙体渗透作用,以及对裂隙缝中破碎体以及充填物长期冲蚀,会造成裂隙体的强度进一步削弱,引起岩体应力场的环境变迁,导致裂隙岩体的渗透变形,增大岩体的渗透性,在两者长期相互作用下,如不采取措施必然会导致较严重的后果.

1 扰动岩体应力场对渗流场影响

扰动岩体重新分布的应力场改变了裂隙的位移及结构面的隙宽,进而影响地下岩体整体渗透特性,应力场对渗流特性的影响主要表现为岩体渗流场渗透系数k的改变[3].

假设裂隙岩体介质统计为规则和均质的单纯裂隙岩体介质,对于布有N个方向裂隙组成的导水裂隙网络,根据渗透性的可叠加原理,将每组渗透裂隙面的渗透张量进行叠加,则有几何参数表达的总渗透张量公式[4-6],裂隙岩体渗透张量K的数学模型为

式中,ei为第i组裂隙的等效水力隙宽;λi为第i组裂隙的平均线密度;nxi、nyi、nzi为第i组裂隙面法向的方向余弦.

在三向应力作用下,垂直于裂隙面方向的正应力对渗流特性起着至关重要的作用[7-8],Snow(1968)、Louis(1974)、Kilsall(1984)等通过试验研究,相继得出了裂隙岩体渗透系数k随正应力的经验关系式.以此为基础,结合大量勘探经验成果,得出扰动岩体应力场改变环境下裂隙岩体渗透系数k随正应力σ变化的经验关系式:

(1)当正应力σ变化范围比较小时,经验关系式为直线衰减形式

式中,k′为σ=σ0时裂隙岩体的渗透系数,λ为待定系数,由试验确定.

(2)当正应力σ变化范围比较大时,经验关系式为负指数衰减形式

式中,k′为σ=σ0时裂隙岩体的渗透系数,λ为待定系数,由试验确定.

2 扰动岩体渗流场对应力场影响

目前水荷载的计算方法忽略了渗流场的影响,而以静水压力的形式表示水荷载.实际上,在任何透水介质中,水荷载应以渗透体积力和渗透压力(面力)表示,是和透水介质中渗流场的分布密切相关的,一定的渗流场分布对应一定的水荷载分布;而渗流场分布的变化也将引起水荷载分布的变化.所以,渗流场对应力场的影响是通过改变岩体外荷载(水荷载)而改变岩体的应力场分布的[9].

若已知某时刻覆盖层内的水头分布函数h(x,z),则可求出某作用面上的渗透水压力u分布为

进而可求出覆盖层内的单位土体所受的渗透体积力f分布为

式中,fx、fz为渗透体积力f在x、z方向上的分力.采用有限元数值方法计算时,将上述求得的渗透水压力和体积力转化为单元的结点外荷载而加以考虑.

3 耦合分析方法及主要步骤

渗流场与应力场耦合模型的建立,关键是两场耦合关系的正确建立,耦合的非线性在于渗透系数是孔隙率的非线性函数.岩体应力场对渗流场的影响在于应力场改变了介质本身的结构,介质渗透性能发生改变;渗流场对应力场的影响在于介质中单元体所受的渗透力和外部水荷载在不断改变,引起岩体介质中应力的不断调整.

运用物理环境耦合方法进行分析,均质岩体渗流场与应力场耦合分析的有限元数值解法及迭代步骤[10]如下:(1)首先假设应力场为0,求出均质岩体的渗透系数,读取流体物理环境,运行FLOT RAN流体模块求解渗流场,得出水头分布h;(2)将水头分布h代入式(4)、(5),求解出渗透水压力分布u,渗透体积力f;(3)读取结构物理环境,加载渗透水压力u及渗透体积力f,按有限元法求解应力场,得出应力场的分布以及位移场分布;(4)将应力场载入模型,代入式(2)或式(3),可得出扰动岩体应力场改变环境下裂隙岩体渗透系数k′分布;(5)回到第1步,重复1～4步,迭代直至满足精度要求,最后得出耦合条件下渗流场的水头压力和孔隙水压力.

ANSYS实现物理环境耦合模拟分析的方法流程图如图1所示.

图1 物理环境法耦合分析流程

4 模型建立

4.1 模型的假设条件

①岩体是各向同性弹塑性体,屈服满足Drucker-Prager强度准则,为等效连续介质[11];②孔隙水渗流服从Darcy定律;③岩体的渗透系数k是动态值,渗透系数与孔隙率呈一定函数关系[12].

4.2 数值计算模型尺寸设计

计算模型取地表面为模型顶面,忽略表层少量土体层的影响,以基岩作为计算模型表面.由于采空区为地表50m以下,离水库水平距离为30m,距水库底垂直距离为30m左右,为了保证计算边界大小,并充分考虑水库渗流场的影响,本次建模取模型尺寸为800m×800m×320m,模型实地所在库区范围如图2所示.将原地面取为模型顶面,采空区平均埋深为57.5m;水库深度为20m,长宽取300m×140m,模型及水库三维网格示意图如图3所示.

4.3 模型岩体参数确定

模拟计算采用的岩体物理力学参数见表1.

表1 岩体物理力学参数表

根据现场地质调查和相关研究提供的岩石力学及渗流试验结果,考虑了尺寸效应和地层构造面的影响,对实验得出的岩石各项参数进行相应的调整和简化.根据采空区现场观测,岩体裂隙开口度小于0.001 m,每米平均裂隙密度为6～8条,裂隙平均间距为0.12m,采空区平均日涌水量为920m3,由公式(1)得到采空区等效渗透张量K,从而得到岩体弹性模型折减系数[13]为0.4,即数值模拟中岩体弹性模量取岩石弹性模量的0.4倍,采用岩体的破坏为塑形破坏(Drucker-Prager屈服准则)和拉破坏.

4.4 模型边界条件的确定

模型区域地应力场最大主应力方向为S82°E,在计算域边界上施加东西向应力,即在模型边界上施加边界法向水平构造应力,应力边界条件及其外加载模型如图4所示.

图4 计算模型应力边界条件

渗流场满足水力边界条件

式中,h(M,t)为时刻t点M的测压管水头值;φs(M,t)为边界上给定的已知测压管水头值;M为边界S1上的点.

在水库底部及侧边加载水力边界条件,按水库水体高度加载水头,以水头压力直接加载.

5 下伏岩体渗流模拟结果分析

取采空区X-560剖面水头进行分析,从图5可以看出,采空区渗流达到稳定时,在采空区区域整体上形成以采空区为中心,水头由中心向四周扩散逐渐变大的渗流区,其中以采空区附近水头最低,区域水头最大值仅1.3m,水头基本以采空区为中心呈梯度变化并趋向稳定,距离采空区近,水头变化较大,距离采空区越远,水头变化就越小,水力梯度也逐渐减小.

图5 空区水头压力X-560剖面分析图

取采空区X-550剖面孔隙水压进行分析,从图6中可以看出,空区形成后,在采空区内部整体上孔隙水压力呈垂直梯度分布,最大值达到8567Pa,良好的水力通道在采空区上方形成,孔隙水压在采空区上方附近呈直线下降趋势.采空区下方孔隙水压相对较大,为3191Pa,呈垂直梯度变化,水平层一定范围内孔隙水压基本一致.

图6 空区孔隙水压力X-550剖面分析图

通过对水库下伏采空区渗流场模拟分析,可得出以下结论:(1)水库下伏采空区区域,上层覆岩除了会产生由于开采引起覆岩移动外,由于孔隙水压力的减小,其产生的附加有效应力同样会导致覆岩产生一定量的位移,应引起重视.(2)采空区上部破坏后的覆岩透水性较好,如果不进行及时治理,水库甚至地表水会成为渗流补给的来源,进一步加大渗流的水力梯度和流速,水在流动的过程中会带走裂缝中的充填体,久而久之,会形成良好的水力通道,在采空区与水库之间形成大量水力通道,一旦通道连通水库底部,必然会产生严重后果.(3)建议采取措施:①水库防渗加固处理:对水库可能存在渗流裂隙带进行及时加固处理,并对水库设立水位监测点进行观测,及时发现水位陡减情况,排除险情.②采空区填充处理:经空区渗水量的监测,虽然渗水量呈稳定趋势,但不排除在水流长期渗透冲刷下渗漏有进一步恶化的趋势,建议对该空区采用高水材料进行填充以消除隐患.

6 结 论

扰动后的岩体渗流特性研究主要表现在岩体应力场与渗流场的相互作用研究,本文采用大型有限元分析软件ANSYS,运用物理环境耦合方法,在考虑扰动应力环境对岩体渗流参数的影响下,提出了渗透体积力的加载方程,即将渗透水压力和体积力转化为单元的结点外荷载而加以考虑,分析了渗流场影响下的岩体应力场,并采用迭代有限逼近的方法,对扰动岩体下伏采空区应力场与渗流场进行了耦合分析,得出了库区下伏采空区渗流模拟剖面水头压力和孔隙水压力,对实际工程处理有一定参考作用.

文中扰动岩体应力环境下的渗流特性的影响分析只考虑了岩体渗流场渗透系数k的变化,而影响岩体渗流特性的因素有很多,综合考虑其他因素情况下,扰动岩体渗流场与应力场的耦合机理还有待进一步的研究.

[1]任长吉,黄 涛.裂隙岩体渗流场与应力场耦合数学模型的研究[J].武汉大学学报:工学版,2004,37(2):8-12.

[2]M asanobu Oda.A n Equivalent Continuum M odel for Coupled Stress and Fluid Flow Analysis in Jointed Rock Masses[J].Water Resources Research,1986(13):98-103.

[3]Yang L Z,Huang T.Laboratory Study of Low Permeability Water Bearing Medium in Environmental Pollution Control[A].San Diego,Geotechnical Laboratory,Lehigh University,1996:263-271.

[4]毛昶熙.渗流计算分析与控制[M].北京:水利电力出版社,1988.

[5]许光祥,哈秋舲.裂隙岩体渗流与卸荷力学耦合作用及裂隙排水研究[M].重庆:重庆大学出版社,2003.

[6]杨立中,周 训.深层地下水渗流的研究[M].成都:成都科技大学出版社,1995.

[7]周维桓.高等岩石力学[M].北京:水力水电出版社,1990.

[8]张金才,张玉卓.应力对裂隙岩体渗流影响的研究[J].岩土工程学报,1998,20(2):19-22.

[9]陈 平,张有天.裂隙岩体渗流与应力耦合分析[J].岩石力学与工程学报,1994,13(4):299-309.

[10]曹阿静,许梦国.裂隙岩体渗流与应力耦合有限元分析[J].矿业研究与开发,2003,23(3):5-8.

[11]张社荣,杨璐玲,钟登华.裂隙岩体渗流场分析及其三维有限元程序设计[J].岩土力学,2005,26(1):46-50.

[12]蒋中明,Dashnor Hoxha.核废料地质贮存介质黏土岩的三维各向异性热-水-力耦合数值模拟[J].岩石力学与工程学报,2007,26(3):493-501.

[13]Ouyang Zhihua,Elsworth D.Evaluation of Groundwater Flow into Mined Panels[J].Int J.Rock Mech.M in Sci,1993,30(2):71-79.