基于单元重构与节点分离的大型地下洞室软弱结构面模拟方法

赵健+肖明+陈俊涛+李冬冬

摘 要：针对在大型地下洞室中软弱结构面难以合理模拟的问题，通过引进“单元重构节点分离”思想，提出了新的非连续软弱结构面模拟方法.将该方法用于实际地下厂房模型中并模拟出非连续结构面，并将该重构模型导入计算程序，与常规连续有限元模型的计算结果进行对比.结果表明，重构模型能够用于数值计算，计算结果在非连续结构面附近呈现非连续性，且计算结果所反映的规律与所认识规律基本一致，能够反映真实的结构面赋存状态.因此，该软弱结构面建模方法实现了结构面的非连续化，为大型地下洞室的复杂地质非连续结构面模拟提供了新的实现途径.

关键词：单元重构；节点分离；结构面；非连续；建模方法

中图分类号：TU45 文献标志码：A

文章编号：1674-2974（2017）03-0134-09DOI：10.16339/j.cnki.hdxbzkb.2017.03.017

Abstract：As weak structural planes are difficult to be simulated reasonably in large underground chamber， “element reconstruction and node separation” was introduced， and the new simulation methodology of discontinuous weak structural planes was proposed in this study. The proposed methodology was used in the actual numerical model of underground powerhouse to simulate the discontinuous weak structural planes. The reconstructed model was then imported into calculation program， and the result was compared with the prediction by conventional continuous finite element model. The comparison shows that the reconstructed model can be used in the numerical calculation， and the discontinuity feature is presented near the discontinuous structural planes. Moreover， the rules in the calculation results are consistent with the general rules， which can describe the actual occurrence mode of the structural planes. Therefore， the simulation methodology of discontinuous weak structural planes achieves the transformation into the discontinuity of the structural planes， and it provides a new way for the simulation of complex geological discontinuous weak structural planes in large underground chamber.

Key words：element reconstruction；node separation；structural planes；discontinuity；modelling

軟弱结构面是力学强度明显低于围岩且一般充填有一定厚度软弱物质的地质界面.在地下洞室开挖与施工中，软弱结构面是控制围岩稳定的关键因素，它直接制约着工程岩体变形、破坏的发生和发展过程，岩体的失稳也往往发生在这类岩体强度最薄弱的部位[1].同时大型地下洞室往往交错纵横，结构布置复杂，软弱结构面与洞室任意切割，在建模前期直接建立结构面单元也较为复杂.因此软弱结构面的正确模拟是分析地下洞室围岩稳定的关键，是确保围岩稳定计算结果合理的前提条件.

目前针对软弱结构面的模拟，主要有以下几种常见的方法：①单元弱化.将结构面穿过的有限单元的力学参数降低而不建立单独的断层单元，对其进行单独赋值处理，以此将结构面单元和周围岩体区分开来.这种方法计算思路较为简便，且容易操作，但是容易得到齿状交错的断层单元，与实际工程状态有差异.②长而薄的有限单元[2].直接在计算模型中依据软弱结构面产状建立有限单元模拟结构面，并用有限元方法进行分析.这种方法能够得到较为准确的、符合实际的平面断层单元，但是前期建模较为复杂，建模效率较低.③隐含断层[3].肖明等采用隐含断层单元的数值模拟计算方法将复杂断层单元隐含在岩体单元中，用一个均质各向异性的岩体当量复合单元模拟断层的结构效应.这种方法能够将断层效果耦合到岩体单元中，力学性质简明，但是需要重新编制有限元程序，不能直接应用到现行的有限元计算程序和软件中.总结方法①②③，可发现目前这几种常见的结构面模拟方法仍然基于变形连续的有限元分析，不能反映由断层所带来的非连续特征，结构面两侧岩体仍然连续，没有断开，与实际情况有出入[4].除了以上几种方法，应用比较广泛的还有Goodman单元[5].Goodman将其提出并应用于岩石力学中作为节理单元，后来被推广应用于土与结构的共同作用、人工块体结构的有限元计算中作为接触面单元.这种方法能够模拟结构面的非连续变形特性，但是单元无厚度，结构面两侧单元容易相互嵌入.

在总结以往研究成果的基础上，本文引进“单元重构节点分离”思想，改进了软弱结构面模拟方法.在不考虑结构面的有限元模型的基础上引入结构面，将结构面穿过的有限单元进行切割重构并生成新单元，同时将结构面两侧共用节点进行分离.这种方法不但使大型地下洞室群中复杂软弱结构面建模简单化，而且能够实现结构面的非连续化，为含有复杂结构面的地下洞室围岩稳定分析提供了一种十分有效的新思路.

1 软弱结构面模拟过程

大型地下洞室软弱结构面的建模可以归纳为以下5个步骤， 如图1所示.

1）不考虑结构面完成单元离散.根据地下结构的布置，考虑各管道和洞室之间的相对位置，建立不考虑结构面的有限元模型，完成单元离散.这一部分的发展已经比较成熟，许多有限元软件都可以完成模型的建立，并且不同软件建立的模型也可以通过编制程序完成相互转换.

2）确定结构面的信息.根据不同的发育程度可将结构面划分成Ⅰ～Ⅴ级，发育程度从区域性的断裂带逐步减小到微小的节理裂隙，由此带来了结构面模拟方法的不同.根据结构面相对于地下厂房的延伸范围，可以分为有限结构面和无限结构面；根据结构面的分布情况，可以分为确定结构面和随机结构面.目前很多学者将结构面模拟为具有圆盘形状的平面，同时通过Monte Carlo随机方法[6]也可以生成具有裂隙特征的结构面信息.

3）采用单元重构技术，完成结构面导入.由于地下洞室存在的结构面数量往往不止一条，因此可以采用分次导入的方式.在不考虑结构面的有限元模型的基础上，导入第一条结构面切割单元并生成重构模型，之后依次在之前的重构模型的基础上，导入剩余的结构面.同时值得注意的是，在第2）步的结构面信息确定中，当某条结构面发育程度较高，结构面间空隙较大且软弱填充物较多时，可以确定结构面间隙的相对距离，并通过平移结构面，对有限元模型进行两次切割，从而模拟出结构面之间空隙的效果，中间的狭长单元可视为软弱填充物.

4）采用节点分离技术，完成结构面非连续化.记录在单元重构中新生成的节点，并通过添加重复节点完成结构面的非连续化.该步骤是区分重构模型与传统有限元连续模型的最根本的技术.在有限元模型变形之前，结构面上重复节点的空间坐标相同，重复节点重合.当结构面两侧单元发生相对变形时，重复节点的空间位置发生相对错动，结构面脱开分离.

5）检查重构单元的几何信息.在完成单元重构和节点分离后，应重点检查重构单元的几何信息，防止出现一些几何形态畸形的单元.可通过计算单元的雅克比矩阵或单元体积来判断单元是否畸形.

以下重点介绍在模拟过程中涉及到的“单元重构”和“节点分离”技术.

2 单元重构与节点分离

2.1 单元重构

关于单元重构方法的介绍详见文献[7-8]，其基本思想是对结构面穿过的单元进行二次划分，生成新单元.首先，采用成熟的网格划分技术，对研究区域进行单元离散，重点研究区域和开挖临空面附近的单元划分相对于其他区域应該较小，以满足研究精度的需要，以此生成不考虑结构面的计算模型；其次，将结构面依次导入计算模型，对结构面穿过的单元进行二次划分，生成满足计算要求的新单元.

为了形象具体地说明该算法，这里给出一个简单易懂的二维实例.首先对该区域进行网格离散，生成4单元9节点有限元模型，如图2所示；然后将结构面1和结构面2导入，对结构面穿过单元进行一次单元重构，单元E1被重新划分为单元E1和E2，单元E2被重新划分为单元E3和E4，单元E3被重新划分为单元E5，E6，E7和E8，单元E4被重新划分为单元E9和E10，如图3所示.

将结构面导入后，在结构面附近会生成许多形态不规则单元，例如图3中的五边形单元E1和E4.目前在二维有限元计算中，比较常见的单元形态为四面体单元和三角形单元，这2种单元形态简单，适应性良好且能够满足工程计算精度需要，而其他形态的单元在有限元计算中往往比较复杂，其形函数很难构造[9-10]，因此对这部分单元需要进行二次单元重构.

将图3中的单元E1和E4添加网格线进行二次单元重构，得到如图4中所示结果.五边形单元E1被重新划分为三角形单元E1和四边形单元E2，五边形单元E4被重新划分为四边形单元E5和三角形单元E6.至此，图2中所有单元被重新切割划分为三角形和四边形单元，重构生成的单元都能进行有限元分析计算.

同理在三维分析中，比较常见的计算单元为八节点六面体单元和四面体单元，其他的复杂形态单元也可通过添加网格面等方式重新划分为四面体单元[11]，图5和图6显示了五节点四棱锥和六节点三棱柱的重构过程，其他复杂单元也可参考下图重构.

结构面导入模型后，由于只在结构面附近进行单元重构，离结构面较远的大部分围岩单元仍然保持原单元形态，因此需要进行重构的单元数量相对整体模型来说较少，计算量不大.且单元重构完成后模型中只存在八节点六面体和四面体2种单元形态，这2种单元形函数表达方式都比较简单[9]，且能够满足工程计算精度要求.

采用该单元重构方法，能够减小含有结构面模型的建模难度，可以在不考虑结构面的有限元模型基础上任意导入结构面.该重构模型能够保留在关键研究部位和开挖临空面附近的细小划分单元，同时能够获得结构面附近的重构单元，这一部分单元往往划分尺寸较小，能够满足研究精度要求.但是用此方法重新生成的计算模型仍然是连续的，结构面两侧岩体仍然共用节点，没有实现节点分离，计算结果不能较好地呈现非连续特征.

2.2 节点分离

由于采用单元重构技术后的计算模型仍然是连续的，因此需要进行节点分离，其基本步骤如下：

1）遍历模型中的所有结构面，存储每个结构面附近的所有节点.

2）查找步骤1）中存储节点所在的所有单元，并分析是否存在单元分布于结构面两侧.若存在，则记录该共用节点和单元.

3）计算共用节点所涉及的单元数量，在相应的单元增加重复节点，该重复节点的空间坐标在发生相对位移前相同.

以下给出了一个二维节点分离的简单示例，如图7所示.结构面附近的节点编号为3， 5， 7， 10， 11， 12， 13和15，将其进行节点分离.除了节点11外，其余7个节点都分离为2个节点，由于节点11处在2个结构面的交界处，因此节点11分离为4个节点.

通过增加共用节点，完成了节点分离，实现了结构面的非连续化.

3 地下洞室断层模拟

本节以镇安水电站地下厂房为例，将单元重构与节点分离技术应用到地下洞室断层模拟中.

镇安水电站地下厂房模型高为622 m，宽为300 m，地下洞室群埋深约为210 m，根据实测地形资料模拟山体表面起伏.镇安水电站地质勘探资料见表1，地下厂房区域断层发育较为充分，主要存在表1所示2组断层，而其他裂隙节理等发育不是非常明显且数量不多，因此主要对以断层为主的结构面进行模拟.

以厂房1#机组为坐标原点建立局部坐标系，以厂房轴向指向4#机组段方向为y方向，以垂直y方向指向下游为x方向.

将断层模拟成具有一定范围的圆盘面[12]，由断层产状和断层圆心点可唯一确定断层的空间位置，同时由断层半径可确定断层的延伸范围.由于断层间裂隙较明显，厚度较大，因此对该模型中的断层采用2个平行面进行切割，模拟出断层厚度，中间切割重构而成的单元可看成是断层的软弱夹泥.

根据表1中的断层产状，将断层导入有限元模型中，经过单元重构和节点分离，完成了断层的模拟，如图8，图9和图10所示.从图中可看出，采用“单元重构节点分离”方法，能够导入任意产状结构面，减小了前期结构面建模的难度，且结构面附近单元相对密集，重构单元形态良好，能够满足工程需要.

4 考虑软弱结构面的重构模型数值计算

为验证所建重构模型能够用于数值计算，将其采用NDDA方法（基于节点的非连续变形分析）[13-14]对第3节中考虑结构面的镇安地下厂房模型进行分析，并与传统连续断层模型的有限元计算结果进行对比.传统连续断层采用真实存在的长而薄的薄层单元来模拟，地下洞室岩体力学参数见表2.

从图9（b）可看出，2条断层都穿过2#机组段，2#机组段受断层切割作用明显，因此选择2#机组段作为典型机组段进行分析.计算结果如图11和图12所示，分别表示了2#机组段采用连续断层和非连续断层的洞周位移和最大主应力分布.规定压应力为负，拉应力为正.

当采用连续断层时，断层附近的应力值有所降低，断层附近应力云图沿断层向下突出，意味着应力值由正常范围向断层带降低[15].其计算结果往往存在一个“影响带”，该“影响带”区域往往大于断层本身厚度.由于结构面的作用，该区域位移增大，应力强度降低，但仍然保持连续变化，该规律从图11（b）和图12（b）可以清晰地看出.而当采用非连续断层后，断层两边的位移和应力状态往往不连续，两侧差异明显，出现齿状和交错的现象[16].

为了更清晰、确切地观察连续断层和非连续断层对洞室稳定影响的区别，选取受断层切割的主廠房上游边墙底部作为研究对象.主厂房上游边墙中下部高边墙效应明显，在洞室围岩稳定分析中是重点研究部位.

从图13和图14可看出，采用非连续断层时，断层单元的位移和应力状态与两侧岩体的位移和应力状态完全不同，呈现出非连续变化[17]，断层与两侧岩体的位移和应力值相比相差较多，且断层和两侧岩体出现错开现象（如图13（a）所示），而采用连续断层时，断层与两侧岩体的位移和应力值连续变化，没有出现位移和应力值的不连续跳跃现象.从数值上看，非连续断层和连续断层相比，断层位移增幅为0.3～0.5 cm，断层最大主应力降幅在2.0 MPa左右.同时可观察到，在采用2种不同模拟断层情况下，“影响带”范围之外的岩体位移和应力状态的变化规律基本一致，没有较大差别.

可以看出，考虑非连续断层后，洞周围岩变形呈现不连续性，计算结果所反映的规律与所认识的经验规律一致，表明含有非连续断层的重构模型可以用于数值计算，能够反映非连续断层对于围岩稳定的影响.

5 结 论

1）基于单元重构提出了以增加共用节点为核心的节点分离技术，完成了软弱结构面的模拟.该方法的突出优势在于能够在传统有限元模型的基础上任意导入软弱结构面，而非前期直接模拟结构面，很大程度上减小了结构面模拟的难度，并且实现了模型从连续到不连续的转换.相比于传统的结构面建模方法，该方法能够合理地模拟模型中结构面位置，正确地反映结构面的实际空间分布，形成非连续结构面，更符合工程实际，是一种全新的、精确的结构面建模方法.

2）以地下洞室为例，将提出的结构面建模方法应用于有限元模型中，并采用NDDA方法进行数值计算，与连续有限元模型的计算结果进行对比.结果表明，当采用非连续结构面时，结构面与周围岩体的位移和应力状态呈现非连续变化规律，其数值结果和连续结构面相比也有变化.计算结果所反映的规律与所认识规律基本一致，更接近工程实际，证明了该非连续模型用于实际工程的可靠性.

由此可以看出，以“单元重构节点分离”为核心的地下厂房软弱结构面建模方法是一种全新的结构面模拟方法，该方法实现了结构面的非连续化，并为结构面的非连续变形分析提供了新的实现途径.

参考文献

[1] 康石磊， 杨峰， 张箭，等. 基于强度折减和上限有限元的椭圆形毛洞隧道围岩稳定性分析[J]. 湖南大学学报：自然科学版， 2015， 42（9）：104-109.

KANG Shilei， YANG Feng， ZHANG Jian， et al. Finite element upper bound analysis of stability of unlined elliptical tunnel based on strength reduction method [J]. Journal of Hunan University： Natural Sciences， 2015， 42（9）： 104-109.（In Chinese）

[2] 邱颖， 孙云贵. 断层力学模型及其数值模拟方法评述[J]. 水利与建筑工程学报， 2004， 2（1）：14-16.

QIU Ying， SUN Yungui. Study on fault mechanics model and numerical simulation methods[J]. Journal of Water Resources & Architectural Engineering， 2004， 2（1）：14-16.（In Chinese）

[3] 肖明， 陈俊涛. 大型地下洞室复杂地质断层数值模拟分析方法[J]. 岩土力学， 2006， 27（6）：880-884.

XIAO Ming， CHEN Juntao. Numerical simulation and analysis method of complex geologic faults in large underground chamber[J]. Rock & Soil Mechanics， 2006， 27（6）：880-884. （In Chinese）

[4] 雷显权， 陈运平. 应用非连续模型研究断层对地壳应力的影响[J]. 中南大学学报：自然科学版， 2011， 42（8）：2379-2386.

LEI Xianquan， CHEN Yunping. Effects of faults on crustal stress using discontinuous model[J]. Journal of Central South University： Science and Technology， 2011， 42（8）： 2379-2386.（In Chinese）

[5] GOODMAN R E， TAYLOR R L， BREKKE T L. A model for the mechanics of jointed rock[J]. Journal of Soil Mechanics & Foundations Div， 2015， 94：637-660.

[6] 方濤， 柴军瑞， 胡海浪， 等. Monte Carlo 方法在岩体裂隙结构面模拟中的应用[J]. 露天采矿技术， 2007（1）： 7-9.

FANG Tao， CHAI Junrui， HU Hailang， et al. Application of Monte Carlo method to simulating the fracture interconnectivity of rock mass[J]. Opencast Mining Technology， 2007（1）：7-9.（In Chinese）

[7] 张雨霆， 肖明， 左双英. 基于单元重构的岩土工程复杂地质断层建模方法[J]. 岩石力学与工程学报， 2009， 28（9）：1848-1855.

ZHANG Yuting， XIAO Ming， ZUO Shuangying. Methodology for modeling of complex geological faults in geotechnical engineering based on element reconstruction[J]. Chinese Journal of Rock Mechanics & Engineering， 2009， 28（9）：1848-1855.（In Chinese）

[8] 张雨霆， 肖明， 丁秀丽，等. 复杂岩石块体识别的单元重构-聚合方法[J]. 岩石力学与工程学报， 2012， 31（3）：507-523.

ZHANG Yuting， XIAO Ming， DING Xiuli， et al. Identification method of complex rock blocks using element reconstruction and aggregation technique[J]. Chinese Journal of Rock Mechanics & Engineering， 2012， 31（3）：507-523. （In Chinese）

[9] 朱伯芳. 有限单元法原理与应用[M]. 2版. 北京： 中国水利水电出版社， 1998：112-124.

ZHU Bofang. Principles and applications of finite element method [M]. 2nd ed. Beijing： China Water & Power Press， 1998：112-124.（In Chinese）

[10]龙述尧， 侯淑娟. 关于无网格方法中点插值形函数的研究[J]. 湖南大学学报：自然科学版， 2004， 31（3）：94-97.

LONG Shuyao， HOU Shujuan. A study on shape functions in the meshless point interpolation method[J]. Journal of Hunan University： Natural Sciences， 2004， 31（3）：94-97. （In Chinese）

[11]卞康， 肖明. 退化单元位移场插值计算方法研究[J]. 武汉大学学报：工学版， 2008， 41（6）：16-20.

BIAN Kang， XIAO Ming. Research on calculation methods for displacement field interpolation of degenerate elements[J]. Engineering Journal of Wuhan University， 2008， 41（6）：16-20.（In Chinese）

[12]李新强， 杨松青， 汪小刚. 岩体随机结构面三维网络的生成和可视化技术[J]. 岩石力学与工程学报， 2007， 26（12）：2564-2569.

LI Xinqiang， YANG Songqing， WANG Xiaogang. Generation and visualization technologies of three-dimensional network of rockmass stochastic structural plane[J]. Chinese Journal of Rock Mechanics & Engineering， 2007， 26（12）：2564-2569.（In Chinese）

[13]SHYU K. Nodal-based discontinuous deformation analysis[D]. Berkeley， USA： Department of Civil Engineering， University of California， Berkeley， 1993：10-20.

[14]TIAN Q， ZHAO Z， BAO H. Block fracturing analysis using nodal-based discontinuous deformation analysis with the double minimization procedure[J]. International Journal for Numerical & Analytical Methods in Geomechanics， 2013， 81（4）：409-415.

[15]曹代勇， 吳李泉， 占文锋，等. 断层应力封闭性的FLAC模拟研究[J]. 煤田地质与勘探， 2005， 33（1）：1-3.

CAO Daiyong， WU Liquan， ZHAN Wenfeng， et al. Application of FLAC in simulation of fault stress seal[J]. Coal Geology & Exploration， 2005， 33（1）：1-3. （In Chinese）

[16]翟晓荣， 吴基文， 彭涛，等. 断层对煤层开采应力分布影响的数值模拟分析[J]. 煤矿开采， 2013， 18（5）：14-16.

ZHAI Xiaorong， WU Jiwen， PENG Tao， et al. Numerical simulation for influence of faults on mining stress distribution[J]. Coal Mining Technology， 2013， 18（5）：14-16. （In Chinese）

[17]黄生文， 司铁汉， 陈文胜，等. 断层对大跨度隧道围岩应力影响的有限元分析[J]. 岩石力学与工程学报， 2006， 25（z2）：3788-3793.

HUANG Shengwen， SI Tiehan， CHEN Wensheng， et al. Finite element analyses of influence of fault on large-span tunnel surrounding rock stress[J]. Chinese Journal of Rock Mechanics & Engineering， 2006， 25（z2）：3788-3793. （In Chinese）