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