1　前言

    　某抽水蓄能电站位于我国东北地区，枢纽建筑物主要由主坝、副坝、输水系统和地下厂房等组成。地下厂房内安装四台单机容量为300MW的水泵水轮发电机组，发电最小净水头403.40m，发电最大净水头447.30m，抽水最大扬程455.90m，抽水最小扬程416.60m；发电最大引用流量4×77.51m3/s,抽水最大流量4×68.0m3/s。厂房洞室系统是一组空间立体交叉的地下洞室群。厂房洞室系统布置以厂房、主变室、尾闸室为核心，三大洞室从上游向下游依次平行布置，厂房各主要洞室的洞轴线方位角均确定为NW311°。厂房从右到左依次为安装间、主机间、副厂房。主厂房洞室尺寸150.50m×24.00m×53.40m，副厂房洞室尺寸24.00m×24.00m×31.50m，主变室洞室尺寸114.55m×20.00m×22.70m，尾闸室洞室尺寸100.0m×7.0m×28.00m，主变室轴线与厂房轴线间距61m，尾闸室轴线与主变室轴线间距41.00m。厂房主机间与主变室之间的岩体开设母线敷设通道，布置了四条垂直厂轴线的母线洞，母线洞洞室断面为8.20m×9.00m，洞室轴线间距22.50m。

   　 厂房系统位于输水隧洞中部的山体内，轴向N W 311°，埋深300m～310m，地下水位高出厂房洞室拱顶290m左右，围岩为新鲜白岗花岗岩，岩质坚硬、完整。岩体中节理不甚发育，主要为走向N10°～30°W，倾向NE或SW，倾角70°～85°和走向N55°～80°W，倾向NE或SW，倾角70°～85°两组，间距一般为1m～2m，大者达5m左右，多呈闭合状态。缓倾角节理只偶有分布，但延伸不长。厂址地段主要通过有f31、f32、f33、f34四条小断层破碎带，其中f33、f34为裂隙状小断层，宽仅3mm～5mm，延伸不长。

    　厂房部位地应力测试结果为：最大主应力12.2MPa，方向N71°W，与厂轴夹角22°，近水平，最小主应力5.7MPa。岩石抗拉强度7.95MPa，抗压强度106MPa。

    2　研究目的及内容

    　地下厂房系统施工期围岩的整体安全稳定问题是控制水利水电工程建设的关键问题之一，国内外大型地下洞室群都要采用地下厂房岩石力学试验、地应力测试、试验洞现场观测、工程岩体分级、有限元数值计算等多种手段，对地下厂房围岩稳定性问题进行综合分析和研究，其中整体围岩稳定弹塑性有限元计算方法已经逐渐成为地下厂房围岩稳定分析的一个必选项目。

    　根据地下厂房系统的结构形式、地质条件和施工运行特点，该地下洞室群的整体围岩稳定计算将重点研究地下厂房等大跨度洞室高边墙、主厂房与主变隔墙、顶拱、洞室交叉口、断层、地质构造带和其它重要部位岩体结构的安全稳定问题，对其应力应变和变形的分布特征、施工成洞条件、整体安全稳定性进行分析评价，进行地下洞室群围岩变形预测研究。

    　计算时重点模拟地下厂房系统的主要结构形式、厂区附近的主要地质构造，考虑初始地应力对围岩稳定的影响，考虑对围岩稳定有利的各洞室开挖方式及开挖顺序，通过选取不同的地质参数进行洞室围岩稳定的对比分析。

    3  研究方法

   　 ANSYS程序功能完善，成果可靠，是国际上应用最广泛、最具有权威性，并通过ISO9001认证的大型有限元计算软件。为用户提供了一百多种单元，所建模型能较全面地反映真实的工作状态，几乎能覆盖所有工程问题。ANSYS程序可研究多种材料非线性，如：塑性、蠕变、粘塑性等。

    　非线性有限元增量形式的基本方程为:

    [k]{Δδ}i={ΔR}i+{ΔRp}+{R}　　　　　　　　　　　　　(1)

    式中：{R}为迭代过程中产生的不平衡力,{ΔR}i为荷载增量,{ΔRp}为非线性等效结点荷载:

    {ΔRp}=Σe∫Ωe[Ce][B]T{Δσp}dv

    {Δσp}=[Dp]{Δε}i-1            　　　　　　　　　　  (2)

    由此可得:

    {δ}i={δ}i-1+{Δδ}i

    {σ}i={σ}i-1+{Δσ}i       　     　　　　　　　　　　 (3)

    由（1）（2）（3），可通过多次迭代，求得非线性有限元的单元应力{σ}、应变{ε}和结点位移{δ}。

    岩石按理想弹塑性材料考虑，使用Drucker-Prager屈服准则，考虑了由屈服而引起的体积膨胀。其等效应力的表达式为:

    σe＝3βσm+                              (4)

    其中: σm＝平均应力或静水压力＝

    {S}＝偏差应力

    β＝材料常数

    [M]＝Mises屈服准则中的[M]

    　上面的屈服准则是一种经过修正的Mises屈服准则，它考虑了静水应力分量的影响，静水应力（侧限压力）越高，则屈服强度越大。

    材料常数β的表达式如下：

                                                   (5)

    材料的屈服参数定义为

                                                 (6)

    其中：C为凝聚力，φ＝内摩擦角（DP材料的输入值）。

    屈服准则的表达式如下：

               (7)

   　 采用“死单元”功能来模拟洞室的开挖过程。挖去的单元将成为“死单元”，程序将其刚度矩阵变成一个很小的值。“死单元”的单元载荷和质量将设为零，因此不对载荷向量生效，其质量和能量也不包括在模型求解结果中。

4 计算模型

    　根据厂房系统的结构和厂区地质资料，建立空间几何模型。模型包括主厂房洞室、主变室、尾闸室、母线洞、尾水洞、交通洞以及四条主要断层构造带。

   　 计算范围：上下游方向取厂房上游边墙到尾闸室下游边墙距离的5倍以上，铅直方向取厂房高度的6倍以上，沿厂房轴线方向取厂房长度的4倍左右。X轴表示上、下游方向，指向下游为正，Y轴表示铅直方向，向上为正，Z轴表示厂房轴线方向，满足右手螺旋法则。

   　 开挖部分的有限元网格如图1所示，四条主要断层f31、f32、f33和f34如图2所示。

 
图1  开挖部分的有限元网格图
 
图2  厂区内四条主要断层示意图

    5 边界条件及载荷

    　厂房系统埋深约300m，由于计算区域顶部未取到山顶，因此将未取的山体重力荷载施加于模型顶部，在计算区域的底部施加Y方向约束。根据初始地应力资料，分别在X和Z两个方向施加侧向压力，两个方向的侧压力系数为0.7、1.0。厂房系统未开挖时，整个结构处于初始地应力场中，承受自重荷载及横向构造地应力。

    　根据厂房系统的施工设计方案，主厂房洞室分7层开挖，主变室分3层开挖，尾闸室分6层开挖。三维弹塑性计算模拟9个开挖载荷步：

    第一步：形成厂房周围的初始地应力场

    第二步：交通洞和尾水洞开挖

    第三步：主厂房第1层与主变室第1层同时开挖

    第四步：主厂房第2层、主变室第2层和尾闸室第1层开挖

    第五步：主厂房第3层、主变室第3层和尾闸室第2层开挖

    第六步：主厂房第4层、母线洞和尾闸室第3层开挖

    第七步：主厂房第5层、尾闸室第4层开挖

    第八步：主厂房第6层、尾闸室第5层开挖

    第九步：主厂房第7层、尾闸室第6层开挖

    各载荷步的开挖高程如图3所示。

 
图3  载荷步示意图

6 计算参数

   　 为了便于模拟，将岩体和断层假定为均匀的连续体。考虑到岩体的粗糙和不均匀性，计算时采用了如表1所示的物理力学参数。

表1  计算采用的物理力学参数

    7 计算结果及分析

    7.1　位移

    　厂房系统的开挖引起周边围岩的各部位产生变形，各洞室顶拱均有不同程度的下沉，同时各洞室底板上抬。随着厂房系统的逐步开挖，顶拱下沉及底板上抬呈增大趋势。整个开挖过程中，主厂房顶拱最大下沉量为6.0mm，底板上抬7.6mm；主变室顶拱最大下沉量5.12mm，底板上抬7.7mm；尾闸室顶拱最大下沉量为2.12mm，底板上抬1.2mm。各洞室边墙在开挖过程中，总的位移趋势是向洞内收敛变形。主厂房边墙向洞内收敛变形值为7.26mm。由于断层切过各洞室，因此边墙与断层的交汇处均有较大位移值。在断层经过处，主变室下游边墙的最大位移达26.9mm，尾闸室上游边墙最大位移值12.77mm。母线洞的开挖对各洞室的边墙位移有一定的影响，位移出现回弹，特别是主变室和尾闸室的上游边墙，位移回弹比较明显。开挖过程中，各洞室顶拱、底板、边墙及端墙的位移如表2所示，单位mm。

    表2　各洞室顶拱、底板、边墙及端墙的位移

    　在＃3机组附近切取一个典型断面。开挖全部完成后，典型断面Ｘ方向的位移等值线如图4所示。由图中可以看出，各洞室边墙均向洞内收敛变形，其中：主厂房5.0mm，主变室6.0mm，尾闸室4.0mm。主变室下游边墙与断层f32交汇处位移较大。

 
图4　典型断面的Ｘ方向位移等值线

  　  随着开挖的深入，各洞室顶拱逐步下沉，底板上抬。开挖全部完成后，典型断面Ｙ方向位移等值线如图5所示。此时主厂房顶拱下沉4.0mm，主变室顶拱下沉4.0mm，尾闸室顶拱下沉1.0mm；而主厂房与尾水洞挖通后底板的上抬位移值稍有减小，由挖通前的6.0mm减为5.0mm。图中明显可见断层对岩体位移的影响，位移等值线在断层经过处不连续。

 
图5　典型断面的Ｙ方向位移等值线

    7.2　应力

    　在自重荷载及横向构造地应力作用下，厂房区域的初始压应力最大值为10.0MPa。交通洞与尾水洞开挖后，尾水洞周围的最大压应力增至16.0MPa。主厂房与主变室开挖后，两洞周围及洞间的压应力增大，主变室底部的压应力增加到26.0MPa。母线洞开挖后，母线洞与主厂房、母线洞与主变室交叉口及母线洞底部压应力较大。开挖全部完成后，主厂房与尾水洞交叉口、母线洞与主变室交叉口处压应力较大，最大值达30.0MPa，如图6所示。

 
图6　开挖完成后典型断面的σ3分布云图

　　厂房区域的围岩初始应力均为压应力。尾水洞开挖后，厂房区域仍未出现拉应力。主厂房与主变室开挖后，两洞室顶部和底部出现拉应力区，最大值1.5MPa。尾闸室开挖后，其底部也出现0.6MPa的拉应力区。至第五载荷步，主厂房底部的拉应力区最大，最大拉应力值达1.4MPa，如图7所示。第六载荷步，即母线洞开挖后，主厂房底部的拉应力明显减小。厂房系统全部开挖完成后，由于各洞室之间的连通，主厂房和尾闸室周围的拉应力已经基本释放。

 
图7　第五载荷步典型断面的σ1分布云图

    7.3　塑性区

    　当交通洞与尾水洞开挖时，尾水洞附近出现塑性区。主厂房和主变室开挖后，在两洞室顶部和底部出现塑性区。尾闸室开挖后，其底部也出现塑性区。随着开挖的深入，塑性区域继续发展。开挖全部完成后，洞室周围的塑性区深度约10m。开挖完成后，典型断面的塑性区分布如图8所示。

    8 结语

    　由于该地下厂房厂址周围的岩体坚硬，完整性好，因此厂房系统的开挖对围岩稳定的影响不大。

    　从位移的空间分布规律上看，断层经过处围岩上的位移等值线均发生间断。由于断层f31和f32同时切过主变室，又在主变室拱肩处交汇，因此对主变室的位移有很大影响，特别是主变室下游边墙的位移受其影响较大。断层f31较宽，又在中部横切尾闸室，因此尾闸室的位移受断层影响亦非常明显。主厂房附近的断层f33和f34较小，且延伸不长，对主厂房的位移有一定影响。虽然断层对洞室的位移有较大的影响，但总体说来顶拱与边墙的位移值均未超过允许值。

    　开挖对厂区围岩的影响非常明显，洞室周围的压应力变化较大，特别是主厂房与尾水洞交叉口、母线洞与主变室交叉口处压应力较大，但远未达到岩体的抗压强度。断层f31和f32在主变室下游拱肩处交汇，使得开挖后该部位产生了7.6MPa的拉应力，未超过岩体的抗拉强度。

    　在断层f33经过处，主厂房顶拱及下游边墙的局部压应力已超过了岩体的抗压强度,此处岩体有可能受压破坏，在施工时应注意预防，及时采取措施。

    　由于开挖施工和断层带的影响，各洞室周围产生了一定范围的塑性区域，在支护设计中应予以考虑。(E-works)