含褶皱巷道围岩破裂碎胀大变形机理FDEM 数值模拟研究

贾永杰

（山西工程职业学院 资源与安全工程系，山西 太原 030009）

0 引 言

当前对含褶皱巷道软弱围岩破裂碎胀大变形机理的研究还存在诸多不足之处，如大都采用连续性数值模拟方法，该方法仅能模拟巷道开挖后围岩的弹塑性连续变形，无法模拟围岩破裂过程、破碎块体间的接触挤压和块体的宏观大运动过程；此外，未考虑地应力侧压系数、褶皱几何形态和岩层特性的影响。因此，本文将采用有限元-离散元耦合数值模拟方法（FDEM） 模拟含褶皱巷道围岩破裂碎胀大变形机制，并研究地应力侧压系数、褶皱几何形态和不同岩层特性的影响。相比于传统的单一连续性或非连续性数值方法，FDEM 能够模拟完整岩石材料受载后的弹塑性连续变形及开裂失效过程，也能够模拟破裂块体间的接触挤压作用及块体的宏观大运动过程，且裂隙网络数目不受限、计算效率较高，模拟结果正确性已经过室内试验和工程现场监测所验证。

1 FDEM 数值模拟基本原理

1.1 FDEM 本构模型

FDEM 基本原理：将岩石材料划分为三角形单元和初始无厚度的四边形节理单元，三角形单元在外荷载或其他荷载作用下仅发生线弹性变形，四边形节理单元则经历峰前线弹性变形→峰后塑性损伤→断裂失效全过程，失效后的节理单元不再进入计算程序并产生一条微裂隙，且两侧的三角形单元由粘结关系转换为接触关系。由于FDEM 基本原理能在诸多文献中轻易找到，本文不再赘述，仅介绍改进后的四边形节理单元本构模型。

本文对四边形节理单元粘结应力本构模型进行了修正，如下：

式中：σn、τ 分别为法向应力和切向应力；o、分别为拉伸/压缩位移和剪切位移；op、ot、sp和st分别为拉伸峰值位移和拉伸极限位移、剪切峰值位移和剪切极限位移；ot和st的值采用Deng 等[18]提出的相对值形式，即相对于网格尺寸，而非绝对位移形式，如此可消除I 型和II 型断裂能GI、GII 对网格尺寸的依赖性，即不同的网格尺寸可采用相同的断裂能；ft、c分别为抗拉强度和粘聚力；z为峰后软化函数。

1.2 层状岩体FDEM 数值模拟原理

可采用断续裂隙网络法（DNF） 和弥散方法模拟（smeared method），而DFN 难以模拟真实层状岩体特性[9]，因而在后续的研究中更多地采用弥散方法[13]，其基本原理概述如下：将层状岩体的层理面通过直接建模的方式表征（图1a），平行于层理面的四边形节理单元岩石强度最低（包括抗拉强度ft、粘聚力c、I 型断裂能GI 和II 型断裂能GII）、垂直于层理面的四边形节理单元岩石强度最高（图1b），与层理面呈其他角度的四边形节理单元强度采用线性计算公式：

图1 层状岩体FDEM 数值模拟原理Fig.1 FDEM numerical simulation principle of layered rock mass

式中：ft,γ、cγ、GI,γ、GII,γ为与层理面呈γ 角度的四边形节理单元抗拉强度、粘聚力、I 型断裂能和II型断裂能；ft,min、cmin、GI,min、GII,min为平行于层理面的四边形节理单元抗拉强度、粘聚力、I 型断裂能和II 型断裂能；ft,max、cmax、GI,max、GII,max为垂直于层理面的四边形节理单元相应参数。

实际上，弥散方法是一种等效方法，即真实岩体若干层理面内的岩石被表征为FDEM 中的两相邻层理面内的网格，这些若干层内的真实岩石宏观力学特性可通过FDEM 四边形节理单元强度参数的线性变化来反映，避免在FDEM 中划分非常微小的层理间距和网格尺寸，提高了建模效率和计算效率。

2 含褶皱岩体破裂碎胀大变形机理FDEM数值模拟

2.1 模型与参数

平直岩层经水平地应力挤压作用后将发生弯曲变形，Fereshtenejad 等[1]指出可采用正弦/余弦函数表示简单的褶皱构造。

为了消除异形断面的影响，本文以TBM 掘进圆形巷道为例。如图2（a） 所示，在本文中，将波长恒定为3.6 m，不同形状的褶皱通过式（4）中的a值体现。将a值分别设定为0（平直层理面）、0.5 m 和1.5 m，建模过程如下：根据公式（4） 获取图2（a） 所示9 个绿色拐点的坐标，在Gmsh 软件中采用样条曲线Spline 命令将这9 个点相连获得一条余弦曲线，同一层面的余弦曲线根据相同的方法绘制，首尾相连获得一条褶皱曲线，将该层的褶皱曲线向上、下平移，确保翼部间距t=1.0 m，由此形成整个褶皱构造数值模型。

巷道开挖位置如图2（a） 所示，圆形巷道直径3.0 m，模型为边长80 m 的正方形，以a=1.5 m为例，其模型如图3（b） 所示，为减小网格数目、提高计算效率，将模型划分为3 个区域[13]：远场区、网格细化区和巷道区（亦称为核心材料区）。褶皱构造仅在巷道区和网格细化区内体现，远场区内不设置褶皱构造，但需确保所有的裂隙网络在网格细化区内扩展。本文中网格细化区为边长约30 m 的正方形，巷道周边网格尺寸h=0.1 m。巷道开挖模拟分为两个阶段：地应力施加阶段和开挖模拟阶段。在地应力施加阶段，根据所需地应力计算出所有节点的相应节点力，而后将节点力反向施加至对应节点上，在该阶段不插入四边形节理单元，仅采用三角形单元，且模型边界自由，系统将产生巨大动能，在粘滞阻尼和临界迟滞阻尼[10]作用下，模型达到平衡，获得所需地应力；在开挖模拟阶段，插入初始零厚度的四边形节理单元，固定模型边界以保持第一阶段获得的地应力；由于四边形节理单元的插入，三角形单元间会发生非常微小的嵌入，模型再次产生较大的动能，但在阻尼作用下很快达到平衡状态，随后可激活巷道开挖模拟程序。

图3 真实岩层褶皱构造及FDEM 数值模型Fig.3 Real rock fold structure and FDEM numerical model

输入表1 参数，对于层状岩体可直接取表1 中的参数；对于各向同性岩体，取表1 中垂直层理面的参数值，该参数的可靠性已在文献[13]中得到了验证。

表1 FDEM 模拟参数[13]Table 1 Parameters of FDEM simulation

2.2 各向同性围岩

在各向同性岩体中，各层岩体物理力学特性相同，且层理面力学参数与完整岩石力学参数相同。将地应力侧压系数λ 分别设定为0.5（水平地应力σh=14 MPa、垂直地应力σv=28 MPa） 和2.0（σh=28 MPa、σv=14 MPa），得到如图4～图5 所示的模拟结果。

根据模拟结果可知，在不同地应力侧压系数和褶皱形态下，巷道均因围岩的破裂碎胀产生严重大变形灾害，其机理可表述为：随着巷道核心材料的逐步软化（在真实巷道掘进中体现为随着掌子面的推进），巷道表面围岩产生径向应力降低、切向应力升高现象，即为卸围压、升轴压的三轴压缩力学模型，当升高的切向应力超过该围压下的岩体强度时，围岩发生X 型共轭剪切破裂并伴随着拉伸断裂，最大切向集中应力向深处的完整岩体转移，直至开挖完毕后切向应力在围岩深处与岩体强度达到极限平衡状态，裂纹不再往更深处扩展；破裂后的块体一方面发生弹性变形恢复，产生体积膨胀，另一方面主剪切带发生滑移剪胀效应、块体间不啮合产生大量空隙，更重要的是破碎块体间由于运动速度的不一致性产生相互挤压导致块体发生翻转大运动并向巷道空间内移动，使巷道断面积急剧缩小，上述过程为围岩体的破裂碎胀大变形灾变过程。

以图4（a） 为例，由于巷道开挖导致径向应力降低、切向应力升高，进而产生的X 型共轭剪切裂隙与牛双建[11]和潘一山[3]等人通过室内巷道开挖模型试验得到的围岩破裂形态和理论剪切滑移线场（图6a） 非常相似；此外，巷道表面破碎块体的翻转大运动形态与潘一山[12]等人通过室内模型试验结果非常相似，如图5（b） 所示。在围岩位移场方面，本文通过FDEM 数值模拟得到的围岩位移场（图6b） 与洛锋等人基于室内模型试验和FLAC3D 数值模拟提出的解析解非常相似（图6a），亦可分为3 个区域：第I 破裂带、第II 破裂带和第III 破裂带，且第I 破裂带呈楔形三角形形态。

图6 巷道围岩X 型共轭剪切滑移解析解及FDEM 数值模拟结果Fig.6 Analytical solution of X-type conjugate shear slip of roadway surrounding rock and FDEM numerical simulation results

根据本节模拟结果可知，对于各向同性岩层且层理面参数与围岩体一致时，褶皱的存在对围岩破裂碎胀大变形影响不大，均呈现以X 型共轭剪切滑移破裂为主且伴随少量拉伸裂隙的对称性破裂碎胀大变形。当然，对于不同的褶皱形态，围岩碎裂形态和裂隙网络形态存在微小差异，这是由于FDEM 数值模拟对网格的依赖性造成的。

3 结 语

本文采用有限元- 离散元耦合数值模拟方法FDEM 揭示了不同地应力侧压系数和不同褶皱形态下的围岩破裂特征和破裂碎胀大变形机理。在不考虑褶皱形成后残余构造应力的影响下，对于各向同性围岩体而言，褶皱的影响较微弱。巷道开挖卸荷引起的切向集中应力超过岩体强度时发生X 型共轭剪切破坏并伴随少量拉伸裂隙，破碎岩块沿着主剪切裂隙带发生滑移剪胀，此外，大量岩块在相互接触挤压下产生向巷道空间的宏观大运动（包括平移和翻转大运动） 并产生大量空隙，造成破碎岩体的体积膨胀并引起巷道断面积的急剧缩小，是破裂碎胀挤压大变形致灾机理。