大型地下洞室复杂地质断层数值模拟分析方法

第 27卷第 6期 岩 土 力 学 Vol.27 No.6 2006年 6月 Rock and Soil Mechanics Jun. 2006 收稿日期：2004-12-08 基金项目：湖北省环境岩土工程重点实验室开放课题支持课题（No. T110304）。 作者简介：肖明，男，教授，主要从事地下工程结构稳定和水力发电结构研究。Email: xm2003@whu.edu.cn 文章编号：1000－7598－(2006) 06―0880―05 大型地下洞室复杂地质断层数值模拟分析方法 肖 明，陈俊涛 （武汉大学 水资源与水电工程国家重点实验室，武汉 430072） 摘 要：根据大型地下洞室中复杂断层对围岩稳定的影响，提出了隐含断层单元的数值模拟计算方法。该方法将复杂断层单 元隐含在岩体单元中，使得岩体单元的剖分不受断层切割的影响。使用这种隐含的断层单元，不但将受复杂断层切割的大型 复杂地下洞室群的有限单元离散简单化，加快分析计算周期，而且能够有效地反映断层对复杂地下洞室围岩稳定的影响。通 过对三峡地下厂房复杂断层结构的分析计算，证明该方法为复杂断层结构的分析提供了一种十分有效的计算方法和思路。 关 键 词：复杂地质；隐含断层；数值分析 中图分类号：TU 451 文献标识码：A Numerical simulation and analysis method of complex geologic faults in large underground chamber XIAO Ming, CHEN Jun-tao (State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China) Abstract: In order to reflect the influence of complex geologic faults in large underground chamber, the authors put forward a numerical simulation and analysis method for the connotative fault element. The method establishes an equivalent element comprising both rock material and fault material, which can decrease the influence of faults incision. By this way, not only the mesh generation of large underground chambers becomes relatively convenient and feasible, but also the influence of complex geologic faults on stability of surrounding rocks is reflected effectively. At last, a sample about underground chamber of Three Gorges Project is given, which testifies that the method is reasonable and provides an effective computation method and train of thought for the analysis of complex geologic faults structure. Key words: complex geology; connotative fault; numerical analysis 1 前 言 随着国民经济建设的发展和水电建设规模的扩 大，出现了一大批正在规划设计和建设的大型和特 大型地下水电站。例如溪洛渡、小湾、龙滩、水布 垭、西龙池、向家坝、张河湾等水电站。受地形和 地质条件的限制，这些水电站所处的地质环境非常 复杂，地下洞室群被众多的断层所切割。大型地下 洞室通常都是洞室纵横交错，遇到复杂断层与洞室 结构相交后，由于断层与地下洞室相交的任意性， 很难正确模拟断层对洞室围岩稳定的影响。在数值 分析中目前都是将断层进行一定的简化后，用长而 薄的有限单元模拟断层，然后用有限元方法进行分 析计算。但这种方法很难直接模拟断层的厚度、倾 角、走向对复杂地下洞室群的影响，而且长而薄的 有限单元数学模型随着长宽比的加大，其误差越来 越大，因此，对复杂的断层影响只能进行一些定性 分析。本文采用隐含断层单元的数值模拟计算方法 将复杂断层单元隐含在岩体单元中，使得岩体单元 的剖分不受断层切割的影响。这种方法不但使受复 杂断层切割的大型复杂地下洞室群的有限单元离散 简单化，而且能够有效地反应断层对复杂地下洞室 围岩稳定的影响，为复杂断层结构的分析提供了一 种十分有效的计算方法和思路。 2 隐含断层的力学参数确定 假定断层穿过如图 1所示的岩体单元，将该单 元变成一个垂直和平行于断层面的各向异性岩体单 元（图 2），可以用一个均质各向异性的岩体当量 单元模拟断层的结构效应。如果有多条断层穿过岩 第 6期 肖 明等：大型地下洞室复杂地质断层数值模拟分析方法 图 1 岩体单元 Fig.1 Rock element 图 2 当量单元 Fig.2 Equivalent element 体单元，该单元可视为由多条岩层和断层相间组成 的复合单元，设第 j层的岩体厚度为 Hj ，则该层的 层厚系数为： /j jH Hφ = ，且 1jφ =∑ (1) 将断层穿过的复合岩体单元用当量单元替换， 该单元应该满足下列几个条件： (1) 当量单元的应力 iσ 和应变 iε（i =1，2，…,6） 与实际的岩体单元的应力 iσ 和应变 iε 满足 1 ( ) di i jV vV σ σ= ∫∫∫ ， 1 ( ) di i jV vVε ε= ∫∫∫ (2) 式中 V为岩体单元的体积，i =1，2，…,6，( ) , ( )i j i jσ ε 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示第 j层岩体应力和应变。 (2) 当量单元的变形能Ue和实际岩体单元的变 形能 Ur相等。即： r eU U= (3) 1r ( )d 2 x x y y xy xyV U vσ ε σ ε τ γ= + + +∫∫∫ " 1e ( )d 2 x x y y xy xyV U y vσ ε σ ε τ= + + +∫∫∫ " (3) 当量单元和实际岩体单元在断层层面边界 上的应力连续和变形协调。即： 应力连续： ( )x j x xjtσ σ= + ，( )y j x yjtσ σ= + ，( )z j zσ σ= ， ( )xy j xy xyjtτ τ= + ，( )yz j yzτ τ= ，( )zx j zxτ τ= (4) 应变协调： ( )x j xε ε= ， ( )y j yε ε= , ( )z j z zjeε ε= + , ( )xy j xyγ γ= , ( )yz j yz yzjeγ γ= + ， ( )zx j zx zxjeγ γ= + (5) 其中 , xj zjt e 表示在第 j 层断层面上的附加应力和 应变。 将式（4）和（5）代入式（2）积分后，进行整 理可得到： 0j xjtφ =∑ ， 0j yjtφ =∑ ， 0j xyjtφ =∑ ， 0j zjeφ =∑ , 0j xzjeφ =∑ ， 0j yzjeφ =∑ (6) 假定岩体单元本身为各向异性岩体，第 j 层的 平行向弹模和泊松比为 1 1, j jE µ ，垂直向弹模和泊松 比为 2 2, j jE µ ，根据层状各向异性岩体的应力和应 变关系[1,. 2]，依据式（4）和（5），则可用当量单 元应力表示出岩体单元第 j层的应变值： 1 2 1( ) [( ) ( ) ] /x j x x xj j y yj j z jt t Eε ε σ µ σ µ σ= = + − + − 1 2 1( ) [( ) ( ) ] /y j y y yj j x xj j z jt t Eε ε σ µ σ µ σ= = + − + − 2 1 1( ) / [ ] /z j z zj z j j x xj y yj je E t t Eε ε σ µ σ σ= + = − + + + 1( ) ( ) /xy j xy xy xyj jt Gγ γ τ= = + 2( ) /yz j yz yzj yz je Gγ γ τ= + = (7) 2( ) /zx j zx zxj zx je Gγ γ τ= + = 根据能量等效，将式（7）代入式（3），可求 出附加应力 , , xj yj xyjt t t ，再代入式（7）进行整理， 与层状各向异性岩体的应力和应变关系式相比较， 则可得出当量单元平行于断层面的弹模和泊松比 1 1, E µ ，垂直于断层面的弹模和泊松比 2 2, E µ ，即： 2 2 1 1 1 1 1 22 2 11 1 11 1 j j j j j j j jj j E E E E φ µ φ φ µµ µ ⎡ ⎤⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎢= − ⎥⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟⎜ ⎟ ⎜ ⎟⎢ −− − ⎥⎝ ⎠⎝ ⎠ ⎝ ⎠ ⎦⎣ ∑ ∑ ∑ 881 岩 土 力 学 2006年 2 2 1 1 1 1 1 1 2 22 1 11 1 11 j j j j j j j j j jj E E Eµ φ φ µ φµ µ µµ ⎡ ⎤⎛ ⎞ ⎛ ⎞ ⎛ ⎞⎢ ⎥= −⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟⎜ ⎟ − −⎢ ⎥− ⎝ ⎠ ⎝ ⎠⎝ ⎠ ⎣ ⎦ ∑ ∑ ∑ 2 2 2 2 2 2 1 1 1 1 2 211 (1 ) (1 ) j j j j j j j j E E E E µ µφ µ µ ⎡ ⎤⎛ ⎞= − +⎢ ⎥⎜ ⎟⎜ ⎟− −⎢ ⎥⎝ ⎠⎣ ⎦ ∑ 2 1 1 1 1 2 22 2 2 11 1 1 11 1 1 j j j j j j j j j jj j j E E Eµ φ φ µ φ φµ µµ µ µ ⎛ ⎞⎛ ⎞= −⎜ ⎟⎜ ⎟⎜ ⎟⎜ ⎟ −− − −⎝ ⎠⎝ ⎠ ∑ ∑ ∑ (8) 按照上式计算当量单元的等效参数，采用各向 异性有限元进行结构分析，则可计算出当量岩体的 应力。 3 隐含断层应力计算和非线性分析 假定断层的走向为α′，倾角为 β ′（向上倾为 正），有限元整体坐标轴 X，Y，Z的 Y轴方向角为 ϕ，设断层面的局部坐标 Z ′为断层面法向矢量，Y ′ 为断层面走向， X ′垂直于Y ′ Z ′平面。按下式则可 求出断层面局部坐标 X ′ Y ′ Z ′的方向余弦： Z ′轴： 3 cos sin( )l β α ϕ= − ， 3 cos cos( )m β α ϕ= − ， 3 sinn β= Y ′轴： 2 cos( )l ϕ α= − ， 2 cos(90 )m α ϕ= + − ， (9) 2 cos90 0n = = X ′轴： 1 2 3l m n= ， 1 3 2m n l= − ， 1 2 3 3 2n l m l m= − 式中 90α α′= ± （断层倾向南取正号，倾向北取负 号），法向矢量的倾角： 90β β ′= − 。 根据局部坐标 X ′ Y ′ Z ′的方向余弦式（9）和应 力转换关系[3]，可以求出断层面坐标的应力转换矩 阵[S]，根据有限元计算的当量单元应力{ }σ ，按下 式可计算沿断层面的法向应力 zσ ′ 和平行断层面的 剪应力 ,yz zxτ τ′ ′ 分量： { } [ ] { }Tn Sσ σ= (10) 当垂直于断层面的正应力 zσ ′ 大于断层的抗拉 强度 R1时，即： z lRσ ′ > (11) 说明断层出现拉裂破坏，此时必须将超出断层面抗 拉强度的应力 n nz z lRσ σ∆ = − 转化成节点荷载，并将 其转移。当平行于断层面的阻滑力小于滑动力，即： n n 2 n 2( ) /z yz zxf c Kσ τ τ+ + < (12) 式中 , , lf c R 为沿断层面的摩擦系数，凝聚力和抗 拉强度，K为沿断层面滑动的安全系数。此时应将 超出抗剪强度的剪力： n 2 n 2 n( )yz zx zf cτ τ τ σ∆ = + − + (13) 转化成节点荷载，并将其转移后进入下一步迭代计 算。迭代计算方法采用三维弹塑性损伤有限元迭代 计算[4]。 4 工程实例分析 本文针对三峡地下厂房断层结构，采用隐含地 质断层方法对地下厂房洞室开挖的围岩稳定进行了 分析。三峡右岸地下厂房装机 6×70 万 kW，地下 厂房由主厂房、出线洞，尾水调压井、引水洞、尾 水洞等洞室群组成，厂房洞室尺寸为 328 m×30 m ×88 m(长×宽×高)，地下厂房的规模较大，开挖 于闪长斜长花岗岩之中，岩体弹性模量为 30 GPa， 泊松比为 0.2，凝聚力 1.7 MPa，摩擦系数为 1.5。 有 6 条断层穿过地下厂房，断层参数为：F20 产状 为 250º∠70º，厚度 1.5 m；F32产状为 250º∠69º， 厚度 0.4 m；F22产状为 250º∠70º，厚度 0.6 m； F84 产状为 350º∠60º，厚度 1.0 m；F100 产状为 354º∠84º，厚度 0.4 m；F205产状为 10º∠78º，厚 度 0.6 m。三维有限元分析计算网格包括了主厂房洞 室、出线洞、尾水洞、引水洞和尾水调压井等洞室， 共剖分了 8 860个 20节点空间等参单元（见图 3）。 图 3 三维有限元计算网格剖分图 Fig.3 3D FEM mesh 对地下洞室进行了毛洞无支护和锚固支护两种 不同工况的分析计算，通过分析比较可看出以下特 征： (1) 在无支护条件下，主厂房洞周的破坏区较 大，开挖完毕围岩的总破坏体积达到57.94×104 m3， 塑性耗散能达到 70.9×105 kN·m。说明在无支护条 件下，断层对围岩稳定影响较大。采用锚固支护后， 整个洞周的塑性开裂破坏区都发展较慢，围岩的总 破坏体积减到 10.2×104 m3，塑性耗散能减到 4.02 882 第 6期 肖 明等：大型地下洞室复杂地质断层数值模拟分析方法 ×105 kN·m，主厂房顶拱的破坏范围很小，基本都 集中在断层出入部位，边墙的破坏区基本限制在 3 m范围内，只是在断层穿过的部位出现了局部开裂 区（见图 4）。说明采用隐含断层网格，能够充分 反应断层对围岩稳定的影响。 图 4 洞室开挖完毕后主厂房洞周破坏区分布 Fig.4 Failure zone distribution (2) 从洞周的应力分布规律看，由于断层的切 割，使得在断层处的围岩应力有较大突变，特别是 在断层与洞室交叉口处的应力集中较为明显（见图 5）。在无支护条件下，断层引起洞室交叉口处的最 大压应力达到 88.2 MPa，采用锚固支护后，应力集 中程度比无支护有较大减小，但在断层与洞室交叉 口处的最大压应力仍然达到 59.5 MPa。说明采用隐 含断层较好地反应了地质断层对围岩应力分布的影 响。 图 5 洞室开挖完毕主厂房第一主应力矢量 Fig.5 First principal stress vector (3) 从断层的正应力分布规律看，垂直于断层 的正应力基本为压应力，其量值一般分布在 －0.46～－9.0 MPa，在断层与洞室的交口处压应力 值较小，只有－0.10 MPa，在引水管与断层交口处 局部出现了 0.08 MPa的拉应力，说明在该处断层可 能张开，对该处的断层支护应适量加强。位于厂房 边墙的断层正应力值较大，其值分布在－2.24～－ 9.03 MPa；位于厂房顶拱的断层正应力值较小，其 值分布在－0.1～－0.84 MPa之间（见图 6）。由此 看出，顶拱断层受力条件较差。 图 6 洞室开挖完毕主厂房断层正应力大小 Fig.6 Normal stress of faults (4) 从沿断层的滑动安全系数分布规律看，在 断层与洞室的交叉口处，滑动安全系数一般都<1.0， 其值分布在 0.3～1.0之间，其他部位的安全系数都 在 2.0～16之间（见图 7），厂房边墙断层的安全系 数值都>1.0，厂房顶拱的断层安全系数较小，其值 在 0.3～1.3之间(见图 7)。说明顶拱断层处的安全稳 定条件较差，对顶拱断层采用一些随机锚索支护加 固是必要的。 图 7 洞室开挖完毕主厂房断层滑动系数 Fig.7 Slide-coeffcient of faults (5) 洞室开挖后，顶拱的位移为 0.3～0.47 cm， 但在断层处位移达到 0.68 cm，开挖完毕后主厂房上 下游边墙的最大位移为 1.88 cm和 1.54 cm，但在断 层处达到 2.68 cm（见图 8）。说明断层对洞周位移 分布有较大影响。 (6) 主厂房顶拱的锚杆应力一般为 13～20 MPa，但在断层处锚杆应力达到 24～25 MPa；洞室 边墙的锚杆应力一般在 20～67 MPa，但在断层处锚 杆应力达到 90～185 MPa（见图 8）。说明断层使 洞周锚杆应力明显加大，加强对断层的支护是必要 的。 883 岩 土 力 学 2006年 图 8 洞室开挖完毕主厂房位移和锚杆应力 Fig.8 Displacement of workshop and stress of anchor bolts 综上分析可以看出，断层对围岩的破坏区、洞 周应力、锚杆应力和洞周位移分布都有较大影响， 采用隐含断层方法能够充分反映断层对地下洞室围 岩稳定的影响，为对断层的支护处理提供了依据。 5 结 论 采用隐含断层单元模拟复杂的断层结构，使得 复杂纵横交错的大型地下洞室群有限单元网格剖分， 不受复杂断层任意切割的影响，这种方法不但使复 杂的地下洞室群建模简单快速化，而且能够有效地 反应断层对复杂地下洞室围岩稳定的影响，缩短计 算周期，为复杂断层结构的分析提供了一种十分有 效的计算方法和思路。 参 考 文 献 [1] 朱伯芳. 有限单元法原理与应用[M]. 2 版, 北京: 中国 水利电力出版社, 1998. [2] 王勖成, 邵敏. 有限元基本方法原理与数值方法[M]. 北京: 清华大学出版社, 1991. [3] 杨桂通. 弹性力学[M]. 北京. 人民教育出版社, 1980. [4] 肖明. 西龙池抽水蓄能电站地下厂房围岩稳定三维非 线性分析[J]. 岩石力学和工程学报, 2000, 19(5): 557 －561. XIAO Ming. Three dimensional stability analysis of surrounding rock of underground houses of Xilongchi Pumped Storage Hydropower Station[J]. Chinese Journal of Rock Mechanics and Engineering, 2000, 19(5): 557－561. 上接第 856页 [5] Chijimatsu M, Fujita T, Kobayashi A. and Nakano M. Experiment and validation of numerical simulation of coupled thermal, hydraulic and mechanical behavior in engineering buffer material[J]. Int. J. for Numerical and Analytical Methods in Geomechanics, 200, 24(4): 403 －424. [6] Chijimatsu, M., Fujita, T., Sugita, Y., Amemiya, K. and Kobayashi A. Field experiment, results and THM behavior in the Kamishi Mine experiment[J]. Int. J. of Rock Mechanics and Mining Sciences, 2001, 38(1): 79 －94. [7] Leiws R W, Schrefler B A. The finite element method in the deformation and consolidation of porous media[M]. New York: John Wiley, 1987. [8] 张玉军. 热—水—应力耦合弹塑性二维有限元程序的 开发与应用尝试[J]. 岩土力学, 2005, 26(2): 169－174. ZHANG Yu-jun. Development and use try of 2D FEM code for coupled thermo-hydro-mechanical elstoplastic analysis[J]. Rock and Soil Mechanics, 2005, 26(2): 169 －174. [9] 张玉军. 一种模拟热—水—应力耦合作用的节理单元 及其数值分析[J]. 岩土工程学报 , 2005, 27(3): 270 －274. ZHANG Yu-jun. A kind of joint element simulating coupled thermo-hydro-mechanical phenomenon and relevant numerical analysis[J]. Chinese Journal of Geotechnical Engineering, 2005, 27(3): 270－274. 884