地下洞室开发需求日益增加，近几年正在进行设计或施工的地下工程数量急剧上升。随着越来越多地下洞室埋深越来越大，施工时面临的问题也越来越多，施工难度也随之增大。例如在高地应力地区，随着开挖进行，洞室围岩经常会出现纵向的劈裂裂纹，对洞室的可靠性造成极大的影响^[1]。Hibino等^[2]利用多种方式在日本多个地下洞室进行监测后发现，大多数埋深较大的地下洞室边墙围岩内部都存在着由劈裂产生的裂隙张开或扩展现象。在中国，研究者也监测到了大量的劈裂裂纹出现在水电站工程地下洞群边墙等部位。这些裂纹的存在极大地影响着地下洞室群的稳定性和可靠性^[3–4]，还可能影响吊车梁的承载能力^[5]及锚杆锚索的可靠性^[6]，甚至可能导致厂房出现渗水^[7]等。

到目前为止，众多学者对裂纹的扩展及围岩的稳定性等相关内容进行了深入的计算研究；基于传统方法提出很多行之有效的本构模型及计算方法^[8–10]，并以此推动了岩石力学的发展。但近年来，地下工程开挖的复杂程度越来越大，所面临的工程难度和问题也越来越尖锐，传统的本构模型不能很好地描述地下洞室围岩劈裂破坏区的力学性质；一些基于传统的模型和方法在计算和处理岩体劈裂问题的时候也存在很多不足与局限。特别是在高地应力地区的地下洞室，当洞室围岩的细小裂隙扩展成为平行裂隙组之后，很多传统研究裂隙的方法或手段就已不再适用，而合适的研究方法和力学模型极大地影响着研究结果的正确与否。因此，还需寻找更加适合描述高地应力环境下洞室围岩劈裂破坏的力学模型及研究方法。

作者拟结合大岗山洞群工程，以现场监测资料为依据深入分析，运用近年新提出的劈裂破坏准则，结合劈裂岩体的力学特性与能量耗散理论建立更符合现场情况的力学模型进行计算，进而推算出劈裂区产生的部位、深度及随着开挖过程的发展规律，并指导类似工程。

1 能量耗散模型的建立 1.1 加卸载条件下的变形模量

通常情况下，脆性岩体受到荷载时，随着作用在岩体上的应力增加，岩体产生的应变也增大。但当应力到达岩体峰值强度时，应变若继续变大，岩体能承受的压应力将会很快跌落到一个较低状态。这是由于当岩体所受到的应力水平到达一定程度后，就会产生许多微裂隙。与此同时，将会通过声能、热能和裂隙表面能等方式损耗许多能量，而这些能量的耗散多是不可恢复的。因此，当岩体到达这一阶段后，即使卸载，应变也不会恢复到初始状态，产生残余变形。这种现象体现在应力–应变曲线上可描述为：岩体卸载阶段的曲线并不是沿加载阶段曲线的路线返回，而是从加载线下部返回，从而与加载阶段曲线构成环状，即为滞回环，如图1所示。

整个加卸载过程中流失的能量总值即为滞回环所围成的面积。为简化问题，将其加卸荷曲线假设为上行和下行的直线，以便于进行近似分析计算，见图2。

处于加载阶段时，试件的弹性模量为 $E$ 。当加载到 $A$ 时，试件包含的总能量为 ${S\!_{AOC}}$ 。开始卸载时，由于部分能量耗散，试件存在残余变形，因此假设卸载时的弹性模量为 $\overline E $ 。当应变恢复至 $B$ 点时，此过程中试件恢复的能量为 ${S\!_{ABC}}$ ，耗散的能量为 ${S\!_{AOB}}$ 。由图2可知：

$\frac{{\overline E }}{E} = \frac{{{{OC}}}}{{{{BC}}}}$

(1)

根据能量守恒可得：

$\frac{{{{OC}}}}{{{{BC}}}} = \frac{{{S\!_{{{AOC}}}}}}{{{S\!_{{{ABC}}}}}}$

(2)

所以：

$\overline E = \frac{{{S\!_{{{AOC}}}}}}{{{S\!_{{{ABC}}}}}}E = \left( {1 + \frac{{{S\!_{{{AOB}}}}}}{{{S\!_{{{AOC}}}} - {S\!_{{{AOB}}}}}}E} \right)$

(3)

引入等效系数 $t$ ，令 $\left( {1 + \displaystyle\frac{{{S\!_{{{AOB}}}}}}{{{S\!_{{{AOC}}}} - {S\!_{{{AOB}}}}}}E} \right) = t$ ，则 $\overline E = tE$ 。

因此，考虑能量耗散作用时有：

${\text{弹性模量}} = \left\{ \begin{aligned}& {\text{加载阶段}},E;\\& {\text{卸载阶段}},\overline E = tE\end{aligned} \right.$

(4)

由上述理论可知，引入能量耗散的实质即为脆性模型在不同阶段使用的弹性模量不同，即在加载时期使用其本身的弹性模量 $E$ ，而使用等效弹性模量 $\overline E $ 进行卸载^[11–12]。

通过以上分析，在使用FLAC3D进行实际计算时，模型的加载阶段可以使用原本构方程不变；在卸载阶段时，由于使用了等效弹性模量 $\overline E $ ，则卸载阶段的本构方程为：

$\left\{ \begin{aligned}& \Delta {{{\sigma }}_{11}} = \overline {{\alpha _1}} \Delta {\varepsilon _{11}} + \overline {{\alpha _2}} (\Delta {\varepsilon _{22}} + \Delta {\varepsilon _{33}}),\\& \Delta {{{\sigma }}_{22}} = \overline {{\alpha _1}} \Delta {\varepsilon _{22}} + \overline {{\alpha _2}} (\Delta {\varepsilon _{11}} + \Delta {\varepsilon _{33}}),\\& \Delta {{{\sigma }}_{33}} = \overline {{\alpha _1}} \Delta {\varepsilon _{33}} + \overline {{\alpha _2}} (\Delta {\varepsilon _{11}} + \Delta {\varepsilon _{22}}),\\& \Delta {{{\sigma }}_{12}} = 2\overline G \Delta {\varepsilon _{12}},\\& \Delta {{{\sigma }}_{13}} = 2\overline G \Delta {\varepsilon _{13}},\\& \Delta {{{\sigma }}_{23}} = 2\overline G \Delta {\varepsilon _{23}}\end{aligned} \right.$

(5)

式中： ${\varepsilon _{11}}$ 、 ${\varepsilon _{22}}$ 、 ${\varepsilon _{33}}$ 分别为在 ${{{\sigma }}_{{11}}}$ 、 ${{{\sigma }}_{{22}}}$ 、 ${{{\sigma }}_{{33}}}$ 主应力空间中对应的主应变分量； ${\varepsilon _{12}}$ 、 ${\varepsilon _{13}}$ 、 ${\varepsilon _{23}}$ 分别对应 ${{{\sigma }}_{{12}}}$ 、 ${{{\sigma }}_{{13}}}$ 、 ${{{\sigma }}_{{23}}}$ 的工程应变分量； $\overline {{\alpha _1}} {{ = }}\overline K {{ + }}\displaystyle\frac{{{4}}}{{{3}}}\overline G $ ； $\overline {{\alpha _2}} {{ = }}\overline K {{ - }}\displaystyle\frac{{{2}}}{{{3}}}\overline G $ ； $\overline K $ 、 $\overline G $ 分别为卸载阶段的体积模量和剪切模量，两者的数值根据现场的岩石情况进行确定。

1.2 加卸载判据

一般来说，基于加卸载概念的加卸载判据主要有两种^[13]。其中，可以使用应力张量第一不变量 ${{{I}}_{1}}$ 和应力偏量第二不变量 ${J_2}$ 的变化作为加卸载判据。根据定义， ${{{I}}_{1}}$ 和 ${J_2}$ 增大时可认为试件为加载状态，而 ${{{I}}_{1}}$ 和 ${J_2}$ 减小时则认为是进行卸载。加卸载判据为：

$\left\{ \begin{aligned}& {{{I}}_1} = \frac{1}{3}\left( {{{{\sigma }}_1} + {{{\sigma }}_2} + {{{\sigma }}_3}} \right),\\& {J_2} = \frac{1}{{\sqrt 2 }}{[{\left( {{{{\sigma }}_1} - {{{\sigma }}_2}} \right)^2} + {\left( {{{{\sigma }}_2} - {{{\sigma }}_3}} \right)^2} + {({{{\sigma }}_3} - {{{\sigma }}_1})^2}]^{\textstyle\frac{1}{2}}}\end{aligned} \right.$

(6)

但在本文假设下，加载和卸载之间没有明显的屈服条件，因此体积模量和剪切模量均为应力和应变张量不变量的函数^[11]。为计算简便，仅通过应力偏量第二不变量 ${J_2}$ 的状态判断模型为加载或卸载状态，即 ${J_2}$ 增加即为加载， ${J_2}$ 减小即为卸载^[14]。这样既简化了计算，又可避免出现 ${{{I}}_{1}}$ 和 ${J_2}$ 不同时增大或减小时的加卸载状态判断问题，同时还能很好地反映 3 个主应力对加卸载的影响^[15]。

根据上述理论，利用FLACE3D的二次开发功能，在C++的编译环境下对模型进行改进：在原有模型中导入加卸载判据；在判定加载阶段时，依旧使用初始的弹性模量 $E$ ，即保留原有的本构模型；在判定为卸载阶段时，使用卸载时的等效弹模 $\overline E $ 代替原弹性模量。最后导出模型，进行 FLAC3D的循环求解计算。

2 横观各向同性弹性体的本构关系

地下厂房开挖时，由于岩体的脆性特征，洞室之间的岩体常常会出现劈裂性的竖向裂缝。通常可以将这种劈裂破坏的发展描述为以下3个阶段：当地下洞室开挖后，开挖区域产生的卸荷将会使洞室围岩出现类似受压的现象，这种现象将会产生洞室围岩的应力集中^[16]，如图3（a）所示。当压应力继续作用于岩体上时，原有裂隙裂纹将会出现扩展，如图3（b）所示。当应力值到达岩体开裂条件时，围岩中之前产生的裂纹会沿最大压应力方向继续扩展；如果压应力继续加载在岩体上，此时由于自由边界的存在，裂纹不再稳定的扩展，而加速扩展至形成尺寸较大的贯通的裂隙组，如图3（c）所示。这种破坏形式下，完整的围岩被竖向裂纹分成了层状围岩。所以在劈裂破坏区，可将岩体简化为层理面垂直于水平面的层状模型，将劈裂破坏区的岩体用横观各向同性模型描述。

横观各向同性弹性体由一系列弹性对称面组成^[17]，对称面将弹性体分割开来，形成了两种不同的力学参数，其中与对称面平行方向上为一种参数，垂直于对称面方向上为一种参数，如图4所示。

若将 $XOY$ 平面假设为弹性对称面，那么横观各向同性体弹性参数应满足如下条件：

$\left\{ \begin{aligned}& {E_x} = {E_y} = E,{E_{\textit{z}}} = {E'};\\& {\mu _{xy}} = {\mu _{yx}} = \mu ,{\mu _{{\textit{z}}x}} = {\mu _{{\textit{z}}y}} = {\mu '};\\& {G_{xy}} = \frac{E}{{2\left( {1 + \mu } \right)}},{G_{y{\textit{z}}}} = {G_{x{\textit{z}}}} = {G'}\end{aligned} \right.$

(7)

横观各向同性弹性体的本构关系可简化为：

$\left[\!\!\!\!\! {\begin{array}{*{20}{c}}{{{{\varepsilon }}_{{x}}}}\\{{{{\varepsilon }}_{{y}}}}\\{{{{\varepsilon }}_{{z}}}}\\{{{{\gamma }}_{{{xy}}}}}\\{{{{\gamma }}_{{{{y}}{{z}}}}}}\\{{{{\gamma }}_{{{{x}}{{z}}}}}}\end{array}} \!\!\!\!\!\right] = {{A}}\left[\!\!\!\!\! {\begin{array}{*{20}{c}}{{{{\sigma }}_{{x}}}}\\{{{{\sigma }}_{{y}}}}\\{{{{\sigma }}_{{z}}}}\\{{{{\tau }}_{{{xy}}}}}\\{{{{\tau }}_{{{{y}}{{z}}}}}}\\{{{{\tau }}_{{{{x}}{{z}}}}}}\end{array}}\!\!\!\!\! \right]= \left[\!\!\!\!\! {\begin{array}{*{20}{c}}{{a_{11}}} & {{a_{12}}} & {{a_{13}}} & 0 & 0 & 0\\{{a_{12}}} & {{a_{11}}} & {{a_{23}}} & 0 & 0 & 0\\{{a_{13}}} & {{a_{23}}} & {{a_{33}}} & 0 & 0 & 0\\{{0}} & {0} & {0} & {{a_{44}}} & 0 & 0\\{0} & {0} & {0} & {0} & {{a_{55}}} & 0\\{0} & {0} & {0} & {0} & {0} & {{a_{66}}}\end{array}} \!\!\!\!\!\right]\left[\!\!\!\!\! {\begin{array}{*{20}{c}}{{{{\sigma }}_{{x}}}}\\{{{{\sigma }}_{{y}}}}\\{{{{\sigma }}_{{z}}}}\\{{{{\tau }}_{{{xy}}}}}\\{{{{\tau }}_{{{{y}}{{z}}}}}}\\{{{{\tau }}_{{{{x}}{{z}}}}}}\end{array}} \!\!\!\!\!\right]$

(8)

式（8）中柔度矩阵 ${{A}}$ 所含参数可用工程弹性常数表示为：

$\left\{ {\begin{aligned}& {{a_{11}} ={a_{22}}= \frac{1}{E},{a_{12}} = - \frac{\mu }{E},{a_{13}} = - \frac{{{\mu '}}}{E}}; \\& {{a_{33}} = \frac{1}{{{E '}}},{a_{44}} = \frac{{2\left( {1 + \mu } \right)}}{E},{a_{55}} ={a_{66}}= \frac{1}{{{G '}}}}\end{aligned}} \right.$

(9)

式中， $E$ 、 $\mu $ 为横观各向同性面（ $XOY$ 平面）内的弹性模量和泊松比， $E'$ 、 $\mu '$ 为垂直横观各向同性面（ ${\textit{z}}$ 轴方向）内的弹性模量和泊松比， $G'$ 为垂直横观各向同性面（与 $XOY$ 平面垂直）内的剪切模量。

因为刚度矩阵与柔度矩阵互为逆矩阵，即 ${{C}} = {{{A}}^{ - 1}}$ ，所以刚度矩阵 ${{C}}$ 中的各个参数可用工程弹性常数表示为：

$\left\{ {\begin{aligned}& {{c_{11}} = \lambda n\left( {1 - n\mu {'^2}} \right)},\\& {{c_{12}} = \lambda n\left( {\mu + n\mu {'^2}} \right)},\\& {{c_{13}} = \lambda n\mu '\left( {1 + \mu } \right)},\\& {{c_{33}} = \lambda \left( {1 - {\mu ^2}} \right)},\\& {{c_{44}} = G'},\\& {{c_{66}} = E/[2\left( {1 + \mu } \right)]}\end{aligned}} \right.$

(10)

式中， $n = E/E'$ ， $\lambda = E'/(1 + \mu )(1 - \mu - 2n{\mu ^2})$ ， $E$ 、 $E'$ 、 $G$ 、 $\mu $ 、 $\mu '$ 为弹性参数。其中： $E$ 、 $E'$ 分别为同性平面和垂直同性平面的弹性模量； $G'$ 为材料剪切模量； $\mu $ 、 $\mu '$ 分别为同性平面和垂直同性平面施力时的泊松比； $G = E/[2(1 + \mu )]$ 为 $E$ 和 $\mu $ 的函数，不是独立的弹性参数^[18]。

劈裂破坏区的岩体被简化为如下模型：将围岩中贯通的劈裂节理面视为弹性对称面，岩体与劈裂面平行方向上，即竖直方向上为一种力学参数，与节理面垂直方向上则具有另一种力学性质，使用橫观各向同性模型描述劈裂破坏区的力学性质。

3 现场的监测工作

大岗山水电站是近年来修建的坐落于大渡河中游的大型水电项目。水电站一共设置4个机组，每个机组主要由平行分布的3个大型洞室构成，分别为主厂房、主变室和尾水调压室。3个大洞室垂直埋深390～520 m，水平埋深310～530 m；厂房总最大开挖宽度约为30.80 m，开挖高度73.78 m，尾水管底板高程920.02 m；主变室长144.00 m，宽18.80 m，高25.60 m，底板高程962.10 m；尾水调压室长132.00 m，宽24.00 m，高75.08 m，底板高程920.02 m。

水电站厂房围岩主要由花岗岩构成，厂房区域的围岩岩体新鲜较完整，其结构主要为块状–次块状，各岩块之间的嵌合紧密，在厂房的洞室区主要有4条较大的软弱结构面穿过。

厂房洞室局部范围内发育有Ⅳ级结构面，主要由宽度约数厘米的小断层构成^[19]。围岩中主要有4条长度约为3～5 m，个别长度大于5～10 m的裂隙存在。水电站主要洞室的周边围岩大部分为Ⅱ、Ⅲ类围岩，围岩整体比较稳定；洞室的局部围岩为Ⅳ、Ⅴ类围岩，稳定性较差，主要为破碎带发育。经过现场实测及后期反分析计算，构造应力和自重应力构成了大岗山水电站地下厂房洞室区域的主要地应力场，主要被构造应力、地形地貌所影响。根据实测和计算得到，大岗山水电站地下厂房区域的最大主应力为26.9 MPa，属于高地应力地区^[20]。

为能更好地监测施工过程中洞室边墙围岩产生的劈裂破坏及变形情况，使用地质钻机在保证钻孔质量的条件下于主厂房与主变室之间的岩柱上钻取4个水平预埋测试孔，每个孔深46.5 m。4个钻孔主要是利于今后使用全景数字钻孔电视观测围岩劈裂情况等。各监测孔位置情况如表1及图5、6所示。

将A1～A4号测孔两两分组，每组中的两个测孔中分别放置一个滑动测微计和一个综合物探探头，用以进行形变电阻率法和数字钻孔摄像法对主厂房和主变室之间的围岩进行岩体劈裂情况的监测。其中，A3和A4测孔位于3^#和4^#机组母线洞之间的岩柱上。A3孔用于放置滑动测微计，A4孔为综合物探测试孔，两孔孔间距为4.5 m。考虑到主变室已基本施工完毕，在钻孔施工及今后的观测过程中所受的施工干扰较小，因此钻孔由主变室上游边墙向主厂房打钻。

监测围岩变形使用的滑动测微计来自瑞士的Solexperts公司，操作便捷，监测精度和可靠性很高。

围岩出现变形时，由于其物理性质会发生变化，例如裂隙张开或合拢时，会导致围岩自身的电阻率变化。所以通过监测围岩的电阻率变化，可以间接得到岩体裂隙的变化过程，这就是形变–电阻率法。

现场电阻率测试使用WDA–1型超级数字电法仪自制两极探头进行。该仪器探头长150 cm，其上有A和M两个电极，另外两个电极B、N布置在无穷远处。监测过程中，由A、B两个电极产生电流，然后由M和N两个电极测试这一过程中的电压，并通过计算得到岩体电阻率的变化。监测过程中须保证探头在推进过程中能与孔壁很好地接触。

数字式全景钻孔摄像系统采用国内生产的技术产品。工程中的测试孔为水平孔，测试时使用多根测杆将全景摄像探头推入测试孔中，并进行视频录制。

根据滑动测微计、形变电阻率法及钻孔电视测试的结果来看，3种测试手段均能反映出由于开挖扰动引起的围岩松弛劈裂破坏现象。由于各种测试手段的测试精度存在差异，且抗现场施工干扰的能力不同，3种测试手段得出的劈裂破坏区深度略有区别，但总体趋势十分吻合。其中，滑动测微计测试精度高，抗干扰能力强，影响测试结果的因素较少。因此可认为滑动测微计的测试结果较为可靠，其他两种测试手段可作为滑动测微计的有益补充。根据监测，可以得出4^#机组劈裂破坏区深度如表2所示。

由监测结果可知，由于主厂房的洞室尺寸较大，因此在开完后围岩的劈裂破坏区较大，主厂房下游围岩又由于主变室与主厂房两个较大洞室的开挖对两个洞室之间岩柱扰动较大，因此主厂房下游边墙围岩的劈裂破坏区深度大于主厂房上游边墙围岩的劈裂破坏区深度。4^#机组主厂房下游边墙围岩劈裂破坏区平均深度为13～15 m，开挖导致的围岩位移约为17～18 mm。主变室上游边墙围岩劈裂破坏区深度为5～6 m，开挖造成的围岩位移约为5 mm。对于大岗山水电站这种大规模的洞室群来说，该程度的围岩变形是在可以接受的范围内，可以认为支护设计能够满足洞室围岩稳定的要求。

4 计算模型的建立 4.1 围岩劈裂破坏产生的应力状态条件及预测公式

从断裂力学角度出发，基于室内试验，李晓静等^[21]运用能量分析及裂纹扩展等方法、劈裂裂纹产生和扩展的应力条件，假设原生裂纹大多为较大角度，利用滑移型裂纹模型推导得出劈裂破坏出现的判据：

${{{\sigma }}_{1}} \!\ge \!\frac{{{K_{ IC}}\sqrt {{\text{π}}L} }}{{L({{sin}}\;\theta {{co}}{{{s}}^2}\theta \!-\! \mu {{si}}{{{n}}^2}\theta {{cos}}\;\theta )}} \!+\! {{{\sigma }}_{3}}\frac{{{\text{π}} \!+\! ({{sin}}\;\theta {{co}}{{{s}}^2}\theta \!+\! \mu {{co}}{{{s}}^3}\theta )}}{{{{sin}}\;\theta {{co}}{{{s}}^2}\theta \!-\! \mu {{cos}}\;\theta {{si}}{{{n}}^2}\theta }}$

(11)

式中： $\mu $ 为劈裂破坏形成裂纹的摩擦因数； $L$ 为形成劈裂缝的长度，m； $\theta $ 为初始裂纹与水平方向的夹角，（°）； ${K_{ IC}}$ 为断裂韧度，MPa·m^1/2。

这一判据主要是用来预测洞室边墙围岩劈裂破坏区的深度。

4.2 计算流程

1）根据FLAC3D自带模型进行计算，确定卸载时岩体的弹性模量，并对原有的摩尔库伦模型进行二次开发。

2）调用自定义模型，用FLAC3D软件进行计算。

3）在分步开挖时，将劈裂破坏判定准则代入，并将判定劈裂的区域模型转换为横观各向同性模型。

4）计算并提取结果。

4.3 计算模型及相关参数

根据4^#机组实际参数建立包括主厂房、主变室和尾水调压室的机组计算模型，将其导入FLAC3D进行厂房洞室群围岩弹塑性计算，如图7所示。

通过现场采集岩石样本进行加卸载试验，进行反演后得到岩石在正常加载情况下各物理力学参数。根据相关参数和试验结果分段反算出等效系数 $t$ ，进而计算得到卸载时的弹性模量 $\overline E $ 。

因此，在进行数值计算时，4^#机组模型的岩体力学参数如表3所示。

计算模型的范围边界：上游边界从主厂房上游边墙洞壁向上游300 m；下游边界从尾水闸门室下游边墙向下游延伸300 m；顶部界面至地面高程；底部界面下至1 120 m高程（尾水调压室底板高程920.02 m）。计算域内有许多大小断层和裂隙密集带，为便于计算，仅模拟对3个大洞室影响较大的断层和裂隙密集带。

5 数值结果分析和工程验证 5.1 围岩劈裂破坏总体分析

根据现场监测情况及反演计算可得到劈裂破坏判据式(11)中主要参数的取值： $ \mu = 0.2,{{\theta }} = 40{{, }}\;{K_{ IC}} = $ $ 0.84\;{{MPa}} \cdot {{{m}}^{1/2}},L = 5\;{{m}}$ 。将参数代入式（11），并将其导入FLAC3D中，按照第4.2节所述流程进行计算，可得到4^#机组每步开挖结束后，围岩劈裂破坏区范围及深度。虽然主厂房分10步开挖，但第9步和第10步开挖范围较小，对劈裂区范围及围岩的变形影响不大，故选取4^#机组主厂房第8步开挖完成后的位移变化结果进行分析。主厂房第8步开挖结束后，4^#机组3个洞室边墙围岩劈裂破坏区范围见图8。

由图8可知：第8步开挖结束后，主厂房下游边墙的劈裂破坏区深度大于主厂房上游边墙的劈裂破坏区深度；主厂房下游边墙劈裂区最大深度位置出现在第8步开挖后主厂房的中部，最大深度约为17.4 m，而主厂房下游洞室边墙劈裂破坏区的平均深度约为13.6 m，主变室上游洞室边墙围岩破裂破坏区平均深度约6.3 m，计算结果与实际监测结果吻合度较高。

5.2 围岩劈裂破坏总体分析

在4^#机组计算模型的主厂房与主变室之间设置与现场监测位置相同的关键点，见图5。分别使用摩尔库伦本构模型、横观各向同性本构模型及基于能量耗散的橫观各向同性本构模型进行分步开挖计算。计算后提取不同开挖步中各本构中各个关键点的位置值，与现场实测数据进行对比，用以研究计算结果的正确性，对比结果如图9～12所示。

由图9～12可知：开挖时由于主厂房洞室尺寸较大，对围岩造成的开挖扰动影响大，故在靠近主厂房下游临空面的围岩位移较大。第8步开挖结束后，最靠近主厂房下游边墙的测点测得的围岩位移最大，为14.52 mm。随着关键点位置与主厂房下游边墙临空面距离增大，计算得到关键点位移值逐渐减小。当测点距离主厂房下游边墙约15 m时，围岩位移减小趋势减缓，大小趋于稳定。当距离主厂房下游边墙约41 m处时，围岩位移开始逐渐增大，在靠近主变室上游临空面约2 m处，围岩变形较大，第8步开挖结束后，该位置测点检测得到位移值为4.08 mm。

对比考虑能量耗散并使用横观各向同性模型、不考虑能量耗散的横观各向同性模型及仅使用摩尔库伦模型3种不同计算本构模型的模拟结果可知，仅使用摩尔库伦模型计算得到的位移曲线较平缓，不能明显表现出关键点距主厂房下游边墙距离增加，位移先减小后增大的变化趋势。当关键点靠近主厂房下游边墙时，使用该本构模型计算得到的位移值远小于实际监测值；关键点在距离主厂房下游边墙临空面距离超过20 m时，计算得到的位移又大于实际监测值。因此，仅使用摩尔库伦模型计算不能很好地反映地下洞室的劈裂破坏现象。不考虑能量耗散并引入横观各向同性模型计算时，计算得到的位移曲线变化趋势与实际监测值相吻合，但关键点距离主厂房下游边墙临空面超过20 m时，位移计算值普遍比实际监测值大。所以仅使用橫观各向同性模型进行围岩稳定性计算时，虽然能反映洞室围岩的位移变化趋势，但与实际监测情况有较大差别，并不能很好地反映实际情况。在考虑能量耗散并使用横观各向同性模型进行计算时，计算结果不但能反映围岩位移变化趋势，且与实际监测值接近，虽仍有误差，但属于可接受范围内。因此考虑能量耗散并使用横观各向同性本构得到的计算结果更能反映实际情况。

对考虑能量耗散并使用横观各向同性模型的计算结果进行研究可以发现，围岩位移随关键点与主厂房下游边墙距离的变化，分别在距主厂房下游边墙15 m处及距主厂房下游边墙41 m处出现两个变形突变点，即4^#机组主厂房下游边墙围岩劈裂破坏区深度约为14 m，主变室上游洞室边墙围岩劈裂破坏区深度约为5 m。位移计算结果显示的劈裂破坏区深度与监测结果相吻合。

5.3 模型计算结果分析

为进一步得到洞室的变形及应力特征，在4^#号机组计算模型的主厂房拱顶、洞肩、上下游边墙中部位置，主变室拱顶和上下游边墙中部位置洞壁及吊车梁洞壁上选取数个关键点用以计算开挖后洞室各个位置的位移情况，如图13所示。在引入能量耗散及橫观各向同性本构的基础上计算各个关键点在分步开挖时的位移值，如图14～18所示。

洞室开挖后，边墙关键点主要是向开挖临空面变形，即洞室拱顶主要出现竖直位移，洞室边墙的水平方向位移比竖直方向上的位移大。

4^#机组主厂房下游边墙关键点位移大于上游边墙相同位置的关键点位移。4^#机组主厂房下游边墙关键点5（吊车梁位置）和关键点7（主厂房洞中）开挖后洞壁出现的位移较大。根据计算，第8步开挖后，关键点5水平位移为29.46 mm，关键点7水平位移为27.80 mm。主厂房拱顶在第1步开挖结束后竖直方向位移增加到7.32 mm，第8步开挖结束后主厂房拱顶竖直位移量达到10.58 mm。开挖后，主变室边墙围岩关键点位移与主厂房边墙围岩关键点位移相似，4^#机组主变室拱顶在第1步开挖结束后拱顶位置的关键点计算得到的竖直方向上的位移值为6.28 mm；拱顶竖直位移随开挖进行不断缓慢增加，到主厂房第8步开挖结束后，主变室拱顶竖直位移为10.06 mm。主变室上游洞室边墙的水平位移大于下游边墙，第8步开挖结束后关键点9计算得到的水平位移为19.1 mm。

6 结　论

1）引入能量耗散原理，并且在此基础上将原有FLAC3D本构进行二次开发，开发了加、卸荷变弹模本构模型，并将此本构与横观各向同性模型共同应用于高地应力作用下脆性岩体环境地下洞室开挖的稳定性分析计算。

2）为大岗山水电站大型地下洞室群进行了现场监测，运用了钻孔电视、滑动测微计及形变电阻率3种观测方法监测了4^#机组主厂房与主变室周边围岩由于开挖扰动引起的劈裂破坏现象，取得了重要现场实测资料。根据现场监测的结果可知，主厂房上、下游边墙围岩均存在劈裂破坏区，下游边墙围岩的劈裂破坏区深度约为13～15 m，平均深度略大于上游劈裂破坏区深度。主变室上游边墙劈裂破坏区深度约为5～6 m；主厂房下游边墙变形为17～18 mm，主变室上游边墙变形为5 mm。

3）根据实际参数建立了4^#机组包括主厂房、主变室及尾水调压室的计算模型，引入能量耗散原理及横观各向同性本构进行厂房洞室群围岩弹塑性计算分析。可知，主厂房下游边墙围岩的劈裂破坏区深度平均值约为13.6 m，主变室上游洞室边墙围岩劈裂破坏区深度平均值约为6.3 m，计算值与现场监测值的差距较小。在 4^#机组计算模型的主厂房、主变室之间岩柱上设置了与现场监测时位置相同的关键点，利用不同本构模型在FLAC3D中进行了计算分析，并提取不同本构模型相同位置关键点位移进行对比。通过对比分析可知，考虑能量耗散并使用横观各向同性模型进行计算得到的关键点位移曲线与实际监测值曲线较吻合，考虑能量耗散并使用横观各向同性模型能较好地模拟高地应力作用下脆性围岩劈裂破坏的情况。

4）分析计算得到的位移曲线可得，关键点距离主厂房下游边墙15和41 m处分别出现了位移突变点，即主厂房下游洞室边墙劈裂破坏区深度约为15 m，主变室上游洞室边墙劈裂破坏区深度为5 m，位移曲线显示与实际监测情况吻合。