石立国 郑 雄 祝 羿

中国建筑第二工程局有限公司 北京 100054

岩体稳定性研究一直是岩土工程领域的基本问题之一，尤其是在众多坡体失稳等自然灾害发生之后，更得到了广泛的关注，目前岩体变形与稳定性已成为岩土工程中重要的研究方向。

为保证工程安全性，通过监测分析并采取相应的预报控制措施，有利于施工过程中的风险控制。

岩土工程稳定性分析常使用的有极限平衡法和数值计算方法[1-2]，对其变形失稳机制的研究，目前多采用野外监测数据分析、物理模型试验以及数值模拟分析等方法[3-4]。有限差分数值计算方法，可适用于比较复杂的边界条件，充分考虑了岩土体的非物质性和不连续性，可以直观地模拟岩体外部形态和内部介质的变化情况。实体内部介质通过有限差分法，能直观地显示实体位移量、应力特征和潜在破裂区等[5-6]，对其稳定性状况以及易破坏风险区进行预先评判，从而可降低施工中发生事故的概率。

本文将结合数值模拟技术与现场监测方法对隧洞口岩体边坡的变形特征与安全状况展开分析和对比，保证施工过程中隧洞与周边环境的安全。

1 工程地质条件

某水电工程地下电站尾水隧洞出口单级坡比为1∶0.3，岩体主要由灰色薄层灰岩、灰色中厚层灰岩地层构成。岩体边坡结构为陡倾横向坡，强卸荷带水平宽一般为10 m，弱卸荷带一般为16～25 m。岩体边坡中部高程950～1 060 m间的缓坡对上部的崩塌块石有缓冲作用，边坡上发育的块体可能影响施工及工程运行安全。右岸尾水出口处岩体岩性复杂，影响范围内的岩体主要由落雪组第八、九、十段组成。整体围岩以Ⅱ、Ⅲ类为主，岩体力学参数较好，边坡表层有浅层强、弱卸荷带，岩体材料参数相对较差。

2 变形监测分析

隧洞出口岩体坡表面布置有多个监测点，选取6个典型代表性的有效测点TP01～TP06，高程依次分别为962、895、866、956、926、867 m。

由6个表面位移监测点的水平（朝向临空方向）变形累计值（图1）可知，近5年来，各部位水平位移均变化不大，最大不超过9 mm，总体的变化规律趋势具有相似一致性，协调性较好。其中，水平位移最大部位为867 m高程的测点TP06，但位于同样低高程的测点TP03位移却相对较小。900 m高程以上的3个测点位移值均较小，变化规律相似，最大不超过5 mm。由此可见，高程高的部位变形相对较小，但水平位移变化与高程并无必然联系。

图1 水平变形累计值

由表面位移监测点的垂直位移累计值（图2）可知，近5年来，各部位的垂直位移变化总体较平稳，在平稳中逐渐增大，变化规律趋势具有相似一致性。其中，垂直位移最大部位为866 m高程的测点TP03，最大沉降达15 mm，位于几乎相同高程的测点TP06沉降也大于其他高程的测点。900 m高程以上的3个测点垂直位移值均较小。

图2 垂直变形累计值

由此可见，高程低的部位沉降普遍大于高程高的部位。

由表面位移监测点的水平位移月动态改变量（图3）可知，近5年来，各部位的水平位移处于动态变化之中，处于0值上下变化，基本在－4～4 mm间起伏波动，总体变化平稳，动态变化规律也具有一定的相似性。高程低的点TP06变化相对最大，但月变化量不超过5 mm。

图3 水平位移月动态改变量

从位移监测数据分析可知，目前位移量值较小，总体变化平稳，无突变、陡增现象，岩体边坡暂处于稳定状态。

3 稳定性计算

3.1 计算方法与本构模型

采用Flac3D中Fish语言自编强度折减程序，对稳定性进行有效计算，其中的参数强度折减为：

式中：F——强度参数折减系数；

c、φ——有效黏聚力和内摩擦角；

c'、φ' ——折减后的有效黏聚力和内摩擦角。

本构模型为各向同性弹塑性模型，采用Mohr-Coulomb（M-C）准则，其力学模型为：

式中：σ1、σ3——最大、最小主应力；

f——屈服函数，且f＞0时，材料处于塑性流动状态，f＜0时，材料处于弹性变形阶段，f=0时，处于弹、塑性的临界状态。

以上即为剪切破坏判据，拉伸破坏判据为：

式中：σt——岩体抗拉强度。

在有限差分程序中，岩体的体积模量、剪切模量分别由下式计算：

式中：K、G——体积模量、剪切模量；

E——弹性模量；

ν——泊松比。

3.2 数值计算模型建立

根据该区工程地质条件，建立区域三维数值计算模型（图4）。

计算模型主要包括尾水出口边坡、尾水隧洞、高导流洞等，坐标原点位于高导流洞出口中心点处，x轴与尾水洞轴线重合，指向尾水洞出口方向为正，y轴从6#尾水隧洞指向4#尾水隧洞为正；z轴与大地坐标系重合，向上为正。模型计算范围为：268.18 m×230.90 m×339.00 m（x×y×z）。共划分单元47 665个、节点51 202个。数值计算共分11期开挖，前2期开挖为高导流洞、尾水隧洞间隔开挖，第3～11期为出口边坡开挖。数值模型各材料参数值如表1所示。

图4 数值计算模型

表1 岩体力学参数

3.3 稳定性计算条件

初始地应力场采用自重应力场，其中侧压力系数根据附近地应力测点实测结果综合比较，选取0.8。采用三维弹塑性有限差分法对整体模型进行稳定性计算分析。计算工况采用以下2种：

1）施工工况：模拟出口边坡的施工开挖与支护过程，该工况分2种计算方案，分别为无支护和施加系统支护。

2）正常运行工况：分析水库蓄水后人工边坡的变形、应力与支护受力情况（水库正常蓄水位为975 m），下游尾水位为847.86 m。

3.4 施工工况结果分析

破坏区体积总量及破坏区深度如表2、表3所示。由隧洞口破坏区分布（图5）可知，无支护和有支护方案，边坡整体破坏区分布规律基本相同，边坡岩体破坏区主要位于尾水出口洞脸部位。有支护方案边坡破坏区范围、深度、破坏总体积较无支护方案有明显减小，其中开裂区深度最大减小55.0%，塑性区深度最大减小63.4%，破坏总体积减小46.2%。边坡破坏区主要集中在洞脸部位，最大破坏区深度达9.56 m，位于5#尾水隧洞附近。5#、6#尾水隧洞出口岩柱破坏区贯穿，深度约7 m，在系统锚杆控制范围内。

由施工期尾水出口边坡主应力分布计算结果（图6、图7）可以看出，相对于无支护条件，锚固支护在一定程度上改善了边坡岩体的受力情况，第1、3主应力均有一定程度减小，且整个坡体应力分布更为均匀。有支护条件下，边坡开挖完成后，与初始地应力相比，受开挖卸荷影响在表层强卸荷带和浅层弱卸荷带坡体第1主应力有所减小，量值范围为－1.18～0.08 MPa，主要为压应力，第3主应力范围为0.06～0.99 MPa。下部新鲜岩体竖直边坡段第1主应力较初始地应力有所增大，量值范围为－3.68～0.08 MPa，第3主应力主要为拉应力，量值范围为0.06～1.30 MPa，局部出现一定程度的应力集中现象。总体看来，边坡应力远小于该部位岩体的抗压、抗拉强度，因此局部应力集中不会引起边坡的整体失稳。

表2 施工期不同计算方案边坡破坏区体积

表3 2种支护条件下开挖完毕典型断面破坏区深度（单位：m）

图5 洞脸边坡破坏区分布

由岩体位移分布（图8）可知，不同的计算方案，边坡岩体的位移分布规律大致相同，从边坡顶部至底部位移逐渐增大。位移较大处均位于落雪组第九段直立边坡处。有支护较无支护方案，边坡变形有明显减小，其中强卸荷部位坡体变形范围为0.5～1.7 mm，相应减小约60.5%，弱卸荷部位坡体变形范围为0.2～5.1 mm，相应减小约50.1%，新鲜岩体斜坡段变形范围为1.2～7.5 mm，相应减小约47.6%，新鲜岩体竖直坡段变形范围为2.1～12.1 mm，相应减小约34.9%。可见，喷锚支护对边坡变形有较好的限制作用，整体边坡支护参数较为合理。

图7 施工期第3主应力分布云图

图8 施工期位移分布云图

3.5 正常运行工况结果分析

正常运行期，区域岩体破坏区分布规律与有支护开挖完毕时基本相同。相对于无支护情况，直立边坡塑性区大部回弹，且大量开裂区调整为塑性区，开裂区深度减小达50%。总破坏体积比有支护条件减小约19.8%。

与施工期相比，边坡上部应力基本相当，下部第1主应力范围为－3.25～0.32 MPa，第3主应力范围为0.06～1.06 MPa。最大压应力为－3.92 MPa，最大拉应力为1.20 MPa，均出现在隧洞出口附近，边坡第1主应力压应力值不大，拉应力值也远低于边坡岩体的极限抗拉强度。

运行期边坡上部斜坡段位移在8.0 mm以内，最大位移为11.1 mm，出现在5#、6#隧洞洞脸处。

运行期边坡锚杆应力范围为0.1～60.5 MPa。相对于施工期，正常运行工况锚杆应力略有减小，最大锚杆应力值为60.5 MPa，比施工期减小约2.3 MPa，最大锚杆应力出现在5#、6#隧洞。边坡锚索应力范围为1 016.9～1 025.4 MPa。总体而言，锚杆、锚索应力水平较低，均在设计强度范围内。综上所述，可以认为，该导流隧洞出口边坡开挖、支护设计合理，岩体边坡整体稳定性较好。

4 结语

1）通过多年实测变形数据分析，隧洞出口边坡岩体变位平稳，总体正常，没有出现位移突变的迹象。

2）通过对比数值分析和实测数据结果，模拟反馈结果与监测数据具有良好的一致性，验证了利用有限差分法开展岩体变形与安全分析的合理性和可行性。

3）通过有限差分数值计算法对隧洞出口岩体边坡稳定性进行数值模拟和分析评价，预测潜在易破坏区，具有较好的效果，且有限差分数值计算方法考虑了岩土体的非线性大变形，克服了传统极限平衡分析和有限元法的不足。