基于横观各向同性模型的片岩隧洞稳定性研究

2019-11-26杨渝南高素芳谢晓康韩绍康

水资源与水工程学报 2019年5期

孙涛，杨渝南，高素芳，刘杰,2，谢晓康，韩绍康

(1.三峡大学土木与建筑学院，湖北宜昌 443002； 2.三峡大学三峡库区地质灾害教育部重点实验室，湖北宜昌443002)

1 研究背景

片岩在深部高地应力区域大量存在，可近似为层厚极小的层状岩体，其变形和强度具有明显的各向异性特性，这一特点又会对工程岩体稳定性产生特殊影响[1]。在室内试验方面，刘志勇等[2]利用MTS815岩石试验系统研究石英云母片岩在加、卸荷应力路径下的力学特性；周新等[3]探讨片岩加速蠕变阶段的轴向应变-应变率-偏应力关系特征；杨以荣等[4]系统研究高围压卸荷条件下石英云母片岩变形破坏特征及能量演化规律；何忠明等[5]通过室内试验和数值模拟分析层状岩体单轴压缩情况下的应力应变响应以及强度特征；于德海等[6]对干燥和饱水状态下的绿泥石片岩进行三轴压缩试验，获取水对试样强度和变形特性的影响规律；刘廷等[7]通过室内三轴试验，研究不同含水率状态下云母片岩的强度和变形特性；王章琼等[8]研究片岩在片理面倾角分别为0°、45°、90°时的变形参数各向异性特征及破坏机制；郑松林等[9]通过片岩在单轴应力作用下的岩石力学特性试验研究,得出岩石试块的强度特征、变形特征与作用应力的大小、作用应力与结构面之间夹角有密切联系的一般规律；张晓平等[10]研究了与片岩裂纹扩展相关的应力门槛值。

针对强度理论和本构关系方面，肖明砾等[11]基于带Kachanov 蠕变损伤因子的蠕变模型，研究丹巴水电站石英云母片岩的三轴蠕变特性及其各向异性特性；张春生等[12]以颗粒流方法为基础,开展不同尺寸、不同结构组成的各向异性特征研究,揭示变形参数和强度参数的变化规律；文献[13-17]针对层状岩体弹性参数的各向异性进行了试验与理论分析，并构建了相应的破坏准则，张学民等[18]研究了岩石材料的变形和强度各向异性性状、各向异性破坏准则和弹塑性本构关系及各向异性特征对层状围岩隧道变形和稳定性的影响；张嘉威[19]根据单轴压缩试验及统计损伤理论,建立了针对单轴压缩状态下的岩石损伤本构模型。

在工程应用方面，已有许多结合大型工程的片岩研究成果，如Nasseri等[20]研究了片岩力学参数和变形破坏特征的各向异性特征；刘杰等[21]、许振华[22]及单治刚等[23]运用大型岩土分析软件ADINA和单轴试验阐明了丹巴石英云母片岩各向异性特征，且岩体的变形和破坏规律主要取决于地应力及岩体的结构特征；在基于隧道等含片岩岩体的大型工程方面，Hwang等[24]研究了各向异性的片理角度对隧道稳定性的影响；Shrestha等[25]指出当隧洞通过高度片状化的岩层时，在塑性变形和剪切区域存在严重的隧道变形。

在上述各方面取得丰硕成果的同时，有以下几方面需进一步深入研究：

(1)不同片理角度下岩体力学参数取值无试验依据。很多学者在进行岩体参数的取值时往往按经验进行参数折减或仅按数值反演来确定，人为主观性太强，参数取值无试验数据支持。

(2)现有研究未系统考虑岩样在不同高度和不同片理角度下的岩样抗拉强度尺寸效应的规律。

(3)现有研究多数采用的计算模拟范围较大，在各边界施加的应力边界效应显著，应力影响范围较小，无法影响到内部局部精细区域，无法对局部区域进行精细的片理模拟和精细应力构造。

(4)现有研究未深入考虑围岩片理与隧洞洞轴不同夹角下隧洞围岩稳定性影响规律并给出最不利工况。

随着我国西部大开发和可持续发展战略的实施，西部区域的引水隧洞工程和交通隧洞工程正不断涌现，针对地下工程中围岩开挖后岩体各向异性显著的研究越来越成为热点。本文基于丹巴引水隧洞工程，通过对该工程中石英云母片岩的各项物理和力学性质进行试验研究，并建立隧洞模型，从位移、塑性区、塑性应变值3个方面对隧洞开挖后围岩稳定性进行评价，对于隧洞开挖加固设计有着重要的理论意义和工程参考价值。

2 工程概况

丹巴水电站位于四川省甘孜藏族自治州丹巴县境内的大渡河干流上，引水隧洞位于大渡河左岸山体中，全长约16 710 m，垂直埋深基本为500～900 m，最大埋深可达约1 220 m，大部分洞室位于石英云母片岩软岩区域，岩石坚硬程度初步分类属较软岩，岩体质量以Ⅳ类为主，成洞条件差，开挖后有可能出现较严重的大变形破坏现象，遇水也易发生软化现象，如图1所示。因此，丹巴水电站引水隧洞软岩成洞问题突出，为枢纽开发的关键技术问题之一。

图1 隧洞工程围岩破坏图

3 片理在不同空间角度下的参数测定

为获得石英云母片岩在高地应力下三轴加卸载力学特性和变形破坏规律，本文开展了室内三轴压缩试验测定不同围压作用下岩样的力学参数。

3.1 室内三轴试验

3.1.1 三轴试验方案试验以RMT-150C岩石力学试验系统为平台，三轴试验采用3组不同片理倾角，4种围压循环加卸载。由于引水隧洞段实测最大主应力为36 MPa，故将围压分为0～10、0～15、0～20、 0～30、0～35、0～40 MPa共6级，轴向和围压同时加载，速率均为0.2 kN/s，每次循环2～4次，轴向应力线性增加直至岩样破坏，不同倾角岩样三轴加卸载试验加载模式如图2所示，部分岩样三轴试验破坏形态图如下图3所示。

图2 不同片理倾角的三轴加卸荷试验

图3 部分岩样三轴试验破坏形态

3.1.2c、φ值各向异性的试验分析根据3组36个岩样的三轴试验成果，对石英云母片岩表现的特征绘制成图，如图4～5所示。

图4 干燥岩样弹性模量、黏聚力和内摩擦角与围压关系图

图5 饱和岩样弹性模量、黏聚力和内摩擦角与围压关系图

由图4～5可知，随着围压的增加，石英云母片岩的弹性模量、内摩擦角和黏聚力均与围压呈线性增大规律，但不同片理倾角岩样的各力学参数对围压的敏感性有所差异。

(1)在干燥条件下，围压对岩样力学参数影响分析：由图4分析可知，①当轴向力与片理面平行时，随着围压的增加岩样参数弹性模量和内摩擦角线性增大速率最快，斜率分别为0.2245 和0.2873；②当轴向力与片理面斜交时岩样黏聚力增大速率最快，斜率为0.0711；③根据斜率对比得出岩样力学参数对围压敏感性排序为内摩擦角>弹性模量>黏聚力。

(2)在饱和条件下，围压对岩样力学参数影响分析：由图5分析可知，①当轴向力与片理面平行时，随着围压的增加岩样弹性模量、内摩擦角增大速率最快，斜率分别为0.135和0.2907；②当轴向力与片理面斜交时，岩样黏聚力增大速率最快，斜率为0.0943；③根据斜率对比得出岩样力学参数对围压敏感性排序为内摩擦角>弹性模量>黏聚力。

(3)饱和石英云母片岩相对干燥状态的弹性模量降低20%～40%，黏聚力降低25%～40%，内摩擦角降低20%-30%。表明水对片岩的弹性模量、内摩擦角和黏聚力有显著影响，在实际工程设计应用中应当考虑水对岩体质量的劣化影响。

3.1.3 各向异性与围压关系分析由于片理的存在，石英云母片岩具有显著的各向异性特征，给出表征各向异性公式(1)和(2)。

(1)

(2)

式中：σ1为最大主应力，Pa；σ3为最小主应力(由于为单轴压缩，故σ3=0)，Pa；k1、k2分别为石英云母片岩的单轴抗压强度和变形模量各向异性参数；E∥，E⊥分别为平行、垂直节理的变形模量。k1和k2等于1时说明岩石各向同性，当其越偏离1时，各向异性越明显。

根据表征各向异性的公式(1)和(2)计算得出各向异性参数k1和k2，根据计算结果得出各向异性程度与围压关系图，如图6所示，在不同围压下对不同片理倾角的石英云母片岩的各向异性进行评价。

由图6可知，在低围压下k1>1且k2<1,各向异性显著，随着围压的增加，k2→1，而后逐渐减小；k1逐步降低，各向异性减小，而后逐步升高，各向异性增加。岩体强度受片理结构影响较大，在实际工程设计应用中应当充分考虑岩样的各向异性特性。

3.2 抗拉强度各向异性劈裂试验分析

3.2.1 试验方案该试验基于笔者自主研发的单个岩样多次劈裂测岩石抗拉强度的试验方法来测定不同片理产状的石英云母片岩抗拉强度。试验分为5组，每组4个岩样(岩样高度分别为2.5、3.5、4.5和5.5 cm)，具体试验方案如表1所示，通过对破坏后岩样进行重粘接，旋转岩样调整加载力与片理交角实现通过单个岩样测不同片理角度下岩石抗拉强度，两次劈裂后岩样如图7～8所示：

表1 劈裂试验方案表

3.2.2 不同高度岩样抗拉强度随片理角度的变化分析由图9可知，石英云母片岩的抗拉强度受控于其片理特征，相同高度下岩样的抗拉强度随片理角度线性增大；相同片理夹角下，岩样抗拉强度随岩样高度递增呈良好线性递增关系。徐燕飞等[26]研究结果显示，各种岩石的劈裂强度都存在尺寸效应，不同岩石劈裂强度尺寸效应主要受岩石的物质成分和结构所影响；硅质砂岩劈裂强度随试样厚度增加呈线性增大；本文研究试验结果与该结论观点相一致。

石英云母片岩的抗拉强度与其试验岩样高度和片理角度有关，为研究不同片理角度下岩样抗拉强度尺寸效应的规律，定义岩样抗拉强度尺寸效应(此处尺寸效应实际指岩样的高径比)系数为ksize=σ5.5/σ2.5，其中σ5.5是高度为5.5 cm岩样的抗拉强度，σ2.5是高度为2.5 cm岩样的抗拉强度，得出岩样抗拉强度尺寸效应系数ksize随片理角度的变化关系图，如图10所示。

图6各向异性程度与围压关系图图7第1次劈裂后岩样图8第2次劈裂后岩样

图9不同高度岩样抗拉强度随片理角度的变化关系图10抗拉强度尺寸效应系数随片理角度的变化关系

由图10可知，不同片理角度下岩样抗拉强度尺寸效应系数ksize均大于2，即相同角度下5.5 cm高度岩样测得的抗拉强度一般为2.5 cm高度岩样抗拉强度的2倍以上，可见岩样的抗拉强度具有显著的尺寸效应，且岩样抗拉强度尺寸效应与片理角度呈线性递减的关系。

3.3 横观各向同性本构模型研究

根据前期室内试验成果和现场围岩破坏特征的分析，石英云母片岩的变形与破坏均表现出显著的横观各向同性特征。

根据横观各向同性理论可知，在局部坐标系下，增量型弹性应力-应变的关系为：

[Δσ′]=[K′][Δε′e]

(3)

式中：σ′为岩体弹性应力；ε′e为岩体线应变, [K′]为局部刚度矩阵，即:

(4)

其中：

式中：E1和μ12分别为各向同性平面的弹性模量和泊松比；E3和μ13分别为垂直与同性平面的弹性模量、泊松比。

通过对不同片理角度的岩样进行室内岩石物理力学试验(图3、7、8)，得到平行片理和垂直片理方向岩样的力学参数，如表2，将其作为片理模型的计算参数，为公式(4)提供参数支撑。

表2 岩体计算参数表

4 计算参数及模型

4.1 计算参数

采用横观各向同性本构模型，在ADINA软件中通过空间两点坐标来定位材料主轴，通过改变空间两点坐标来调节材料轴的方位，如图11，材料主轴z方向取垂直于片理参数，x和y轴均以平行片理参数选取，在测点地应力实测与计算值(表3)的基础上得出空间任意角度的岩体片理模型的计算参数，实现各种片理类型及角度的模拟。

4.2 引水隧洞软岩段地质体大区域三维模型

沿着引水隧洞方向建立三维模型，如图12所示，模型顺引水洞方向取10 358.215 m，垂直引水洞方向取2 066.367 m，根据地形，最大相对高度3 750 m，选取实测地应力测点位置为ZK2和ZK3，ZK2和ZK3测点剖面位置如图13所示。

采用水压致裂法测取测点ZK2和ZK3的地应力实测值，通过调节所建立大区域三维模型边界应力的大小，反演分析得出与实测值相接近的反演计算值，如表3。

4.3 主应力构造精细应力场方法的提出

4.3.1 大区域模型反演分析存在的不足采用以上大区域模型反演分析可较好地得出引水隧洞大区域的整体地应力场，但大区域模型反演分析存在以下不足：①隧洞开挖各工况下围岩片理角度在不断改变，无法对局部区域进行精细的片理模拟和精细应力构造；②大区域模型只是在各边界面上施加应力，应力边界效应显著，应力影响范围较小，无法影响到内部局部精细区域。因此本文对大模型局部区域建立小模型，进行精细计算模型应力场构造。

表3 测点地应力实测与计算值 MPa

建立小模型后，若直接采用测得各正应力、剪应力在各边界上施加应力，各边界面均有3个不同的应力，共18个应力，构造的应力场空间受力复杂，无法对局部区域进行精细应力场精准构造，空间单元体的受力示意图如图14所示。

4.3.2 精确主应力构造精细应力场的方法本文基于此提出一种新的主应力构造精细应力场方法，提出的《一种精确主应力构造精细应力场的装置及方法》(201811096241.5)发明专利已公示[27]，且该方法在丹巴工程中得以成功运用。采用实际精确的主应力方式进行构造，首先根据ZK2测点实测地应力的6个应力分量σx、σy、σz、σxy、τyz、τzx可得出3个不变量I1，I2，I3。

I1=σx+σy+σz

(5)

I2=σxσy+σyσz+σzσx+2τxy+2τyz+2τzx

(6)

I3=σxσyσz+2τxyτyzτzx-2σxτyz-

2σyτzx-2σzτxy

(7)

由公式(5)～(7)推导可得出主应力方程为：

σ3-I1σ2-I2σ1-I3=0

(8)

运用Matlab软件，可求出三元一次方程组(8)，得到各主应力值大小，根据主应力值的大小，由公式(9)～(10)联立求解后按公式(11)计算主应力倾角αi和方位角θi，测点ZK2各主应力、方位角和倾角如表4。

(9)

(10)

公式(9)中li、mi、ni(i=1、2、3)分别为主应力σi对于X、Y、Z轴的方向余弦。

(11)

选取局部小区域，根据ZK2测点各主应力倾角和方位角建立隧洞模型，在模型边界上施加与之垂直的相应面力，其大小为3个主应力值，即为精细计算模型应力场构造，如图15所示。该方法用实际精确的3个主应力进行精细应力场精准构造，大大优化了空间多个应力构造应力场的复杂性并提高了精准度。

表4 测点ZK2各主应力、方位角和倾角计算值

根据表4，建立尺寸为240 m×240 m×240 m的立方体隧洞模型，隧洞断面为圆形，直径为16 m，片理层面产状取N65°、WNE∠80°，岩体强度参数和变形参数参照表2进行取值，对开挖后隧洞围岩稳定性进行评价。

5 隧洞开挖后围岩稳定性评价

5.1 空间多角度片理与隧洞轴向夹角关系说明

由于施工现场开挖各工况下围岩片理角度在不断改变，使得围岩片理与隧洞轴向夹角也在不断改变，施工现场实地考察结果显示，片理与洞轴方向水平面夹角的度数分布集中在45°～90°之间，局部地段也可能出现其他类型及各种角度穿切情况。由此设定了3类片理走向类型及0°、15°、45°、75°、90°片理进行模拟计算，如图16所示，假定洞轴方向为x方向，定义第1类45°片理面为与平面ACHF平行，第2类45°片理面为与平面ABHG平行，第3类45°片理面为与平面BCGF平行，其他角度亦按照此进行分类，以此对空间多角度片理与隧洞轴向夹角的关系进行说明。

5.2 空间多角度片理产状下隧洞围岩稳定性分析

本文以洞长36 m，洞径8 m的试验洞为例，研究洞室围岩的稳定性随片理产状的变化规律，从而预测最不利片理隧洞工况和隧洞失稳破坏可能性最大的位置。通过调整材料主轴来实现3类片理走向及0°、15°、45°、75°、90°片理进行模拟计算，下面主要从应力、位移、塑性区3个方面对开挖后隧洞围岩稳定性进行分析，计算云图中位移和应力数据均为以开挖隧洞为中心点的径向位移和径向应力，其方向以指向隧洞圆心为正，反之为负。

图11空间任意角度岩体片图12引水隧洞软岩段大区域图13ZK2和ZK3所在y向

理模型示意图三维模型图x=1671m处的位置图

图14空间单元体的受力示意图图15主应力构造精细应力场构造图图16片理的假设走向

5.2.1 应力分析隧洞开挖后洞室围岩处于受压状态，各工况下洞室周围应力分布基本一致，主要集中于洞室的左下方拱脚和右上方拱肩，水平方向应力均较小，没有明显的变化规律。

由计算结果得出不同片理工况下应力变化趋势图，如图17所示。

由图17分析可知，不同片理角度工况下围岩应力随着片理角度总体趋势是先递增而后递减。第2类和第3类不同片理产状的隧洞围岩的应力变化规律较为一致，均在片理角度为15°时，隧洞围岩应力值达到最大，分别为25.81和23.07 MPa，第1类片理变形规律略有差异，在片理角度为45°时，隧洞围岩应力值达到最大，为24.14 MPa，表明工程中3类片理走向类型在片理角度分别为15°、15°、45°时为最不利工况，围岩应力对隧洞稳定性影响较大，应重点加固并监测该工况下隧洞围岩的变形情况。

5.2.2 位移分析根据计算结果统计得出3类片理角度位移云图对比表，如表5所示。

由表5分析可知，第1类和第2类不同片理产状的隧洞围岩的变形规律较为一致，其最大变形多发生在左侧拱肩和右侧拱脚，第3类片理不同片理倾角下，变形规律略有差异，45°时最大变形位于拱顶和拱底。工程中应针对变形较大的区域重点加固监测。

表5 3类片理角度位移云图对比表

由计算结果得出不同片理工况下位移值随片理角度变化图，如图18所示。

由图18分析可知，3类片理走向类型均在片理角度为90°时，隧洞围岩位移值达到最大，分别为10.88、10.88、9.72 cm；表明工程中3类片理走向类型均在片理角度为90°时为最不利工况，隧洞围岩会出现较大变形或失稳，应着重考虑隧洞开挖的安全性。

5.2.3 塑性区分析根据计算结果统计得出3类片理角度塑性区云图对比图，如表6所示。

对表6各工况进行分析可知，开挖卸荷后3类各角度片理均出现较大塑性区，塑性区开展部位均类似，其主要集中于隧洞左侧拱底和右侧拱肩处。多数片理角度工况下，隧洞开挖后左侧拱底塑性区开展深度较右侧拱肩塑性区略大，仅在第1类45°和第3类45°工况下，右侧拱肩塑性区深度大于左侧拱底。

由计算结果得出不同片理工况下围岩塑性区深度和面积变化趋势图，如图19和20所示。

表6 3类片理角度塑性区云图对比表

图17不同片理工况下应力变化趋势图图18不同片理工况下位移变化趋势图

图19不同片理工况下围岩塑性区深度变化趋势图图20不同片理工况下围岩塑性区面积变化趋势图

由图19～20分析可知，不同片理角度工况下的围岩塑性区深度及面积差异很大，但总体趋势是随着片理角度的增加围岩塑性区深度和面积先递增而后递减。

(1)围岩塑性区深度及面积分析：①第1类片理角度隧洞开挖后，45 °时塑性区开展深度最大，为8.2 m，15°时塑性区面积最大，为139.7 m2；②第2类片理角度隧洞开挖后，45 °时塑性区开展深度最大，为8.3 m,且45°时塑性区面积也最大，为200.1 m2；③第3类片理角度隧洞开挖后，45 °时塑性区开展深度最大，为9.0 m，75°时塑性区面积最大,为188.5 m2。

(2)从围岩塑性区深度和面积分析可得出：①3类片理均在45°时塑性区深度最大，为最不利片理角度；②不同类片理塑性区面积无明显规律性，第1类在15°时为最不利工况，第2类在45°时为最不利工况，第3类则在75°时为最不利工况。