考虑互层状岩体接触状态的地下洞室围岩稳定分析

赵 勐, 肖 明, 陈俊涛 ,左双英

(1. 武汉大学 水资源与水电工程科学国家重点实验室, 湖北 武汉 430072; 2. 武汉大学 水工岩石力学教育部重点实验室, 湖北 武汉 430072; 3. 贵州大学 资源与环境工程学院, 贵州 贵阳 550025)

层状岩体在复杂地质构造作用下,由于具有多种优势节理面,其力学行为和破坏模式与一般岩体相异,表现出明显的各向异性[1],如琅琊山抽水蓄能电站、乌东德和白鹤滩水电站等.层状岩体的结构依其所赋存地质环境的多样性而具有不同形式,王思敬等[2]将其划分为3种主要类型:互层状结构、薄层状结构和夹层状结构.在层状岩体中进行地下工程建设,尤其对于软硬互层状岩体,由于其既具有层状岩体变形和强度性质的各向异性特征,又受到岩体软硬互层的影响,当岩层倾角较大时,结构面发育明显,在构造应力作用下,层间剪切、挤压破碎带较为常见.因此,针对互层状岩体,考虑不同岩性介质相互接触作用对岩土体发生剪切滑移破坏的影响,能更好地描述岩体破坏模式,具有重要的现实意义.

目前国内外学者关于层状岩体模拟主要有两类方法,一种是将层状岩体分离成层间岩块和层面,分别采用不同模型来模拟这两种介质,黄书岭等[3]分别对母岩和节理面建立不同的屈服准则,提出层状岩体多节理本构模型,描述其变形和强度各向异性;徐国文等[4]采用基于PFC-FLAC的连续-离散耦合算法,对岩石和层理面用不同本构模型模拟,并运用于层状软岩隧道围岩的稳定性分析中;另一种是将节理面等效考虑进岩体中,把层状岩体看成是横观各向同性体,利用层状弹性体系力学理论分析其变形和强度性质的各向异性特征.该方向的研究主要归结为以下两类:

1) 基于张量函数形式的强度准则.Pietruszczak等[5]引入微结构张量aij,将各向同性准则映射到各向异性中;钟世英等[6]将微结构张量参数引入到Jaeger针对岩体沿节理面滑移破坏提出的M-C准则中,并将其运用到白鹤滩水电站的柱状节理分析中;李良权等[7]提出一种改进的各向异性Hoek-Brown强度准则,并进行参数敏感性分析;但上述方法引入参数过多且确定麻烦.

2) 引入各向异性函数的强度准则.该法主要将层面对岩体破坏的影响因素及材料力学参数的各向异性考虑进各向同性的屈服准则中,推导岩体强度与相关参数之间的函数关系.Jaeger[8]认为抗剪强度是破坏面与主应力夹角的函数,并给出相关的表达式;Zhang等[9]在文献[8]的基础上提出具体的层状岩体黏聚力和内摩擦角的表达式;韩昌瑞等[10]建立相应强度参数方程,引入到M-C屈服准则中,并运用到深埋长隧道锚固支护中.

本文针对互层状岩体采用“隐含层面”的建模技术,在计算中通过输入层面方位角实现应力、位移和荷载在层状法向坐标系和整体坐标系之间转换[11].对于软、硬互层交界面处,由于层间剪切错动带发育,需考虑其应力及位移在层面两侧的不连续性[12],但常规连续有限元模型难以解决这个问题,本文引入FLAC3D专门模拟岩土体接触面力学行为的界面单元[13],分析计算中可能出现的多种接触状态并提出判别修正方法.针对接触面边界节点对发生的嵌入问题,基于Desai薄层单元的嵌入控制方法[14],改进界面单元边界节点对发生嵌入时的迭代算法.以巴基斯坦阿扎德帕坦水电站互层状岩体为例,应用本文的研究方法,分析在构造应力作用下,不同岩性介质相互接触作用对互层状岩体层间剪切挤压带和地下洞室的影响,以期对复杂互层状岩体中地下洞室围岩稳定分析进行有益探索.

1 层状岩体弹性本构关系

在地质构造作用下,层状各向异性岩体由于优势层理面的存在,具有明显的变形各向异性.一般将其等效处理成横观各向同性介质,在平行于层面方向是各向同性的,在垂直于层面方向具有与层面不同的弹性主方向.依据层状弹性体系力学理论[15],层状法向坐标系下层状各向异性介质的弹性矩阵如式(1)所示:

(1)

E1/2(1+μ1);G2=E1E2/[E1(1+2μ1)+E2];

由于局部的层状法向坐标系通常与计算的整体坐标系不一致,需要通过三维坐标系下的坐标转换矩阵R和应力转换矩阵S将层状法向坐标系Ox′y′z′下的刚度矩阵Κ′和弹性矩阵D′转换到整体坐标系Oxyz下的刚度矩阵Κ和弹性矩阵D[11]:

K=RTΚ′R;

(2)

D=STD′S.

(3)

2 互层状岩体接触系统数值模型

对于互层状岩体,在水平向构造应力和竖直向自重应力联合作用下,层间剪切、挤压破碎带发育,此种破坏特性在陡倾角岩层下更为常见,互层状岩体介质在外荷载作用下将难以保持整体连续变形.本文基于FLAC3D界面单元接触理论[13],推导出4种不同接触状态的判别准则及修正方法,对接触边界节点相互“嵌入”时的迭代方法进行改进,从而解决互层状岩体层间非整体连续变形问题.

2.1 接触面本构模型

本文的计算模型通过在互层状岩体层间接触部位设置FLAC3D中提供的界面单元的方法来模拟其接触特性.界面单元是无厚度接触面单元,附属于三维实体单元的网格面,由一系列三节点的三角形平面单元构成.因此,不像等效连续模型一样层间网格共用同一节点,接触面模型可以模拟互层状岩体层间滑移和分离行为[13].

未施加外荷载时,接触面处于黏结状态,其本构模型如图1所示,接触面的法向和切向接触力分别由下式确定:

(4)

(5)

图1 接触面单元的本构关系示意图

2.2 接触边界条件

外荷载施加前,互层状岩体接触面上接触节点处于“点对点”的接触方式,接触点对间存在黏结力.如图2所示,对于发生接触的两种介质M1,M2,相应的接触边界上的点对i和i′需要满足接触边界条件,包括位移接触条件和接触力边界条件.

对于接触面两侧的物质,在法向上假定其不会相互嵌入，即满足法向接触条件;在切向上当接触面接触时,若无相对滑动产生,则点对还应满足切向位移协调条件,如式(6)所示:

n·(ui′-ui)+d0≥0 ;

(6)

(7)

式中:n为接触节点对的法线方向矢量;d0为接触节点对i和i′的初始法向间隙;t代表与法线方向矢量n对应的切线方向矢量.

图2 接触示意图

发生接触时,接触面两侧节点对的法向接触力Rn1,Rn2和切向接触力Rs1,Rs2需满足以下条件:

(8)

式中，c,φ为接触面的黏聚力和内摩擦角.

2.3 接触面屈服准则

2.3.1 拉伸破坏

对垂直于接触面的拉伸破坏情况,先计算沿接触面法向局部坐标下的应力分量:

σ′=STσ.

(9)

根据单轴拉伸强度准则,由局部坐标下的应力分量校核接触面发生拉伸破坏时的应力状态,当

(10)

说明接触面的局部法向应力超过极限抗拉强度Rt,接触面节点脱离目标面,产生缝隙且切应力和正应力变为0.

2.3.2 剪切滑移破坏

对平行于接触面的破坏情况,根据Coulomb抗剪强度准则校核沿接触面是否发生滑移:

(11)

式中:f为接触面的摩擦系数;k为接触面滑动安全系数.此时沿接触面方向将发生剪切滑移破坏.

2.4 接触状态判别准则及修正方法

根据接触面边界上节点对的位移相互关系和受力状态,结合Coulomb抗剪强度准则和法向抗拉强度准则,分析节点对可能存在的4种接触状态:黏结、滑移、张开和嵌入状态.接触面上的节点荷载和该节点所处的接触状态有关,在弹塑性非线性分析时,节点对发生滑移时超出抗剪强度的切向接触力和发生张开时超出抗拉强度的法向接触力需进行转移[11].由于黏结接触无需对接触力进行修正，故对滑移、张开和嵌入这3种接触状态的判别准则及修正方法详述如下.

2.4.1 滑移接触状态

对于接触面边界节点对发生滑移时,需对下式进行判断:

(12)

(13)

若节点对同时满足式(12)和式(13)条件,则判定节点对发生滑移.此时节点对的受力情况如下式所示:

(14)

此时,沿接触面方向节点对将发生剪切滑移,应将超出抗剪强度的切向力

ΔRsi=Rsi-(cA+Rnitanφ),i=1,2

(15)

进行转移,将其转换成等效节点荷载,并将所有接触面单元的上述节点荷载累加，作为下一时步的增量荷载,同时该组节点对的法向接触力和切向接触力修正如下:

(16)

2.4.2 张开接触状态

对于接触面边界节点对发生张开时,需对式(12)中变量H值进行判断.

若H>d0,则判定节点对发生张开,此时节点对间不传力,法向刚度及切向刚度参数均设置为0,即此时步中该节点对的法向和切向接触力均修正为

(17)

直至该节点对重新闭合.需注意的是，对于发生开裂破坏的单元,在上述处理过后,直接标记分离状态,不再进行额外处理.

2.4.3 嵌入接触状态及调整

对上述两种接触状态进行判断及节点力修正时,还需进行节点对位移判别,以检查接触面两侧节点是否发生嵌入.对发生嵌入的节点对,基于Desai薄层单元的嵌入控制方法[14],对其位移和节点力进行修正,使相互侵入的节点退回到黏结或滑移接触状态,具体步骤如下:

1) 遍历每一接触面两侧节点对,计算式(12)中变量H值;

2) 若H≥0,则节点对未发生嵌入,继续遍历下一组节点对;

3) 若H<0,则节点对发生嵌入,应对节点对的位移进行修正.计算α值如下:

(18)

式中，λ为控制嵌入判别精度的百分数,通常取1%[14],当节点对嵌入值小于λd0时,忽略不计.

(19)

将求得的节点修正荷载ΔQ按等值反向的原则施加在接触面两侧的节点对上,更新节点接触力,重新计算节点变形,得到新的点对法向位移,依此循环迭代直至满足平衡条件.

2.5 接触判别计算流程

在进行互层状岩体层间接触计算时,先假定初始时步内接触节点对无相对滑动,均处于黏结接触状态.基于接触边界条件,求解当前时步下节点对法向、切向的接触力及位移,判别节点的接触状态并修正.在本文所提出的接触判别迭代计算中,需注意的是,对于节点处于张开状态,在进行接触力修正时,认为节点先退回至黏结接触,再达至其余接触状态[16].具体接触判别计算流程见图3所示.

3 工程实例及分析

3.1 工程概况

阿扎德帕坦水电站位于巴基斯坦的Jhelum河上,为该河段水电开发中的一级.引水发电系统位于河流左岸,电站采用地下厂房“一洞室”方案,仅保留主机洞室,主变电站及开关站均布置在户外,主厂房内安装4台机组,机组间距为40.0 m,主机洞室开挖尺寸长×宽×高为175.6 m×24.9 m×60.15 m.

阿扎德帕坦水电站工程区地质条件复杂,厂房区基岩岩性为砂岩与非砂岩类，呈互层状分布,且具有层状各向异性特征.地层主要为单斜构造,岩层总体走向NE10°,倾向NE,倾角65°～75°(地表较缓,厂房部位较陡)；受构造应力影响,层间剪切、挤压破碎带较为常见.

图3 接触判别计算流程

3.2 计算模型和计算条件

综合分析该水电站的地形、地质特征,考虑到互层状岩体发育明显,本文在计算时主要选取穿越2号、4号机组段的Ⅳ类岩层与周围的Ⅲ类岩层共同组成互层状岩体为研究对象,计算模型均采用八节点六面体单元进行离散,共剖分了266 104个单元和278 775个节点,其中开挖单元26 784个,如图4、图5所示.三维计算坐标x轴顺水流向与厂房纵轴线方向垂直,y轴与厂房纵轴线方向重合,z轴与大地坐标重合,3个方向的计算范围分别为170,160,236.80 m.根据设计院提供的资料,锚固支护型式及参数见表1,共分7期进行开挖计算分析,开挖分期模型见图5.

图4 互层状岩体三维有限元模型

图5 开挖支护模型

表1 洞室支护型式及参数

三维初始应力场根据设计院提供的实测地应力反演分析得到,侧压力系数取kx=1.3,ky=1.4,kz=1.1.互层状岩体砂岩与泥岩材料物理力学参数取值见表2.开挖过程中通过释放开挖单元应力并将其施加到开挖临空面周围单元来模拟开挖释放荷载,该地下厂房在开挖过程中采用喷锚支护等支护措施[17].锚杆和锚索的模拟采用FLAC3D中的Cable单元[13],能有效模拟锚杆对围岩的加固效果.

计算工况分两种:①采用横观各向同性模型分析互层状岩体对开挖地下洞室的影响;②采用本文改进的界面单元模拟互层状岩体层间相互接触作用对开挖地下洞室的影响.

表2 材料物理力学参数

3.3 计算结果分析

为了论证采用本文改进的界面单元模拟互层状岩体层间剪切、挤压破坏的合理性,本文分别采用横观各向同性模型和含界面单元模型进行数值模拟,选取主厂房洞室对比分析锚固支护分期开挖下洞周围岩的变形和应力分布规律.

3.3.1 围岩变形分布规律

对于工况①,采用不考虑接触的横观各向同性模型进行分析计算(见图6a),厂房开挖完毕,洞周围岩变形在3.00～9.00 cm,最大变形量出现在泥、砂岩层间过渡处10 m范围内,边墙位移从上至下先增大后减小.由于将层面等较大的节理等效考虑进岩体单元中,围岩层间相对位移为0,互层状岩体在岩层面两侧位移变化连续,没有出现位移值的不连续跳跃情况.

对于工况②,采用本文考虑多种接触状态的含界面单元模型(见图6b),厂房开挖完毕,洞周围岩变形在5.00～13.00 cm,相对于工况①增大约2.00～4.00 cm.从图6b可以看出,互层状岩体在层间出现了错动位移,取主厂房顶拱、边墙和底部不同监测点进行分析,主厂房洞室各监测点处错动位移随开挖分期的变化如图7所示. 随开挖进行,顶拱的位移逐步发生回弹,错动位移值略有减小,最终稳定在0.48 mm附近,主厂房洞室其余各部位的错动位移值逐渐增大.尤其对于上下游边墙侧,在开挖至第3～4期时,高边墙逐渐形成,错动位移值增加较为明显,开挖完毕时,该位移值约为4.91～5.25 mm,局部最大值出现在边墙侧中下部(见图7),该处高边墙效应明显,是围岩稳定分析中重点支护部位.对比两种工况,结果表明对于互层状岩体,当考虑层间接触作用后,围岩变形要比等效连续模型的略大,且导致厂房边墙侧产生的相对变形要比顶拱和底部大,层间接触面处于滑移状态.

图6 主厂房洞周位移分布

图7 主厂房洞室围岩各监测点相对位移变化过程

3.3.2 围岩应力分布规律

采用不考虑接触的横观各向同性模型进行模拟时，主厂房洞室开挖完毕，围岩边墙侧第一主应力不断减小，在边墙一定范围内形成了应力松弛区，量终第一主应力稳定分布在-4.65～-0.96 MPa，第三主应力分布在-1.23～0.20 MPa.由图8a可知，互层状岩体在层间剪切带附近主应力值有所降低，降幅在1.0 MPa左右，该附近应力云图沿接触面向下突出，但是仍然保持连续变化，未产生突变.

工况②采用本文考虑多种接触状态的含界面单元模型进行计算分析,开挖完毕后,主厂房洞室边墙侧第一主应力分布在-6.68～-1.05 MPa,第三主应力分布在-2.03～0.41 MPa,相较于工况①的计算结果,应力增加约为10.53%～43.66%,上下游边墙侧局部拉应力超过岩体抗拉强度,导致部分围岩开裂破坏,且主要出现在层间剪切带附近(见图9b).表明对于互层状岩体,泥、砂岩软硬间错分布,软岩穿过部分厂房边墙处应力较大,受力情况复杂,层状岩体层间相对滑移对厂房开挖时应力扰动影响较大.

图8 主厂房洞周围岩第一主应力分布

图9 主厂房洞周围岩第三主应力分布

工况②采用改进的界面单元接触算法,在层间接触带两侧应力的变化不连续,软岩与两侧硬岩的应力值相差较为明显,呈现交错分布现象[18].从图8b和图9b可看出,在互层状岩体层间过渡处,第一主应力值出现明显跌落,应力降幅在1.4 MPa左右,第三主应力降幅在0.6 MPa左右.

从图8和图9可以看出,采用改进的接触面单元模拟互层状岩体层间的接触状态,洞周围岩位移和应力呈现非连续变化,能够比较好地模拟互层状岩体层间剪切、挤压破坏情况,与工程实际更为吻合.

4 结 论

1) 基于FLAC3D界面单元接触理论,改进界面接触面模型算法,针对界面单元接触边界发生黏结、滑移、张开和嵌入等接触状态建立了相应的判别准则和修正迭代算法,充分考虑互层状岩体层间相互作用特点,能有效模拟构造应力作用下在互层状岩体中开挖地下洞室时发生的非线性剪切滑移破坏现象.

2) 将本文所建模型应用于地下洞室的开挖计算中,与传统有限元的横观各向同性模型对比,结果表明,当考虑多种接触状态时,岩体层间过渡处应力和位移呈现非连续变化规律,且量值也进一步增大.在构造应力和开挖荷载作用下,互层状岩体层间易发生相互错动,且在上下游边墙侧中下部较为明显.本文模型可运用于实际工程计算，为互层状岩体层间非连续性接触问题的模拟提供了一种有效的实现途径.