围岩隧道蠕变塌方的离散元模拟研究

Discrete Element Simulation of Tunnel Creep Collapse in Surrounding Rock

摘要: 本文以某公路隧道为工程背景，采用离散元数值模拟方法，系统探讨了开挖扰动和蠕变效应对裂隙围岩稳定性的影响，揭示了塌方灾害的时空演化规律。通过PFC3D软件建立了隧道蠕变模型，模拟了隧道围岩在开挖扰动和长期蠕变作用下的应力分布和变形特征。研究结果表明，隧道开挖后围岩应力分布发生显著变化，尤其是在隧道周围的围岩区域，应力集中现象明显。通过应力腐蚀模型，模拟了材料在长期应力和腐蚀环境耦合作用下的力学响应，有效捕捉了岩石在长期载荷作用下的变形与破坏特征。本研究为隧道工程支护设计和长期安全运营提供了科学依据，对实际工程具有重要的参考价值。

Abstract: Based on the engineering background of a highway tunnel, the influence of excavation disturbance and creep effect on the stability of fractured surrounding rock is systematically discussed by using discrete element numerical simulation method, and the spatio-temporal evolution law of collapse disaster is revealed. The creep model of tunnel is established by PFC3D software, and the stress distribution and deformation characteristics of tunnel surrounding rock under excavation disturbance and long-term creep are simulated. The results show that the stress distribution of surrounding rock changes significantly after tunnel excavation, especially in the surrounding rock, the stress concentration is obvious. The stress corrosion model is used to simulate the mechanical response of the material under the coupling action of long-term stress and corrosion environment, and effectively capture the deformation and failure characteristics of the rock under long-term load. This study provides scientific basis for tunnel engineering support design and long-term safe operation, and has important reference value for practical engineering.

1. 引言

隧道工程在交通基础设施建设中扮演着重要角色，但在其施工过程常面临复杂地质条件带来的挑战。常常在开挖隧道过程中会遇到岩爆、渗流、管涌等现象发生[1]-[4]，其中一个困扰地下工程施工的难题就是围岩随时间的变化发生的蠕变现象。裂隙围岩作为隧道工程中常见的地质构造，其不仅受到开挖扰动的影响，还受到长期蠕变效应的控制，其稳定性直接影响隧道的安全性与经济性。近年来，国内外学者针对隧道围岩的蠕变行为和塌方机理开展了大量研究。例如，2015年，马乾天等[5]基于颗粒流理论(PFC)对块石胶结充填体进行了蠕变模拟研究，分析了灰砂比对块石胶结充填体蠕变破坏过程的影响；张学朋等[6]通过PFC程序建立岩石二维离散元虚拟单轴蠕变试验模型，研究了颗粒粒径、摩擦因数和颗粒刚度比对岩石蠕变特性的影响，得出了颗粒流Burgers模型能够有效模拟岩石蠕变行为并为岩石流变的细观机制研究提供有益参考的结论；2018年，潘卫东等[7]基于离散元方法(PFC3D)和伯格斯模型，对深部巷道围岩的蠕变特性进行了数值模拟研究，通过虚拟蠕变试验标定了粉砂岩的微观参数，并与室内试验数据对比验证了离散元方法的有效性，得出了该方法能够较好地模拟岩石蠕变行为并为深部巷道蠕变问题研究提供参考的结论。2020年，王瑶等[8]基于应力腐蚀理论和离散元法(PFC 5.0)，通过编写考虑颗粒破碎的粗粒土蠕变程序，分析了接触模型和应力腐蚀模型细观参数对蠕变的影响规律及其敏感性，并基于室内试验结果进行了参数标定，得出了接触模型中颗粒摩擦系数和胶结键粘结强度以及应力腐蚀模型中激活应力、材料参数β₁和β₂对蠕变影响显著的结论。2023年，骆祚森等[9]基于PFC中的平行黏结模型(PBM)，考虑水–岩作用对黏结的弱化和材料特性随时间的变化，提出了水–岩作用下砂岩蠕变的离散元模拟方法，研究了砂岩在水–岩作用下的蠕变损伤机制，得出了水–岩作用导致砂岩微裂纹扩展加速、剪性破坏增强以及胶结能降低的结论；2023年，黄亚军等[10]基于PFC3D离散元模拟和FLAC3D数值分析，研究了软弱夹层蠕变特性对矿山边坡稳定性的影响，分析了软弱夹层倾角、厚度和层数对复合岩样强度和变形的影响，得出了考虑蠕变特性时边坡安全系数降低约10%的结论，为矿山边坡稳定性分析提供了更安全可靠的参考。2024年，李庆文等[11]基于PFC2D和伯格斯–平行黏结混杂接触模型(BPB)，对砂质泥岩在不同加轴压卸围压条件下的蠕变特性进行了细观研究，分析了其应力–应变行为、裂纹演化和能量变化规律，得出了轴–围增卸比对蠕变失稳、裂纹发育和能量演化有显著影响的结论。这些研究为本文提供了重要的理论基础和参考。为研究本文结合某公路隧道的工程背景，采用离散元模拟方法，探讨开挖扰动对裂隙围岩稳定性的影响，揭示塌方规律，为工程支护设计提供参考。2024年，Zhang等[12]提出了一种基于平行键和Kelvin-Voigt接触模型的蠕变模型，对梯度蠕变加载及其产生的裂纹进行全过程的模拟，有效地反映了蠕变过程中裂纹的扩展，并反映了蠕变过程中的最终破坏模式。

为深入研究隧道裂隙围岩在开挖扰动和蠕变效应共同作用下的变形及塌方规律，本文以某公路隧道为工程背景，采用离散元数值模拟方法，系统探讨开挖扰动和蠕变效应对裂隙围岩稳定性的影响，揭示塌方灾害的时空演化规律，为工程支护设计和长期安全运营提供科学依据。选取某公路隧道为研究对象，该隧道围岩以裂隙发育的软岩(如泥岩、砂岩)为主，裂隙密度和方向随机分布。施工过程中出现明显蠕变变形，导致局部塌方。

2. 数值模拟

2.1. 模型建立

选取某公路隧道为研究对象，该隧道围岩以裂隙发育的软岩(如泥岩、砂岩)为主，施工过程中出现明显蠕变变形，导致局部塌方。为研究隧道埋深对隧道自稳的影响，本文采用改变重力方向的固定荷载来模拟不同埋深下隧道的蠕变失稳机理。隧道模型建立使用PFC3D5.0软件，该软件是一款以非连续介质力学模拟为目标，采用离散元法作为基本理论的三维分析程序，特别适用于因不连续界面导致变形和破坏现象的机制性研究。

Figure 1. Discrete element 3D model of tunnel local structure

图1. 隧道正面结构离散元三维模型

本构模型以实际地质参数为基础建立，模型几何尺寸设定为长 × 高 = 20 m × 30 m。如图1所示，隧道采用标准马蹄形断面，断面跨径实测16 m，侧墙垂直高度8 m，拱顶为半径8 m的1/2圆弧结构。根据地应力实测数据，模型顶部施加垂直均布荷载2 Mpa，侧向采用水平压力系数K = 1.0的对称加载模式，即水平向压力与垂直荷载等值。模型底部采用全自由度固定约束，双侧边界设定法向位移约束，允许竖直方向自由位移。

整体建模的流程为：首先建立空间区域，模型的高度和宽度通常为隧道直径的3~5倍，以减小边界效应对计算结果的影响。随后在区域内分布坐标种子以固定生成颗粒的随机点位；在区域内建立刚体墙，调整局部坐标之后塞入特定粒径的颗粒，基于容器球形颗粒填充球形颗粒，此时根据ball property (fish命令)赋予球密度和球体刚度和摩擦系数相关颗粒属性。Cycle (运行)至颗粒间重叠量近乎于0后，删除不在容器覆盖范围内的小球，并对整个容器内存在的实心颗粒进行预先挤压，将颗粒向中心轴施加压力。

2.2. 蠕变模型的建立

Potyondy和Cundall [13] [14]首次提出平行黏结模型(Parallel-Bond Model)，用于模拟岩石颗粒间的胶结作用。该模型的核心假设是，在颗粒接触点处存在一个与两颗粒中心连线垂直的截平面。此截平面为摩擦界面，仅能传递法向力和切向力，但不能承受力矩。同时，平行黏结模型引入了黏结界面，该界面不仅可以传递力和力矩，还能够抵抗颗粒间的相对旋转。

应力腐蚀效应的数值模拟通常基于颗粒离散元方法(如PFC3D)，旨在研究材料在长期应力和腐蚀环境耦合作用下的力学响应。该模型的核心机制是通过动态调整颗粒间平行粘结模型的黏结强度来模拟腐蚀过程。具体实现方式是借助软件内置的FISH语言，将颗粒接触半径(pb_rmul)作为时变参数，使其随蠕变应力和时间的累积逐步减小，从而表征黏结强度的退化规律。

在围岩蠕变力学行为的研究中，这一模型能够有效捕捉围岩在长期荷载作用下的变形与破坏特征[15]。本研究采用应力腐蚀模型构建数值试样，并基于试验数据对模型参数进行优化，以确保模拟结果的可靠性。图2展示了该模型的示意图。

Figure 2. Simulation of contact creep of spherical particles

图2. 球颗粒接触蠕变模拟

PSC裂纹扩展速度模型公式：

$\frac{\text{d}\overline{D}}{\text{d}t}=\{\begin{array}{l}0,\overline{\sigma }<{\overline{\sigma }}_a\\ -{\beta }_{1}e{\beta }_2{}^{\left(\frac{\overline{\sigma }}{{\sigma }_c}\right)},{\overline{\sigma }}_a<\overline{\sigma }<{\overline{\sigma }}_c\\ -\infty ,\overline{\sigma }>{\overline{\sigma }}_c\end{array}$ (1)

式中：$\overline{D}$ 为平行键粘结平均直径；t为时间；$\overline{\sigma }$ 为颗粒间平行键粘结最大拉应力；${\overline{\sigma }}_a$ 为颗粒间应力腐蚀阈值；${\overline{\sigma }}_c$ 为颗粒间抗拉强度阈值；${\beta }_1$ 、${\beta }_2$ 为材料系数；

应力腐蚀模型通过调整颗粒间接触半径(半径乘子)的变化来模拟材料在不同应力条件下的力学行为。其工作原理可分为以下三个阶段：

${\overline{E}}_k=\frac{1}{2}\left(\frac{{\overline{F}}_n^2}{{\overline{k}}_{n}A}+\frac{{\Vert {\overline{F}}_s\Vert }^2}{{\overline{k}}_s\overline{A}}+\frac{{\overline{M}}_t^2}{{\overline{k}}_s\overline{J}}+\frac{{\Vert {\overline{M}}_b\Vert }^2}{{\overline{k}}_n\overline{I}}\right)$ (2)

式中：${\overline{E}}_k$ 为应变能；${\overline{F}}_n$ 、${\overline{F}}_s$ 分别为法向和切向力；${\overline{k}}_n$ 、${\overline{k}}_s$ 分别为法向和切向强度；${\overline{M}}_t$ 、${\overline{M}}_b$ 分别为切向和法向的力矩；$\overline{I}$ 、$\overline{J}$ 分别为接触处惯性矩和极惯性矩；$A$ 为平行粘结键的接触面积(PFC^2D中接触面积为$2\overline{R}t$ ，$t=1$ ；PFC^3D中接触面积为$\pi {\overline{R}}^2$ )。

3. 参数选取和标定

本次模拟的是岩石隧道在恒定荷载下的蠕变变化，模型的瞬态细观参数标定是根据“试错法”，类如不同法向应力状态下的剪应力破坏值。通过建立PFC离散元模型，调节颗粒间粘结强度，形成损伤微裂缝。在PFC中，颗粒与颗粒之间采用平行粘结接触，模拟土颗粒–土颗粒之间的接触，颗粒细观参数见表1。

Table 1. Model parameter table

表1. 模型参数表

参数/单位	数值
颗粒最小半径/mm	4.0
颗粒最大半径/mm	5.0
颗粒半径比值	1.25
颗粒弹性模量/GPa	9
颗粒刚度比	1.5
颗粒摩擦系数	0.5
墙体摩擦系数	0
粘结接触抗拉强度平均值/GPa	4
粘结接触切向强度平均值/GPa	4

4. 模拟结果

根据上述参数建立基于球颗粒之间的微观离散元模型，根据球颗粒的粘结作用分析宏观自重应力及水平应力稳定荷载下，对该隧道的开挖造成的流动影响。图3为隧道内部的球单元颗粒在x方向的速度矢量示意图。受到上部自重应力荷载的加载，隧道围岩沿应力小的分布方向坍塌。图4为隧道开挖初始状态下的应力分布图，应力在整个岩石结构中传递，并在围岩附近显著集中。这种应力集中现象有助于强化围岩的原始结构，提高其整体稳定性。图5展示了隧道开挖完成后的应力分布变化情况。与开挖前的状态相比，开挖后岩体中的应力分布进行了重新调整，尤其是在隧道周围的围岩区域，应力分布发生了显著变化。通过图4的力链图可以看到，红色较粗的接触力链到了开挖后就削弱了很多。这是由于应力随着蠕变时间的增长发生而产生了断裂，导致力的传递形成了局部破碎。将本该传递的荷载因破裂重新分配至围岩其他支撑链上，这种应力的局部化现象使得围岩内部产生微裂缝，并在蠕变过程中扩展为贯穿性裂缝，最终形成宏观断裂带。隧道拱顶区域由于重力作用，应力集中更为显著，导致该区域的颗粒位移速度显著增加。

在隧道开挖的初期，应力的变化是一个动态重新分配的过程，通常可以分为以下几个阶段来描述：

1) 开挖瞬间的应力集中：隧道开挖后，开挖面附近的岩体被移除，导致原有应力平衡被打破。由于开挖面的出现，岩体失去了支撑，应力迅速集中在开挖面附近。这种应力集中现象主要发生在开挖面周围，表现为局部应力显著增加，尤其在开挖面的轮廓线附近，可能形成应力集中带。

2) 应力的重新分配：随着开挖的进行，开挖面附近的应力集中逐渐向周围岩体传递。岩体内部的应力开始重新分布，应力从开挖面附近向深层岩体转移。这一过程中，开挖面附近的应力逐渐减小，而深层岩体中的应力逐渐增加，最终达到一个新的应力平衡状态。

3) 深层岩石的应力分布：在应力重新分配完成后，开挖面的应力集中现象显著减弱，应力分布趋于均匀。深层岩石成为主要承载区域，承担了开挖过程中转移的应力。这种应力分布的变化反映了岩体的自稳能力和应力传递特性，同时也为后续支护设计提供了重要依据。

Figure 3. Vector motion direction of tunnel inner wall particles

图3. 隧道内壁颗粒矢量运动方向

Figure 4. Stress distribution at the beginning of tunnel excavation

图4. 隧道开挖初期的颗粒接触力分布

Figure 5. Particle contact force distribution after tunnel excavation

图5. 隧道开挖后的颗粒接触力分布

图6展示的是通过离散元模型(DEM)模拟的岩石试样在单轴压缩条件下的应力–应变关系曲线。通过对曲线的分析，可以将其划分为以下几个典型阶段。在初始弹性阶段(0.00~0.05)，应力随应变的增加呈线性增长，表明岩石试样主要表现为弹性变形。此阶段反映了岩石在低应力水平下的刚性特征，变形速率较低，且符合胡克定律。进入非线性强化阶段(0.05~0.30)，随着应变的增加，应力增速逐渐放缓，曲线呈现非线性特征。这一阶段反映了岩石内部微裂纹的萌生和扩展，导致应力传递路径发生变化，岩石逐渐从弹性状态向塑性状态过渡。在峰值强度阶段(0.30~0.35)，岩石试样达到峰值应力后，承载能力开始下降，曲线出现拐点。这一阶段标志着岩石内部微裂纹的贯通和宏观破裂面的形成，岩石结构整体失稳。随后的应变软化阶段(0.35~0.50)中，应力随应变的增加而逐渐降低，表明岩石进入软化阶段。此阶段反映了岩石在宏观破裂后的残余强度特性，变形速率加快，岩石结构逐渐破坏。在最后的残余稳定阶段(0.50~0.60)，应力变化趋于平稳，岩石试样的残余强度基本保持不变。此阶段反映了岩石在经历宏观破坏后仍保持一定的承载能力，变形趋于稳定。

Figure 6. The displacement of the tunnel under pressure changes in the direction of gravity stress

图6. 自重应力方向的隧道受压下位移变化

总体而言，该应力–应变曲线揭示了岩石试样在单轴压缩条件下的变形与破坏特征，反映了岩石从弹性变形、微裂纹扩展、峰值强度、应变软化到残余稳定的全过程。这一演化规律为理解岩石的力学行为和强度特性提供了重要依据，对岩土工程的设计与安全评估具有重要的参考价值。

摩擦力是控制围岩隧道塌方机制的核心力学参数之一，其作用贯穿于围岩从弹性变形到塑性破坏的全过程。基于离散元模拟(PFC3D)中不同摩擦系数(μ = 0.3~0.6)的定量分析结果(图7位移响应)。颗粒间的摩擦力f_t由法向力f_n和摩擦系数μ决定：

f_t ≤ μf_n

其中，𝜇 = tan(ϕ)，ϕ为颗粒间内摩擦角。当剪切力超过μf_n时，颗粒发生相对滑移。摩擦力通过与颗粒间的相互作用，延缓了围岩的蠕变失稳过程。通过调整颗粒级配或注浆加固，可以提高μ和f_n，从而增强围岩的长期稳定性。

粘结强度是岩土材料抵抗颗粒间分离或滑移的临界力学参数，当外部应力低于粘结强度时，材料主要表现为弹性变形，变形是可恢复的，且不会随时间显著发展。然而，当外部应力超过粘结强度时，材料内部的黏结结构开始逐渐破坏，导致塑性变形的产生。故在低粘结强度时，图7(b)中0.10 Gpa仅仅出现了蠕变的第一第二阶段，岩石试样并未出现破坏，应变随时间变化逐渐变小。对比0.25 Gpa的粘结强度下，应变率在蠕变第二阶段下显著上升，此时隧道已经出现了明显失稳垮塌的变形。

(a) (b)

Figure 7. Influence of interparticle friction and bond strength on tunnel creep

图7. 颗粒间摩擦力和粘结强度对于隧道蠕变的影响

图8所示是加载应力为2 MPa时隧道开挖剖面的损伤演化过程。损伤导致岩石弹性模量的降低，当单元完全损伤时以黑色表示，图中损伤的演化代表了裂纹的衍生过程。损伤的演化过程本质上是裂纹从微观到宏观的衍生过程，揭示了岩石在应力作用下的破坏机制。微裂纹主要围绕着开挖面的顶部和底板展开，特别是底部的颗粒由于破碎之后呈现向上隆起的状态，反映了应力集中区域的破坏特征。通过分析损伤演化过程，可以更好地理解隧道开挖后的力学行为，并采取针对性的工程措施，确保施工安全和结构稳定性。

5. 结论

1) 通过PFC3D建立了演示隧道蠕变模型，数值模拟结果表明，隧道开挖后围岩应力分布发生了显著变化，尤其是在隧道周围的围岩区域，应力集中现象明显。这种应力集中现象有助于强化围岩的原始结构，提高其整体稳定性。

Figure 8. Tunnel creep time compression microcrack distribution map

图8. 隧道蠕变时间受压微裂纹分布图

2) 隧道围岩在开挖扰动和蠕变效应共同作用下，表现出明显的变形和塌方规律。受到上部自重应力荷载的加载，隧道围岩沿应力小的分布方向坍塌，揭示了塌方灾害的演化规律。

3) 通过应力腐蚀模型，模拟了材料在长期应力和腐蚀环境耦合作用下的力学响应，有效捕捉了岩石在长期载荷作用下的变形与破坏特征。

4)通过对颗粒间摩擦力和粘结强度对岩石隧道重力方向坍塌应变分析，得出应对长期隧道纵向荷载下，提高围岩中岩石微观缝隙的摩擦力大小和粘结强度能有效遏制蠕变从第二阶段发展至第二阶段。

5) 理解隧道开挖初期应力的集中与重新分配过程，对于评估围岩稳定性、预测变形趋势以及合理设计支护方案具有重要意义。通过监测应力变化，可以及时发现潜在风险，采取针对性措施，确保隧道施工的安全与稳定。