摘要： FEM是目前许多行业数值计算的主要方法。本文结合二滩地下厂房洞室群的围岩稳定分析，主要讨论了进行大型有限元建模及计算时的相关问题，并对洞室群的围岩稳定进行了力学分析，给出了围岩整体是稳定的结论。
关键词： FEM 围岩稳定 相关问题
许多工程分析问题，如固体力学中位移场和应力场分析、热学中的温度场分析、流体力学中的流场分析等，都可借助于计算机数值模拟技术来求出它们的数值解，其中有限单元法（Finite Elements Method，以下简称 FEM）是工程中最常用的数值模拟技术之一。本文结合二滩水电站地下厂房洞室群的围岩稳定分析，对使用有限元软件进行大规模数值模拟计算时的相关问题进行探讨，给出了建议。
1 建模与有限元网格划分 1.1 工程概况二滩地下厂房洞室群设置在左岸下游的山体内，上覆岩体厚度为250~350m，水平埋深300m，洞轴线方向NW6°。洞室群所在区域主要为正长岩和蚀变玄武岩，围岩岩体有80％质量达到A,B级。地下厂房洞室群由三大洞室和母线洞以及尾水洞组成。其中主厂房：长×宽×高＝280.29m×25.5m×65.78m；主变室：214.9m×18.3m×25.0m；尾调室：203.0m×19.8m×69.8m。主厂房主要采用喷锚支护，在主副厂房和尾调室高边墙采用了一定数量的预应力锚索。地下洞室群水平区域布置如图1。
1.2 建模 一般的有限元计算，70%~80%的时间花在建立模型上面，因此选择合适的建立模型的方法是很重要的，否则将有可能事半功倍。
二滩洞室群位于地下，结构复杂，操作起来比较困难。为此把大模型分解为许多很小的模型，划分单元以后再把它们拼装成一个大的模型；并且不能把它们合成一个整体，否则将无法划分单元。具体的作法是：先建立洞室群的模型，然后在长、宽、高方向上分别向外延伸，形成一个长×宽×高=404.2m×283.7m×185.3m的巨大岩体，用该岩体挖去洞室群，就是所需要的模型。由于有限元软件的自适应网格划分功能只能划分规整的模型，如六面体或规则的棱柱体（如三角棱柱体），而地下厂房洞室群则是纵横交错，极其不规则的，建模难度可想而知。经过不断的修改和试验，最终建立的模型共有3647个体组成。
尽管如此，本文模型与洞室的布置相比，还是进行了一些简化，主要有：（1）出于单元形状不致过分奇异以及单元数量方面的考虑，忽略了尾调室和主变室之间的夹角，将它们的轴线设为平行。（2）由于仅考虑三大洞室及其之间的围岩变形与破坏，忽略了进水洞和两个尾水洞。（3）忽略了洞径较小并且仅有局部影响的施工支洞、交通洞、电缆斜井和排水廊道等。
建模过程中，需要注意的有以下几点：
1)对模型内部操作（如空洞）时，可以先生成空洞部分，再和外部模型体组合成所需要的模型。
2)对于任何体的某一个面，都应当是凸多边形（即图形必须在任何一边或其延长线的一侧），且其边数不应超过4，否则将会出现拓扑错误，无法划分单元。
3)合并重复的实体图元，如相同位置的两个关键点或线，以节省存储空间。
4)不能出现扭曲的面，但允许是弧面。
1.3 网格划分
此工程模型中间是地下厂房的洞室群，为此，先给洞室群部分划分单元，然后给厂房邻近的模型体划分单元，这里要注意的一点是要保证它们的接触部分的节点能够吻合。最初采用的是四面体单元，之所以使用该种单元，主要是它的自适应划分能力较强，比较适用于这种复杂的模型。但是当划分网格完成后，单元数量过多（达10万多个），并且出现一些非常小的尖角（小于15度），对应的矩阵难以收敛，故更改为6面体8结点单元。又做了一些简化，并调整单元长度，最终划分单元34224个，结点数量为37632个，计算规模仍然是十分巨大的（划分单元后整体模型如图2，洞室群部分模型如图3）。在划分单元时应注意的问题：
1)必须要保证所划分的单元和邻近体的单元能够吻合，即使不能完全吻合，也要使误差控制在一定的范围内。
2)划分单元完成后，应合并整个模型的所有节点。否则，由于模型是各个独立的体，加荷载计算后，在体与体的结合处将会出现开裂。
3)压缩关键点、线、面、单元、节点等编号，以提高运算效率，节省机时。
4)在应力集中的地方单元网格密度适当加密；反之，则可以加大划分单元的边长，以达到不致太多的影响计算精度，又提高了计算速度。
地下洞室的开挖过程是地应力逐步释放的过程。在本工程中，主要采用应力函数法和有限元分析法相结合[1]。具体做法如下：
第1步：在计算模型边界的某些点处试探施加给定位移得到参考计算结果；然后采用应力函数法在有限元模型的范围内进行回归。回归计算中，顾及上述参考结果以及实测测点应力的最大、最小和平均值，在计算模型的边界上设定若干点的应力分量的值，并使这些边界条件同时参与回归计算；进而得到在建模范围内地应力的第1步近似结果。
第2步：利用第1步的结果计算出在建模边界上所有单元结点的应力，采用类似边界荷载调整法的做法把它们作为远场外荷载施加到有限元模型的边界上。此时由于边界荷载采用的是满足弹性力学中变形协调方程和平衡方程的应力场所求得的，边界上的正应力和剪应力可以共同保证构成平衡力系，并且解决了求剪应力场所遇到的困难。利用上述远场边界荷载，借助有限元程序计算建模范围内的地应力，然后将此结果在实测地应力的测点上与实测地应力进行比较，并且考虑在洞室群范围内地应力的分布情况结合工程概念考察回归计算结果的合理性。根据所得到的趋势，调整在第1步计算中施加在边界上的有限各点的力学边界条件，再次回归计算内部各点的地应力，并由此获得全部边界上由远场位移所引起的所有应力分量。
重复以上两步的分析计算，直到得到满意的应力场为止。
在本次计算中，采用双参数准则中的Drucker-Prager准则来判断复杂应力状态下岩体的破坏。该准则是Mises准则的推广，因此又被称为广义Mises准则[2]。该准则认为，处于弹性区的应力满足式(1)时，则认为岩体处于破坏的临界状态。
（1）
其中：
（2）
（3）
C,φ分别为凝聚力和摩擦角；Rc为单轴抗压强度。
当岩体材料处于弹性阶段时，其本构关系为：
（4）
其中：
为弹性矩阵，
分别为应力矩阵和应变矩阵。
当材料进入塑性以后采用增量理论的本构模型。即认为总应变增量
可以表示为弹性应变增量
与塑性应变增量
之和： 
。
在施工期间，二滩水电站地下厂房的1号尾水洞和2号尾水洞之间曾发生过岩爆及塌方，故取2号机组所在的横断面作为破损区的示意图[3]。由图4可以看出：
1）横断面上破坏主要集中在主厂房下游边墙1号和2号两个母线洞之间，对应的两个尾水洞之间，以及主变室上游边墙对应的两个母线洞之间的部位，这与该处的地应力较大及结构复杂有关。
2）母线洞底板破坏的位置集中在母线洞中部到主厂房一侧。
围岩的垂直位移分布规律十分明显：顶拱向下，底板向上。这是由于开挖引起地应力释放所造成的。围岩总位移的分布规律则基本上是按1号尾水洞到6号尾水洞的方向依次递减（见图5）。各主要洞室的位移分布如下：
1）主厂房上游边墙位移的最大值出现在边墙的中部，下游边墙水平位移的最大值出现边墙上部和顶拱位置。恰好在此位置上，岩爆所造成的破坏最严重。
2）主变室高度较低，约为主厂房高度的2/5，但是其最大位移值却大于主厂房。
3）尾调室的最大位移发生在1号和2号尾水洞之间。
与破损区图相对照，可以看出是否破坏以及破坏程度与绝对位移的大小不一定有确定的对应关系。最明显的是在主厂房的上游拱肩的位移比下游拱肩的位移为小，但是下游拱肩并没有破坏，而上游拱肩则发生了破坏，这主要是由于下游拱肩的刚体位移较大，形变位移却不大的原因。
由以上可看出，地下厂房洞室群的围岩整体是稳定的，破损区是局部的，并且由于该岩体的完整性较好，故可通过适当喷锚处理。围岩的破坏从1号尾水洞到6号尾水洞依次递减，故应加强1号和2号尾水洞之间岩体的稳定性观测，遇到问题及时处理。
图1二滩地下厂房区域水平布置示意图
|
|
图4 2号机组所在面破损区等值线示意图 图5主厂房下游边墙总位移示意图
3.3 计算时应注意的问题
在计算前，可以编制计算文件，计算时直接调入该文件执行即可。例如要模拟开挖的过程就可以编制相应的步骤，在计算不同的开挖步骤时，只需调入不同的文件执行即可。这样就减少了人工干预，提高了计算效率。在计算后，若出现模型开裂的情况，则说明在不同的小模型之间尚未完全合并，需调整合并图元的容差，直至所有小模型合并成一个整体为止。
4 结束语 利用有限元软件进行数值模拟计算已成为当前许多行业数值计算的主流，本文结合具体工程，讨论了在运用此类软件进行大规模数值计算时的相关问题，主要有以下几点：
1）建立模型时，整体考虑选择合适的模型和建模方法，尽量符合工程原型且要便于计算分析。
2）把大模型分解为小模型，做到化繁为简。
3）计算时，可编制计算步骤文件，在不同的载荷步调用不同的文件即可。
参考文献：
[1] 张有天,胡惠昌.地应力场的趋势分析.水利学报.1984,(4).
[2] 徐秉业,刘信声.应用弹塑性力学.北京:清华大学出版社,1995.513~514.
[3] 二滩水电站地下厂房洞室群模拟开挖过程的三维非线性有限元分析.昆明理工大学电力工程学院,2000.