渗流-应力耦合下深埋引水隧洞变形稳定性分析

王克忠，唐雨蔷，李伟平，孔令民，秦绍坤

(1.浙江工业大学 建筑工程学院，浙江 杭州 310014；2.浙江省交通规划设计研究院 浙江 杭州 310014；

3.日照市供水工程管理处，山东 日照 276800)



渗流-应力耦合下深埋引水隧洞变形稳定性分析

王克忠1，唐雨蔷1，李伟平2，孔令民3，秦绍坤1

(1.浙江工业大学 建筑工程学院，浙江 杭州 310014；2.浙江省交通规划设计研究院 浙江 杭州 310014；

3.日照市供水工程管理处，山东 日照 276800)

摘要：隧洞开挖过程中，地下水渗流作用引起围岩应力重分布，其对围岩稳定性的影响不容忽视.依托日照市沭水东调引水隧洞工程，基于渗流-应力耦合理论，利用MIDAS-GTS有限元分析程序,对有无渗流作用下隧洞的开挖过程进行仿真分析，得出了隧洞开挖过程中围岩孔隙水压力分布及围岩变形规律.计算结果表明：开挖使一定范围内的围岩孔隙水压力降低，并导致围岩内孔隙水压力呈环形带状分布;渗流-应力耦合下围岩位移空间分布规律与未考虑孔隙水压力时相似，耦合作用对拱顶的竖向位移影响最大.

关键词：引水隧洞；开挖；渗流-应力耦合；孔隙水压力；位移

水在岩体中流动会改变岩体的原始应力状态，同时，岩体应力状态的变化又会影响岩体中水的流动特性,两者相互作用产生渗流-应力耦合现象[1-2].隧道工程施工过程中，渗流-应力耦合作用主要体现在隧道开挖后应力重分布所引起的原始应力场及渗流场的改变.目前，国内外学者对隧道岩体内的渗流-应力耦合分析开展了大量的理论研究，靳晓光等[3]在对越江隧道的渗流-应力耦合分析中，考虑了开挖过程中渗流边界的变化；刘仲秋等[4]对锦屏二级水电站深埋引水隧洞施工到运行的全过程进行了渗流-应力耦合数值模拟分析，并考虑了围岩和衬砌结构渗透系数及变形特性的动态变化；李术才等[5]利用研制的新型流固模型试验系统，结合新型流固耦合相似材料，对海底隧道施工过程中洞壁压力和围岩位移场、渗流场等的变化规律进行了研究.纪佑军等[6]利用渗流-应力耦合数学模型，借助Comsol模拟隧道开挖中不同工况下围岩应力及地面沉降情况.黄涛等[7]借助含水裂隙岩体中地下水渗流场与地应力场之间的耦合作用，对裂隙围岩隧道涌水量预测计算的确定性数学模型方法进行了研究.已有的研究表明，岩体中的渗流场和应力场耦合原理及过程极其复杂，现有的理论还不能完全解释实际工程中所出现的渗流应力及其破坏现象[8-12].笔者依托日照市沭水东调引水隧洞工程，利用MIDAS-GTS有限元分析软件建立数值模型，结合开挖过程中渗流边界的变化，确定合理的渗流-应力耦合数值模拟方法，获取孔隙水压力分布规律，分析开挖过程对渗流场的重分布及其演变规律，进一步获取围岩位移场的演化规律及其特征点的位移值.

1渗流-应力耦合计算原理

在MIDAS-GTS软件中，模型采用有限单元法求解渗流-应力耦合问题[13-14]，其求解步骤如下：

1) 利用渗流控制方程求得水头函数H=H(x, y, z).在计算模型中，假定围岩为均质各向同性连续介质，渗流方式为稳定渗流，服从达西定律，其控制方程为

(1)

式中：H为各计算点的总水头；kx，ky，kz分别为x，y，z各方向的初始渗透系数；Q为流量.

2) 在渗流域内，利用前面得到的水头函数，计算相应的渗透压力和渗透体积力.由水力学原理可知渗透体积力的计算方式为

(2)

(3)

式中：f为渗流体积力的大小；γw为水的容重；fx，fy，fz分别为x，y，z渗透体积力在方向的分力；Jx，Jy，Jz分别为x，y，z单元在方向的水力坡降.

3) 在有限元分析中，将得到的渗透体积力以等效结点荷载的方式作用于岩体上，求解岩体的位移场和应力场.

4) 将前一步计算得到的正应力值代入等效连续介质受荷载作用下的渗透系数矩阵中.一般情况下，渗透系数的变化设为应力的指数函数，因此，等效连续介质受荷载作用下的渗透系数矩阵为

(4)

5) 重复1)～4)步的计算，直至计算结果达到程序中给出的相应允许值.

2工程概况

日照市沭水东调工程引水隧洞工程起点位于店子集镇工业园东北、何家庄村北；终点位于三庄镇西北、省道335北的三庄河左岸，隧洞全长约18.368 km，采用无压过流，洞身为城门洞形，毛洞断面尺寸为4.1 m×3.3 m加半圆拱，半圆拱直径4.1 m.整个隧洞分为入口段、洞身段以及出口段，其中大埋深段主要位于洞身段的东段，累计长度达6 km，最大埋深420 m左右，属于深埋小断面隧洞.隧洞洞身段沿线岩石分布大致可分为三段：西段岩性以安山质角砾石、长石砂岩和凝灰质长石砂岩为主；中段为白垩系与侵入岩群混合区域，岩性有长石砂岩、石英二长岩等；东段主要为侵入岩群，主要岩性为不等粒角闪石英二长岩、二长花岗岩和二长闪长岩等.岩体完整，节理裂隙较发育，裂隙面较平直光滑，连通性较好.以岩石强度、岩体完整性程度、结构面状态、地下水和主要结构面产状等五项因素总和的总评分为基本依据，围岩强度应力比为限定判据，该岩体完整段围岩分类为II类，较完整段围岩分类为III类，破碎带及构造断裂带影响带围岩分类为IV类，构造断裂带围岩分类为V类.

3计算模型

3.1模型建立

在数值模拟中，假定岩体为均质各向同性连续介质，模型范围：选取输水隧洞大埋深段，取其中100 m进行数值计算.该段隧洞穿越地层主要以花岗岩为主，地下水赋存、裂隙较发育，属弱-中等透水，围岩分类为II类.三维模型坐标原点位于模型中心剖面上输水隧洞横剖面顶拱的圆心处，该处实际高程为101.3 m，利用地形线生成曲面，用以模拟该段实际地表情况，确定该段最大埋深处实际高程为275 m.取隧洞轴线在水平面上的投影线为y轴，指向下游为正，x轴为水平面上垂直隧洞轴线方向，z轴以竖直向上为正.模型范围：-40 m≤x≤40 m，-50 m≤y≤50 m，z向原点以上取至地表，原点以下取30 m.采用MIDAS-GTS进行建模计算，模型为四面体剖分单元，单元数为55 938个，节点数为10 460 个.有限元分析本构关系采用弹塑性模型，计算采用Mohr-coulomb屈服准则，模型按全断面法开挖，每步开挖进尺5 m，共分20 个步骤开挖，采用锚喷支护，且支护滞后一个开挖施工步进行.MIDAS-GTS采用钝化单元来实现隧道的开挖，并采用激活单元实现支护的模拟.所谓钝化单元，即是此单元在该次计算中不发挥作用，相反，激活单元就是在计算中运用此单元.整个模拟过程：1) 初始地应力平衡，设置位移清零，得到初始应力场.2) 初始渗流分析，得到未开挖前的渗流场.3) 开挖一个施工步，即钝化相应的开挖土体单元.4) 开挖下一个施工步，并激活上一个施工步的支护单元，计算采用板单元模拟混凝土喷层，桁架单元模拟锚杆.应力计算的同时进行渗流计算，以此模拟渗流与应力的耦合作用.计算中不考虑开挖的时间效应.5) 重复以上步骤，直至施工结束.计算模型如图1所示.

图1　计算模型网格剖分图Fig.1　The net segmentation image of calculation model

3.2边界条件及参数确定

根据实际情况，并结合静力学中相关连续介质模型有关应力边界、位移边界和混合边界的规定，设定边界条件如下：底面限制z向位移，左右侧面限制x向位移，前后面限制y向位移，上边界为自由面.渗流边界定义中，模型的左右两侧设置竖直方向的节点水头高度为80 m，底部边界为不透水边界，假定隧洞开挖前节点水头以下的围岩饱和，以上的围岩为非饱和.具体力学参数见表1，2.

表1　岩体计算参数

表2　支护计算参数

4计算结果及分析

4.1围岩孔隙水压力分布特征

图2　隧洞周围各点孔隙水压力Fig.2　Pore pressure of feature point around tunnel

隧洞开挖改变了围岩的初始应力场,与此同时产生开挖松动区，并导致围岩损伤，节理裂隙进一步扩展.进而影响原岩内孔隙水压力的分布.利用MIDAS-GTS，对模型隧洞节点孔隙水压力进行监测，并将其开挖前后的孔隙水压力进行对比分析.对比分析结果见图2，未开挖时，地下水以静水压力的形式保持平衡，水平方向上孔隙水压力分布较均匀，且孔隙水压力随着深度增加而增加.开挖后，岩体内渗流场发生变化，隧洞周围较大范围内形成的孔隙水压力场重分布，隧洞掌子面及洞壁孔隙水压力为0 kPa.计算和监测均显示，靠近隧洞掌子面及洞壁的孔隙水压力梯度较大.图3为开挖后围岩孔隙水压力的分布图，其中等值线上孔隙水压力的数值单位为kPa，由图3可知:围岩孔隙水压力分布沿隧洞轮廓呈环形带状分布.

图3　开挖后岩体孔隙水压力等值线图Fig.3　The contour map of pore pressure after excavating

4.2围岩位移分布特征

图4(a)是不考虑孔隙水压力时的围岩竖向位移等值线图，图4(b)是考虑孔隙水压力时的围岩竖向位移等值线图，图5(a)是不考虑孔隙水压力时的围岩水平位移等值线图，图5(b)是考虑孔隙水压力时的围岩水平位移等值线图.如图4所示，未考虑渗流-应力耦合作用时Y=0断面拱顶的沉降值为2.18 mm，底部隆起值为2.97 mm.考虑耦合作用时，拱顶沉降值为2.97 mm，增大了0.79 mm，底部隆起值为2.35 mm，减小了0.62 mm.图5中，考虑耦合作用的水平位移比未考虑时要小，减小量约为1 mm.由图5可知：隧洞在不考虑孔隙水压力的影响下开挖支护完成后，隧洞底部的隆起值大于拱顶沉降值.当考虑渗流场的影响时，拱顶的沉降值增大，底部的隆起值与水平位移值均有一定程度的减小.产生这种现象的原因：隧洞开挖后，隧洞底部岩体因为卸荷作用而产生向上位移，但考虑渗流作用时，因孔隙水压力降低产生的固结沉降将使底板的隆起值小于不考虑渗流时的隆起值.这说明，隧洞的变形并不完全受应力影响，围岩的固结变形对隧洞变形同样有一定程度的影响.

图4　Y=0断面竖向位移等值线图Fig.4　The contour map of vertical displacement on plane Y=0

图5　Y=0断面水平位移等值线图Fig.5　The contour map of horizontal displacement on plane Y=0

为进一步确定渗流作用对围岩的影响，在数值计算过程中选取隧洞拱顶正上方不同的特征点进行沉降值监控，各监测特征点及对应坐标见表3.

表3　监测节点位置信息

通过计算，得出节点随开挖推进的竖向沉降值.从图6中可以看出：随着开挖的推进，围岩拱顶的竖向位移从最初的增大到逐渐趋于稳定，其中，节点845为距开挖区域较近的特征点，三个特征点中其沉降值最大，在考虑渗流-应力耦合作用时，其沉降值有所增大，但增幅较小.而节点493，该特征点距离开挖区域较远，三个特征点中其沉降值最小，但考虑渗流-应力耦合作用时，其增幅最大.由图6可知：拱顶上方距离开挖区域越远的点沉降值越小，受应力耦合作用影响越大.

图6　开挖过程中围岩竖向位移与开挖进尺关系曲线Fig.6　The relationship between excavation progress and vertical displacement

5结论

隧洞开挖破坏了岩体的初始应力场，降低了岩体的孔隙水压力，且越靠近开挖临空面，孔隙水压力梯度变化越大.开挖后的围岩孔隙水压力分布沿隧洞轮廓呈环形带状分布；在渗流-应力耦合作用下，围岩的变形空间分布规律与未考虑孔隙水压力的情况相似，隧洞拱顶沉降值有所增大，而隧洞围岩水平收敛值和底部隆起值均有不同程度的减小，孔隙水压力产生的固结沉降对围岩的位移有很大的影响；对比开挖后的孔隙水压力分布图与位移分布图，从中可知，开挖对孔隙水压力的影响范围大于开挖对位移的影响范围；从孔隙水压力对围岩沉降值的影响来看，隧洞拱顶沉降值受孔隙水压力影响最大，且拱顶上方距离开挖区域越远的点，沉降值越小，受渗流-应力耦合作用的影响越大.因此，在深埋小断面引水隧洞中，拱顶的位移稳定性较差，应加强防护.

参考文献：

[ 1]何川，谢红强.多场耦合分析在隧道工程中的应用[M].成都：西南交通大学出版社,2007.

[2]刘仲秋,章青.岩体中饱和渗流应力耦合模型研究进展[J].力学进展，2008,38(5):585-600.

[3]靳晓光,李晓红,张燕琼.越江隧道施工过程的渗流-应力耦合分析[J].水文地质工程地质,2010,37(1):62-67.

[4]刘仲秋,章青.考虑渗流-应力耦合效应的深埋引水隧洞衬砌损伤演化分析[J].岩石力学与工程学报,2012,31(10):2147-2153.

[5]李术才,宋曙光,李利平,等.海底隧道流固耦合模型试验系统的研制及应用[J].岩石力学与工程学报,2013,32(5):883-890.

[6]纪佑军,刘建军,程林松.考虑流-固耦合的隧道开挖数值模拟[J].岩土力学,2011,32(4):1229-1233.

[7]黄涛,杨立中.渗流与应力耦合环境下裂隙围岩隧道涌水量的预测研究[J].铁道学报,1999,21(6):75-80.

[8]王克忠,李仲奎.深埋长大引水隧洞三维物理模型渗透性试验研究[J].岩石力学与工程学报,2009,28(4):725-731.

[9]王克忠,王玉培,林峰,等.平行双隧道盾构法施工地表沉降仿真计算研究[J].浙江工业大学学报,2013,41(3):300-308.

[10]王克忠,程青云,王玉培,等.粉质砂性土地基中盾构区间隧道开挖过程数值计算研究[J].浙江工业大学学报,2012,40(5):587-590.

[11]张春阳,曹平,靳瑾,等.金川矿区深部巷道围岩流固耦合稳定性数值模拟[J].科技导报,2013,31(33):31-36.

[12]喻军,刘松玉.隧道洞口边坡变形控制与数值分析[J].浙江工业大学学报,2012,40(1):101-105.

[13]李新平,代翼飞,胡静.某岩溶隧道围岩稳定性及涌水量预测的流固耦合分析[J].山东大学学报,2009,39(4):1-6.

[14]魏纲,郭志威,魏新江,等.软土隧道盾构出洞灾害的渗流应力耦合分析[J].岩土力学,2010,31(s):383-387.

(责任编辑：陈石平)

Seepage-stress coupled analysis on the deformation stability of deep diversion tunnels

WANG Kezhong1, TANG Yuqiang1, LI Weiping2, KONG Lingmin3, QIN Shaokun1

(1. College of Civil Engineering and Architecture, Zhejiang University of Technology, Hangzhou 310014, China;2.Zhejiang Provincial Institute of Communications Planning,Design and Research ,Hangzhou 310014,China;3.Water Supply Management Office of Rizhao, Rizhao 276800, China)

Abstract:The stress redistribution of surrounding rock caused by groundwater seepage in the process of tunnel excavation has significant influence on the stability of surrounding rock. Based on the diversion tunnel project of water transferring in Rizhao, a numerical modeling analysis is conducted on the excavation process of diversion tunnels with and without seepage using the finite element analysis software MIDAS-GTS with the coupling theory of seepage and stress. The redistributions of pore water pressure and deformations in surrounding rock during the process of tunnel excavation are obtained. The calculation results show that the pore pressure around the tunnel decreases owing to excavation and the distribution shape is annular and zonal. The displacement distribution in surrounding rock is similar for coupling and uncoupling cases but in the former a larger influence is exerted on the vertical displacements of the vault.

Keywords:diversion tunnel; excavation; coupling of seepage and stress; pore pressure; displacement

收稿日期：2015-10-23

基金项目：清华大学水沙科学与水利水电工程国家重点实验室资助项目(sklhse-2014-C-02)；浙江省科技厅公益基金资助项目(2012C21105)

作者简介：王克忠(1965—) 男，山东冠县人，教授，研究方向为地下工程，E-mail：wkz@zjut.edu.cn.

中图分类号：TV672+.1

文献标志码：A

文章编号：1006-4303(2016)02-0207-05