水电站地下厂房结构抗震计算方法探讨

张志国 ，杨 阳 ，牟春来，肖 明

（1.长江设计院，武汉 430010；2.武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072；3.武汉大学 水工岩石力学教育部重点实验室，武汉 430072）

1 引 言

我国西南地区兴建了一批大型水电站地下厂房洞室群。该地区属于地震多发带，地震烈度一般在VII 度以上，地下厂房洞室的抗震性能直接关系到生成设备的正常运行和工作人员的生命安全。对汶川大地震震中区域映秀湾、渔子溪和耿达等水电站地下厂房洞室群的震后调查表明[1]，（1）该地下结构较地面结构有较强的抗震特性，高山峡谷地区水电站厂房应优先布置于山体内；（2）在高达XI 度的震中区域，3个地下厂房均未出现整体围岩失稳，而厂房衬砌等结构却出现明显的开裂和失效，厂房结构是地下厂房洞室抗震设计的关键。

我国《水工建筑物抗震设计规范》[2]水工地下结构部分未对地下厂房建筑物的抗震设计方法作出明确说明，而是建议“在计入结构和围岩相互作用的情况下进行专门研究”。国内各设计单位一般根据工程的实际情况，分别采用拟静力法、反应谱法、动力时程法进行地下厂房结构的抗震计算[3-6]。从3种算法的计算理念及发展过程来看，动力时程法较反应谱法更能合理地模拟围岩与结构的非线性作用过程，而反应谱法较拟静力法更能合理地模拟结构本身的频率特性。拟静力法理论简单，应用时间长，有助于很多具有丰富设计经验的工程师对计算结果的把握；反应谱法是一种频域求解方法，很难准确模拟地下洞室中围岩对混凝土结构抗震特性的影响，而这恰恰是地下洞室结构抗震分析的特点所在。动力时程法是一种时域求解方法，虽然能模拟地震的全过程，但其在实际工程设计中应用时间较短，如地震波荷载时程曲线的确定、地震波荷载的输入、人工边界、围岩与结构的接触、计算结果的分析评价等，均处于快速发展和研究阶段。采用不同的抗震计算方法，可能会对水电站地下厂房结构抗震设计带来一定的影响。

本文采用拟静力法、反应谱法、动力时程法分别对位于汶川地震震中区域的映秀湾地下厂房结构进行动力计算，并将计算结果与震害调查相对比，探讨了3 种计算方法的特点，分析其在水电站地下厂房结构抗震设计中的适用性，以期能对工程设计人员在3 种抗震计算方法的应用上提供些许帮助。

2 3 种方法在地下洞室结构抗震计算中的特点

2.1 拟静力法

地下结构抗震计算的拟静力法一般采用波动场应力法，首先根据指定的波动方程求解介质波动应力场，再考虑波动应力对地下结构的影响。地壳中传播的地震波主要分为纵波和横波，其中纵波产生正应力，横波产生剪应力。复杂的地震波经过傅里叶变换，可表述成多组简谐波的叠加。根据规则的简谐波形，从理论上可推导出弹性介质体中波动应力场[7]：

式中：ρ为密度；kc为峰值加速度放大系数；g为重力加速度；Cp、Cs分别为纵波和横波波速；T0为卓越周期。

设地震波与有限元模型坐标系X、Y、Z 轴的夹角为α、β、γ，则模型任意单元形心波动应力可表示为

式中：

地震波施加到结构上的地震荷载为

式中：[B]为有限元单元应变矩阵。

从波动场应力法的求解过程可以看出，（1）对于地震波荷载，其只针对单一简谐波，并考虑峰值加速度和卓越周期，而未考虑频谱、持时等地震波荷载特性；（2）其将整个地下厂房工程区看为质点，工程区内任意点的波动场应力状态均由式（1）求解，未能反映洞室群围岩波动场的分布特征；（3）其地震动荷载转换为等效静荷载，未能反映地下结构自振特性。但该方法综合考虑了峰值加速度荷载在卓越周期内对地下结构的持续作用，很多工程师认为地下结构拟静力法的计算结果是较动力时程法偏大的，工程依此设计是偏安全的。

2.2 反应谱法

反应谱法有多种计算形式，《水工建筑物抗震设计规范》推荐采用振型分解反应谱法。该方法的理论是：结构物可以简化为多自由度体系，多自由度体系的地震反应可以按振型分解为多个单自由度体系反应的组合，每个单自由度体系的最大反应可以从反应谱求得。

式中：｛φ｝i为结构第i 阶振型，γix、γiy、γiz分别为振型组合系数；uix、uiy、uiz为频率为ωi的单自由度结构在地面加速度、、作用下的反应，可以根据设计位移反应谱或加速度反应谱直接简化计算。

从振型分解反应谱法的计算理论上看，其考虑了结构自身振型特性，通过设计反应谱考虑了地震波形的频谱及幅值特性，但考虑到地下厂房结构的特殊性，在结构自振分析中需考虑围岩对结构振型的影响。目前在众多地下工程反应谱分析中，围岩与结构的作用模拟主要可分为3 类：（1）围岩对地下结构全约束；（2）围岩对地下结构不约束；（3）围岩对地下结构“弱约束”，其中“弱约束”是指通过在围岩与结构间设置1 层材料参数较低的垫层单元，模拟围岩与结构的相互作用。上述3 种方法均很难全面反映围岩与结构的相互作用以及地震作用下围岩对结构的影响。因此，如何考虑围岩与结构的相互作用，是地下洞室结构抗震分析中采用反应谱法计算的关键。

2.3 动力时程法

拟静力法考虑了地震波的幅值特性，反应谱法考虑了地震波的幅值、频谱特性，动力时程法[9]作为一种时域分析方法，则全面考虑了地震波的幅值、频谱和持时3 大特性。该方法将结构抗震、地基抗震等问题统一归结为近场波动问题。从波动学的角度，分析结构作为地震波传播载体的受力特征，通过对工程区波动微分方程进行直接积分求解，模拟地震波在介质体中的传播过程，真实反映地震灾变过程；对地下结构的动本构、围岩与结构的接触算法、地震波的入射、边界条件的设置等提出了较高的要求，但无法给出结构对不同频率波的反应，不利于设计人员对地下厂房结构自振特性的把握。

3 映秀湾地下厂房结构震害调查

2008年四川省汶川大地震，具有震级高（M8.0级）、震源浅（距地表14 km）、破坏性强、波及范围大的特点[10]，对我国西南地区的水电工程造成了较大影响。映秀湾水电站是岷江干流上9个梯级电站之一，上游电站为太平驿水电站，下游为紫坪铺水利枢纽。该站共有3 台发电机组，1971年发电，采用地下厂房的布置形式，埋深约150～200 m。洞室群由主副厂房、主变洞、母线洞、引水洞、尾水洞以及交通洞组成。主厂房的尺寸（长×宽×高）为52.8 m×17.0 m×37.2 m，主变洞的尺寸为59.4 m×7.2 m×27.9 m。洞室群围岩主要为花岗岩及花岗闪长岩，变形模量为10 GPa，泊松比为0.25，黏聚力为2.18 MPa，内摩擦角为41.6°，抗拉强度为1.97 MPa。主厂房洞室采用全断面混凝土衬砌，其中顶拱衬砌厚1.0 m，边墙衬砌厚0.5 m，吊车梁立柱1.0 m×1.5 m，吊车梁为钢梁。发电机层楼板厚0.2 m，横梁0.2 m×0.5 m。

根据汶川地震烈度分布，映秀湾水电站的地震烈度高达XI 度，距震中约8 km，是距离震中最近的水电工程之一，其洞室结构实际发生了一次“原型破坏试验”。将其作为分析对象，研究水电站地下洞室结构的抗震分析方法，具有一定的实践基础。

汶川地震后，笔者曾对映秀湾地下厂房进行震害调查，并参与了后续的厂房加固工作。主厂房结构破坏主要表现为以下几方面。

（1）厂房顶拱衬砌部分在震前就发现多条裂缝，震后部分裂缝宽度及长度有所增加（震前设有测缝计及裂缝标识），并新增加2 条沿上游下游边墙长2～3 m，宽0.5～2 mm 的裂缝（见图1），一条顺机组纵轴线方向长约5 m，宽1 mm 的裂缝。

图1 顶拱开裂Fig.1 Top arch fracture

（2）各洞室两侧边墙衬砌均有明显的纵向裂缝，裂缝延展较长，局部达到5 m 以上，裂缝宽度肉眼可见。

（3）发电机层、水轮机层楼板沿对角出现闭合裂缝，上覆地板隆起（见图2）。发电机层支撑横梁交接处出现表层剥落和局部裂缝。

（4）发电机风罩及机墩裂缝较多。以1#机组风罩为例，裂缝长度一般为0.85～1.70 m，宽0.5～1.8 mm，缝深3.12～7.00 cm。检查表明裂缝均为V形开口，在缝底部尖灭或闭合（见图3）。

（5）各廊道交口处混凝土边墙出现明显的表层剥落、错动开裂，局部钢筋弯曲外露（见图4）。

（6）灾后采用基于环境激励的结构无损动态检测，认为厂房吊车梁立柱及牛腿结构完整，未受损。

图2 发电机层楼板Fig.2 Generator floor

图3 发电机风罩Fig.3 Generator wind cover

图4 洞室交口处Fig.4 Tunnel crossing

映秀水电站地下厂房结构在汶川地震中的震损形式主要是表层脱落和闭合裂缝，局部衬砌结构裂缝延伸较长，宽度肉眼可见，破坏明显。结构破坏形式以拉裂破坏为主，未出现明显的压裂破坏。整体厂房结构在加固后，尚能正常运行。映秀电站的震害调查表明，地下厂房较地面结构有较好的抗震性能，在XI 度的大地震中地下结构尚未全面垮塌失稳，确保了生产人员的生命安全。但对于大型水电站洞室群的建设，其规模非映秀水电站可比，其洞室结构抗震性能依然不容忽视。

4 计算分析及比较

4.1 计算条件

为模拟“汶川”地震中映秀湾水电站主厂房地下结构的动力响应，建立了包括厂房顶拱、边墙衬砌、机墩、梁板柱等在内的3个机组段地下结构有限元模型（见图5、6），共划分96 744个六面体8结点单元，其中Y 轴为厂房轴线，自1#机组指向3#机组为正，X 轴垂直于Y 轴顺水流指向下游，Z 轴与大地坐标系重合。计算混凝土物理参数按照C25选取，抗拉强度取1.25 MPa，抗压强度取12.5 MPa。

图5 有限元整体计算模型Fig.5 Finite element model

图6 厂房混凝土结构模型Fig.6 Concrete structure model

动力计算前首先采用三维弹塑性损伤有限元程序模拟主厂房洞室的开挖过程，将开挖后的围岩扰动应力场作为动力计算的初始地应力场，同时考虑围岩压力、机组运行荷载、自重等因素，计算地下结构在静荷载作用下的应力状态。3 种抗震方法采用统一的结构初始应力状态。动力计算中，围岩采用基于Mohr-Coulomb 屈服准则的弹塑性损伤本构[11]，混凝土结构均采用弹性本构，采用应力强度准则判别破坏，使3 种抗震计算结果有统一的评比指标。为了保障结构动静应力对比的一致性，将反应谱法与时程法求解的应力值按照规范折减0.35，与拟静力法的应力结果相对比。3 种抗震方法具体计算如下。

（1）拟静力法计算条件

采用拟静力计算，需确定地震波的卓越周期 T0和峰值加速度放大系数 kc=αmax/g。根据叶建庆等[11]对汶川地震时P 波和S 波卓越周期的研究，本次计算中取P 波和S 波卓越周期分别为0.15 s 和0.35 s。kc一般可以根据《水工建筑物抗震设计规范中“烈度与地面最大加速度值关系表”确定。但映秀湾水电站的地震烈度为XI 度，超出现有规范量值范围。为保障3 种抗震计算方法荷载量值的一致，本文根据图8 中 Ax(t)、Ay(t)、Az(t)三向加速度时程曲线，计算合成加速度时程 A(t)：

计算可得t=13.06 s时地表最大峰值加速度为10.14 m/s2。同时规范规定，在基岩面以下50 m 及其以下部位设计地震加速度代表值可取为规定值的50%。本文计算中取 kc=0.5×10.14/9.81=0.52。计算程序采用自主开发的地下洞室三维弹塑性损伤有限元计算程序[13]。

（2）反应谱法计算条件

规范规定除重力坝、拱坝外，其余混凝土建筑物设计反应谱最大值的代表值 βmax取2.25。由卧龙波地震最大峰值加速度为10.14 m/s2，可以绘制设计加速度反应谱（见图7）。映秀地下厂房岩体主要为花岗岩，围岩以II 类为主，属I 类场地，反应谱特征周期Tg取0.2 s。计算程序采用ANSYS 的动力响应分析模块。

图7 地震波设计反应谱Fig.7 Design response spectrum

（3）动力时程法计算条件

卧龙强震台是距震中最近的测站，其所测地震波峰值加速度较周边几个测站都大。本次计算选取卧龙地震波EW、NS、UD 三向峰值较大的20 s 加速度时程作为地震动荷载。由于映秀湾地下厂房轴线方位角为108°，需对实测地震波曲线进行坐标转换以及滤波和基线校正，最终可得到施加于模型X、Y、Z 三向的地震波加速度时程（见图8）。计算程序采用自主开发的地下洞室三维弹塑性动力显式有限元计算程序[11]。

图8 汶川地震卧龙台实测加速度时程曲线Fig.8 Measured acceleration time history curves of Wolong in Wenchuan earthquake

4.2 3 种方法计算结果分析

由拟静力法计算得到的映秀湾水电站厂房结构应力分布情况如图9 所示。主厂房结构第一主应力以压应力为主，最大压应力为5 MPa 左右，在C25混凝土的抗压强度范围内，拱座、发电机层楼板梁、风罩与楼板结合部混凝土出现明显的应力集中。主厂房结构第三主应力以拉应力为主，水轮机层以上顶拱以下的边墙衬砌、发电机风罩和吊车梁牛腿部位混凝土拉应力超过1.25 MPa，超出C25 混凝土的抗拉强度范围，最大拉应力达到3.2 MPa 左右。

图9 拟静力法中第一、三主应力图（上游面,单位:MPa）Fig.9 Quasi-static method contour maps of σ1 and σ3 (upstream face,unit:MPa)

根据反应谱法计算得到的厂房结构应力分布情况见图10。由图可以看出，结构应力主要集中在发电机层楼板、横梁、风罩等部位，这些部位最大压应力在0.5 MPa 左右，在C25 混凝土的抗压强度范围内，但拉应力普遍超过1.25 MPa，超出C25 混凝土的抗拉强度范围。下部大体积混凝土、顶拱、边墙衬砌等部位拉、压应力均较小。

根据时程法计算得到的厂房结构应力分布情况见图11。由图可以看出，厂房结构第一主应力量值以压应力为主，并且大部分区域压应力在3 MPa 之内，局部梁板柱交口、风罩交口等部位压应力达到9 MPa，在C25 混凝土的抗压强度范围内结构未出现压裂破坏。厂房结构第三主应力以拉应力为主，水轮机层以上结构整体拉应力较大，大部分区域位于1.2～2.8 MPa 之间，尤其是边墙、风罩、结构交接处等部位拉应力超过1.25 MPa，超出材料的强度承载范围，但吊车梁牛腿部位混凝土拉应力较小，结构未出现拉裂破坏。

4.3 3 种方法计算结果对比

取主厂房结构2#机组段不同特征部位的第一、三主应力值，对比分析3 种算法计算结果并绘制成曲线如图12 所示。图中，特征部位1 指水轮机层楼板；2 指水轮机层边墙；3 指发电机风罩；4 指发电机层楼板；5 指发电机层边墙；6 指吊车梁牛腿；7指拱座；8 指顶拱。

图10 反应谱法第一、三主应力图（上游面,单位:MPa）Fig.10 Response spectrum method contour maps of σ1 and σ3 (upstream face,unit:MPa)

图11 时程法中第一、三主应力图（上游面,单位:MPa）Fig.11 Time-history method contour maps of σ1 and σ3 (upstream face,unit:MPa)

图12 3 种方法第一、三主应力对比图Fig.12 σ1 and σ3 contrast figure in three methods

由3 种方法中计算结果的第一应力量值看，结构应力状态主要为受压，拟静力法计算结果较反应谱法和时程法都大，这是因为拟静力法认为在卓越周期内结构作用地震加速度始终为峰值加速度。从结构压裂破坏的角度看，拟静力法计算结果是偏安全的。从第三主应力量值看，在顶拱以下部位拟静力法计算结果均较反应谱法、时程法大，结果偏安全，但在厂房顶拱部位应力量值较小，较时程法计算结果偏危险。反应谱法计算结果主要表现出发电机楼板、发电机层横梁等结构薄弱部位应力量值较大，而与围岩直接接触的结构应力相对较小。可见，反应谱法对地下结构动力响应薄弱部位的确定有其固有优势，其对薄弱部位的应力计算量值是偏安全的。时程法计算结果除顶拱部位外较拟静力法小，除发电机楼板、横梁等薄弱部位外较反应谱法大，这是因为时程法根据地震波荷载曲线计算，较拟静力法采用固定的峰值加速度，地震荷载相对细化，同时时程法真实反映结构的波动过程，较反应谱法采用各阶振型的最不利组合，更能精细地反映结构的动态受力过程。

4.4 3 种方法计算结果与震后调查对比

根据特征部位第一、三主应力量值判断3 种算法得到的主厂房特征部位破坏情况，并与震后调查结果进行对比见表1。表中，压裂表示特征部位压应力超过混凝土极限抗压强度，拉裂表示特征部位拉应力超过混凝土极限抗拉强。

表1 3 种方法计算结果与震后调查对比Table 1 Results of three calculation methods contrast with the investigation after earthquake disaster

由表1 可知，拟静力法、反应谱法和时程法计算结果均表明主厂房结构特征部位未出现压裂破坏，而是出现不同程度的拉裂破坏，这与震灾调查中“厂房结构未出现压裂破坏”与“结构破坏形式以拉裂破坏为主”较吻合。

拟静力法计算结果显示，水轮机层楼板、水轮机层边墙、发电机风罩、发电机层楼板和发电机层边墙混凝土拉裂，与震灾调查相吻合，但拟静力法显示吊车梁牛腿部位混凝土拉裂，拱座及顶拱未开裂，与灾后调查中吊车梁牛腿完整无损，拱座以上顶拱出现多条裂缝不相符。反应谱法计算结果显示，发电机层楼板、横梁拉裂，吊车梁牛腿未开裂，这与震害调查中“发电机层楼板裂缝，上覆地板隆起，支撑横梁局部裂缝，吊车梁牛腿完整无损”等实际情况相吻合，但其他特征部位均为开裂的结果与震后调查不符。时程法计算结果显示，除吊车梁牛腿外，其他特征部位混凝土均拉裂，这与特征部位震后调查完全吻合，尤其是时程法得到的发电机层楼板第三主应力均在1.25 MPa 以上，且等值线呈现X形，与震后楼板“沿对角出现闭合裂缝”较吻合，这是另外两种方法无法做到的。

总体来看，拟静力法计算结果与震害调查基本吻合，应力量值总体较反应谱法和时程法偏安全，但其对顶拱、牛腿等局部结构存在一定的不合理之处。反应谱法计算结果集中反映结构薄弱部位的地震响应情况，但在顶拱、边墙衬砌、水轮机大体积混凝土等与围岩直接接触的结构地震响应分析时有较大误差。时程法精细反映了地下结构的波动过程，计算结果与震害调查吻合较好。

5 厂房结构抗震计算方法探讨

拟静力法、反应谱法和时程法是水电站厂房结构抗震计算的3 种主要方法，通过数值分析及其与震后调查对比表明，3 种计算方法均能在一定程度上反映动力荷载作用下厂房结构的破坏特性，为抗震设计提供理论依据。

拟静力法理论简单，应用时间长，计算效率高，且基本能反映厂房结构动力破坏情况，但实质上是在结构物表面均匀施加了1 组伪地震荷载，无法精细反映顶拱、吊车梁牛腿等特性结构物在地震中的受力特性。反应谱法考虑了结构自身振型特性，能反映结构薄弱部位的地震响应情况，但其未能全面反映地震过程中围岩与地下结构的相互作用关系，仅表现出围岩对结构的强约束，致使顶拱、边墙衬砌、下部大体积混凝土等地震响应较弱，应力较小，对厂房结构破坏情况模拟效果相对较差。时程法综合考虑了地震波的幅值、频谱和持时三大特性，精细反映了地下结构的波动过程，能准确模拟厂房结构的动力破坏情况，但在实际工程设计中应用时间较短，各种计算理论仍不完善，计算耗时长，尚未得到已建工程震害实例的检验。

3 种抗震计算方法各有利弊，在抗震计算中可综合考虑时间、经济以及抗震工作者的经验等因素选择合适的计算方法。动力时程法作为一种能真实反映厂房结构破坏情况的计算方法值得大力研究和推广，对安全级别高、规模大或处于高地震烈度区的水电站地下厂房结构，建议采用动力时程法进行抗震计算。

6 结 论

（1）映秀湾地下厂房结构未出现压裂破坏，而出现了不同程度的拉裂破坏，3 种算法计算结果与该震害调查相吻合，并能在一定程度上反映动力荷载作用下厂房结构的破坏情况，可为抗震设计提供理论依据。

（2）3 种抗震算法各有利弊，拟静力法计算结果偏安全，但其对局部结构计算有一定不合理之处。反应谱法能展现结构的薄弱部位，便于抗震加固，但在与围岩接触的结构分析中存在较大误差。时程法能精细反映地下结构的地震破坏情况，计算结果与震害调查吻合较好，但其计算过程繁琐，耗时长，并且理论不成熟。

动力时程法作为一种优秀的抗震计算方法，值得研究与推广，建议对安全级别高、规模大或处于高地震烈度地区的地下厂房结构采用该方法进行计算，便于指导其抗震的经济设计。

[1]宋胜武,蒋峰,陈万涛.汶川地震灾区大中型水电工程震损特征初步分析[J].四川水利发电,2009,28(2):1-7.

[2]中国水利水电科学研究院.SL 203-1997 水土建筑物抗震设计规范[S].北京:中国电力出版社,1998.

[3]高峰,石玉成,严松宏,等.隧道的两种减震措施研究[J].岩石力学与工程学报,2005,24(2):222-229.GAO Feng,SHI Yu-cheng,YAN Song-hong,et al.Study of two shock absorption measures in tunnel[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(2):222-229.

[4]孙铁成,高波,叶朝良.地下结构抗震减震措施与研究方法探讨[J].现代隧道技术,2007,44(3):1-5.SUN Tie-cheng,GAO Bo,YE Zhao-liang.Discussion on anti-seismic and seismic-relieving measures and corresponding research methods for underground structures[J].Modern Tunnelling Technology,2007,44(3):1-5.

[5]章伟.近年日本的隧道施工与抗震设计[J].世界地震工程,2006,22(2):68-71.ZHANG Wei.Current construction and earthquakeresistant design for tunnels in Japan[J].World Earthquake Engineering,2006,22(2):68-71.

[6]王明年,崔光耀.高烈度地震区隧道设置减震层的减震原理研究[J].土木工程学报,2011,44(8):126-131.WANG Ming-nian CUI Guang-yao.Study of the mechanism of shock absorption layer in the supporting system of tunnels in highly seismic areas[J].China Civil Engineering Journal,2011,44(8):126-131.

[7]曹善安.地下结构力学[M].大连:大连工学院出版社,1987.

[8]朱伯芳.有限单元法原理与应用（第三版）[M].北京:中国水利水电出版社,2009.

[9]ANIL K,CHOPRA.Dynamics of structures:Theory and applications to earthquake engineering[M].Upper Saddle River:Prentice Hall,2011.

[10]晏志勇,王斌,周建平,等.汶川地震灾区大中型水电工程震损调查与分析[M].北京:中国水利水电出版社,2009.

[11]张志国,肖明,陈俊涛.大型地下洞室地震灾变过程三维动力有限元模拟[J].岩石力学与工程学报,2011,30(3):509-523.ZHANG Zhi-guo,XIAO Ming,CHEN Jun-tao.Simulation of earthquake disaster process of large-scale underground caverns using three-dimensional dynamic finite element method[J].Chinese Journal of RockMechanics and Engineering,2011,30(3):509-523.

[12]叶建庆,苏金蓉,陈慧.汶川8.0 级地震动卓越周期分析[J].地震研究,2008,31(增刊):498-504.YE Jian-qing,SU Jin-rong,CHEN Hui.Predominant period of the Wenchuan Ms 8.0 earthquake motion[J].Journal of Seismological Research,2008,31(Supp.):498-504.

[13]俞裕泰,肖明.大型地下洞室围岩稳定三维弹塑性有限元分析[J].岩石力学与工程学报,1987,6(1):47-56.YU Yu-tai,XIAO Ming.Three dimensions elasto-plastic finite element analysis for the surrounding rock stability of large-scale underground openings[J].Chinese Journal of Rock Mechanics and Engineering,1987,6(1):47-56.