深埋构造混杂岩隧道围岩热-水-力耦合响应

《TUNNELLING AND UNDERGROUND SPACE TECHNOLOGY》：Thermo-hydro-mechanical coupling response of surrounding rock in deep-buried tectonic mélange tunnel

　　

**摘要**

构造混杂体中深埋岩隧道的热-水-力（THM）行为受到强非均质性、渗流和地温的控制，但其耦合影响机制仍认识不足。本研究基于中国西南某深埋岩隧道的地质剖面、岩性分布及现场条件，建立了工程尺度二维非均质连续体模型，用于完全耦合的THM过程分析。通过系统的数值研究，探讨了岩体非均质性、渗流、地温、渗流驱动力及原始岩温对应力重分布、变形模式和渗流特征的影响。结果表明，非均质性主导力学和水力响应，通过引发岩性界面附近的显著应力重分布和渗流局部集中。渗流从根本上重塑了竖向位移分布，增加了拱顶、拱肩和侧墙的沉降，同时抑制了拱底和起拱点的隆起。较高的地温通过热致刚度退化、岩石基质热膨胀及地下水黏度降低，显著放大了变形和达西速度。此外，达西速度与渗流驱动力和原始岩温呈准线性增长关系，其中原始岩温的影响更为显著。这些发现深化了对构造混杂体隧道中工程尺度THM响应模式的认识，并为识别高渗流风险区域提供了实用指导。

**1. 引言**

位于极深位置的深埋隧道不断受到高地温、高渗流压力和高地应力引起的热-水-力（THM）耦合效应影响，导致极其复杂的力学行为。构造混杂体是由长期板块相互作用形成的复杂变形地质体，以岩性多样、混合强烈和高地应力为特征（Festa等，2010；Ogata等，2021）。其断面围岩分布和力学行为表现出极端复杂性，具有显著的空间和特性非均质性（Hsü，1974；Malistani和Aydan，2024；Peng等，2025）。因此，深埋构造混杂体隧道的围岩极易受到开挖扰动的影响，可能引发大变形、局部地下水涌出和高温度危害，最终导致重大财产损失和人员伤亡（Dhang，2016；Fan等，2019；Kucuk，2011）。遗憾的是，目前对深埋构造混杂体隧道灾害形成机制的研究仍严重不足，其THM耦合响应机制仍不清楚（Fan等，2019；Kahraman等，2008；Peng等，2025）。亟需采用先进方法表征此类条件下围岩的力学行为。

在深岩工程研究中，现场监测（Chen等，2024；Gui等，2022；Shen等，2008；Wagner，2019）和实验室试验（Liu等，2025；Liu等，2023b；Wang等，2024a）是表征和分析深埋隧道围岩力学行为的常用方法。然而，现场监测由于人力/材料成本高、数据获取周期长和环境干扰显著，在实现深岩体力学行为的动态分析方面面临重大挑战。实验室试验则受限于试件制备复杂、测试时间长和设备依赖度高，严重影响了深岩体力学行为的表征效率。近年来，随着计算机软硬件技术的快速发展，数值模拟方法，尤其是有限元方法，已被广泛应用于深岩工程领域（Malatip等，2012；Peng等，2023；Wang等，2021；Ying等，2025）。该方法将连续体离散化为单元和节点的模型，通过迭代数值求解来近似解析求解不可行且力学复杂的问题（Becker等，2023；Fries和Belytschko，2010；Zhang等，2023）。与现场监测和实验室试验相比，它具有更大的适应性、灵活性和效率，特别是对于涉及复杂几何形状的多场耦合分析（Beuchler和R?sch，2023；Chen等，2015；García等，2011）。因此，它非常适合研究深埋构造混杂体隧道的THM耦合响应，有助于揭示潜在机制和行为模式，并为安全施工提供理论和技术支持。

近年来，地下工程建设的快速发展促进了深埋隧道围岩力学行为的广泛研究（Duan等，2024；Guo等，2023；Xiaoxu等，2024）。主要研究重点包括损伤演化机制（Wang等，2024b；Xie等，2016；Yang等，2014；Yang等，2017）、能量转化原理（Li等，2018；Li和Weng，2016；Liu等，2024）、结构面效应（Feng等，2019；Hu等，2025；Liu等，2026；Manouchehrian和Cai，2018）、特殊地质影响（Fan等，2024；Lisjak等，2015）、创新本构模型（Ren等，2025；Zhang等，2017）和先进数值方法（Tian等，2024；Tian等，2025；Wu等，2026）。此外，在隧道围岩多物理场耦合分析方面也取得了显著进展。Liu等（2018）利用COMSOL Multiphysics，基于变形方程和孔隙介质流动方程，对中国江水源隧道的超大断面渗流条件进行了水-力耦合模拟。他们分析了围岩中的变形分布、应力场特征和达西速度模式，为隧道施工提供了具体技术指导。Zhang等（2020a）和Zhang等（2020b）纳入Mohr-Coulomb关联流动法则和Hoek-Brown准则，并显式考虑围岩软化特性，对圆形隧道进行了水-力耦合模拟。他们的分析侧重于压缩和剪切应变作用下隧道围岩的塑性力学行为和软化行为，最终揭示了孔隙水压力对墙位移、残余半径和塑性区半径的影响机制。Gang-jian等（2022）通过COMSOL Multiphysics开发了专门针对岩溶隧道围岩的应力-破坏-渗流耦合模型。他们进行了数值模拟，研究了岩体非均质性对应力分布、应变能密度和破坏模式的影响，为岩溶隧道灾害预测提供了理论依据。Jia等（2023）开发了专门针对非Darcy渗流条件下隧道的水-力耦合模型。该模型利用中国二滩水电站引水隧洞进行了验证，重点分析了高水压引水隧洞中衬砌-围岩系统和注浆加固区的力学特性。Huang等（2025）建立了软岩隧道的湿度迁移-裂隙耦合模型，模拟分析了湿度扩散对隧道围岩变形和破坏机制的影响。该研究调查了湿度扩散范围、砂砾层厚度和倾角对围岩湿度时空分布、变形特性和裂隙扩展的影响。

现有研究表明，目前关于隧道围岩的多场耦合研究主要集中于基于断面层理和现场地质条件的水-力耦合模拟。相比之下，同时考虑地层特征和实际工况的THM耦合分析仍然极为有限。此外，构造混杂体作为古海洋演化过程的关键档案，其形成机制是地质分析中的基本研究重点（Du等，2019；Leah等，2022；Li等，2015）。目前研究通常采用现场调查和岩石学分析，对全球典型构造混杂体区域的形成机制进行地质考察（Du等，2019；Zhou等，2024）。迄今为止，关于构造混杂体力学性质的研究仍极为有限，作为隧道围岩的深埋构造混杂体力学特征仍不清楚。这些发现强调了进行深埋构造混杂体隧道THM耦合数值模拟的关键必要性，尤其需要重点研究地质剖面、岩性分布和现场工程条件下温度、渗流和变形的相互作用行为。此类研究将建立坚实的理论框架，并提供实用的技术方案，以确保此类复杂地质条件下的施工安全。

本研究针对中国西南某深埋构造混杂体隧道，开发了考虑该隧道地质剖面、岩性分布和现场工程条件的非均质稳态THM耦合数值模型。模型代表支护安装前的开挖后状态，用于研究围岩的应力分布、变形模式和渗流行为，以及岩体非均质性、渗流、地温、渗流驱动力（SDF）和原始岩温（VRT）对围岩响应的影响。本文其余部分组织如下。第2节介绍THM建模框架，包括控制方程、地质模型和材料参数。第3节描述数值建模过程及模型的半定量验证。第4节展示模拟结果。第5节讨论关键因素对围岩力学和水力响应的影响。第6节总结本研究的主要结论。附录A展示了比例系数的敏感性分析。附录B展示了岩壁温度的敏感性分析。

**2. 材料与方法**

本节介绍用于构建构造混杂体隧道THM数值模型的方法学框架。介绍了在COMSOL Multiphysics中实施的控制方程和多物理场耦合策略，包括固体力学、达西定律和孔隙介质传热之间的相互作用。随后，基于现场隧道断面建立层状结构地质模型，并构建多岩性数值模型以表征非均质岩体。此外，定义了主要岩性的力学、水力和热学参数，并引入附加比例系数以考虑工程尺度下岩体的等效力学和水力性质。还纳入了温度相关的弹性模量函数以表征围岩的热退化行为。

**2.1. THM耦合控制方程**

围岩的THM响应在连续体框架内建模，该框架结合了孔隙固体基质的力学平衡、达西流动与质量守恒以及孔隙介质传热。THM耦合在COMSOL Multiphysics中通过孔隙弹性多物理场耦合（固体力学和达西定律物理接口）和热膨胀多物理场耦合（固体力学和孔隙介质传热物理接口）实现。水力场和热力场还通过定义达西定律物理接口中的温度相关流体属性以及包含由达西速度驱动的对流热传输在孔隙介质传热物理接口中进行耦合。

孔隙介质的准静态动量平衡以散度形式表示，如公式（1）所示（Qarinur等，2020）。

(1) ?·?? + ??V = ??

其中?·为散度算子，??为总应力张量，??V为体力向量（包括重力）。运动学遵循线性和小应变假设，如公式（2）所示（Liu等，2018）。

(2) ?? = 1/2 [??? + (???)?]

其中?为梯度算子，??和??分别为总应变张量和位移向量。在本研究中，计算得到的位移（20–80 mm，即0.020–0.080 m，如第4.1节图7(c)和(d)所示）相较于模型域的特征长度和隧道跨度（分别为200 m和12.82 m，均在第3.1节给出）较小，分析聚焦于工程尺度下应力、变形和渗流的重分布而非大变形运动学。因此，小应变公式被认为适用于本研究。本构行为以热-孔隙弹性形式表示，如公式（3）和（4）所示（Qarinur等，2020）。

(3) ?? = ??:(?? ? ??th) ? ??B·??·??

(4) ??th = ??T·(?? ? ??ref)·??

其中:为双点积算子，??th为热应变张量，??B为Biot系数，??为达西定律物理接口中计算得到的孔隙压力，??T为热膨胀系数，??ref为参考温度（取第3.2节中的TR），??为单位张量。??为弹性刚度张量，其模量包含温度相关弹性模量??(??)（定义于第2.3节公式（12））和弹性模量折减系数kE（定义于第2.3节）。??为孔隙介质传热物理接口中计算得到的温度。

围岩的塑性屈服由 Mohr–Coulomb 准则控制，如公式 (5) 所示（Shen 等，2011）。

(5) $f(\sigma) = \tau - \sigma_n \tan\varphi - c$

其中，$\tau$ 为剪应力，$\sigma_n$ 为法向应力，$\varphi$ 为内摩擦角，$c$ 为粘聚力。

此外，渗流采用包含重力的达西定律进行建模，如公式 (6) 所示（Jiang 等，2022）。

(6) $\mathbf{q} = -\frac{k}{\mu_f(T)} \left[ \nabla p - \rho_f(T) \mathbf{g} \right]$

其中，$\mathbf{q}$ 为达西速度（也称为达西通量或比流量），定义为单位总截面积的体积流量（m/s）。$k$ 为渗透系数（乘以在第 2.3 节中定义的渗透系数放大系数 $k_k$）。$\rho_f(T)$ 和 $\mu_f(T)$ 分别为依赖于温度的流体密度和粘度（详见第 2.3 节）。$T$ 为通过孔隙介质传热物理接口计算的温度。$\mathbf{g}$ 为重力加速度（m/s2），方向沿负 y 轴。重力的水力驱动效应通过 $\rho_f(T) \mathbf{g}$ 项引入。

流体质量守恒方程如公式 (7) 所示（Tao 等，2023）。

(7) $\frac{\partial (\phi \rho_f)}{\partial t} + \nabla \cdot [\rho_f(T) \mathbf{q}] = Q_m$

其中，$\phi$ 为孔隙率，$Q_m$ 为源/汇项。

与体积变形相关的流固耦合通过孔隙弹性多物理场耦合进行处理。对于稳态渗流分析，忽略瞬态储水项 $\frac{\partial (\phi \rho_f)}{\partial t}$，公式 (7) 简化为公式 (8)。

(8) $\nabla \cdot [\rho_f(T) \mathbf{q}] = Q_m$

此外，孔隙介质中的热量传递由公式 (9) 控制（Nasir 等，2014）。

(9) $(\rho C_p)_{\text{eff}} \frac{\partial T}{\partial t} + \rho_f C_{p,f} \mathbf{q} \cdot \nabla T = \nabla \cdot (k_{\text{eff}} \nabla T) + Q$

其中，$Q$ 为热源项，$C_{p,f}$ 为流体的比热容。对流热传输由渗流驱动，此处采用达西速度 $\mathbf{q}$ 表示，与孔隙介质传热物理接口中所采用的一致。$(\rho C_p)_{\text{eff}}$ 和 $k_{\text{eff}}$ 分别为有效体积热容和有效热导率，均采用体积分数（混合规则）近似计算，分别由流体相和固体相的贡献组成，如公式 (10) 和 (11) 所示（Sun 等，2017）。

(10) $(\rho C_p)_{\text{eff}} = \phi \rho_f C_{p,f} + (1 - \phi) \rho_s C_{p,s}$

(11) $k_{\text{eff}} = \phi k_f + (1 - \phi) k_s$

其中，$C_{p,s}$ 为固体比热容，$k_s$ 和 $k_f$ 分别为固体和流体的热导率函数。

该模型通过三种主要机制显式考虑了热-水-力耦合。流固耦合通过孔隙弹性来表示，其中孔隙压力通过有效应力项（公式 (3)）贡献于总应力。该耦合通过连接固体力学和达西定律物理接口的孔隙弹性多物理场耦合实现。此外，热力耦合通过热膨胀引入，其中热应变（公式 (4)）影响围岩的变形和应力状态。该机制通过连接固体力学和孔隙介质传热物理接口的热膨胀多物理场耦合实现。另外，水热耦合通过允许渗流场中存在依赖于温度的流体性质来考虑，其中流体密度和粘度随温度变化（公式 (6)），在达西定律物理接口中实现。此外，传热方程中包含由达西速度驱动的对流热传输（公式 (9)），在孔隙介质传热物理接口中实现。

基于上述控制方程，热-水-力相互作用通过孔隙弹性耦合、热膨胀和依赖于温度的流体性质来实现。本公式假设渗透系数恒定且为达西流，温度场在对应于稳定通风的稳态条件下表示。此外，本研究未考虑显式失效准则。

**2.2 横截面地质模型**

层理结构的数字化模型基于中国西南部某深埋构造混杂岩隧道的实际横截面建立。隧道几何形状和岩性分布来源于工程图纸和现场开挖记录，据此识别了围岩中层理和多种岩性的原位分布。在 AutoCAD 中通过映射岩性界面和层理边界来重建横截面层理结构。生成的 CAD 几何随后导入 COMSOL Multiphysics 以建立数值模型。

在几何重建过程中，为降低数值复杂度和计算成本同时保留主要结构特征，对局部不规则边界进行了平滑处理。考虑到隧道轴线几何形状和加载条件沿纵向变化缓慢，采用二维平面应变假设以捕捉围岩的主要横截面响应。横截面层理模型如图 1 所示。选取了一个具有代表性地质横截面，其中构造混杂岩的层理结构和岩性界面清晰可辨，用于几何重建和岩性建模。

在本研究中，采用具有代表性的地质横截面构建被调查的深埋构造混杂岩隧道的工程尺度数值模型。所选截面保留了主要层理结构和岩性界面，被认为适用于捕捉岩体的主要横截面非均匀性。在连续介质有限元框架内，岩性非均匀性通过空间变化的材料属性来表示，而岩性界面不作为显式力学不连续面建模。因此，无法显式再现界面滑动、张开、块体运动和裂隙控制流等不连续行为，局部响应应被解释为非均匀连续介质中的集中效应，而非真实的界面机制。

**2.3 材料参数与缩放因子**

在等温条件下，矿化云母石英岩 X、矿化云母石英岩 Y 和石榴石-矽线石-黑云母片岩的材料参数详列于表 1。矿化云母石英岩表现为各向异性参数 X 和 Y，代表同一岩层的方向变异性。均质岩参数表示这三种岩性的算术平均值，仅用于数值基准测试。密度、弹性模量、泊松比、粘聚力、内摩擦角、孔隙率和渗透系数来源于实验室试验结果和现有工程数据。矿化云母石英岩（X 和 Y）的热导率、比热容和热膨胀系数选取石英岩报告范围内的平均值（Eddemani 等，2024；Liu 等，2023a；Marciniszyn 和 Sieradzki，2013；Sedara 等，2022）。对于石榴石-矽线石-黑云母片岩，热导率、比热容和热膨胀系数作为其三种组成矿物报告参数范围的等比例平均值计算得出（Chen 等，2022；Ridley 等，2020；Xiang 等，2022）。

表 1 构造混杂岩材料参数。

| 岩石 | 密度 (g/cm3) | 弹性模量 (GPa) | 泊松比 | 粘聚力 (MPa) | 内摩擦角 (°) | 孔隙率 | 渗透系数 (10?2? m2) | 热导率 (W/(m·K)) | 比热容 (J/(kg·K)) | 热膨胀系数 (10??/K) |
|------|-------------|--------------|--------|------------|------------|--------|-------------------|----------------|-----------------|-------------------|
| 矿化云母石英岩 X | 2.69 | 24.60 | 0.25 | 40.15 | 36.00 | 2.8% | 1.00 | 3.50 | 900 | 7.00 |
| 矿化云母石英岩 Y | 2.69 | 41.30 | 0.10 | 30.00 | 33.50 | 3.5% | 3.00 | 5.50 | 1000 | 9.00 |
| 石榴石-矽线石-黑云母片岩 | 2.83 | 32.60 | 0.20 | 21.40 | 31.50 | 4.5% | 10.00 | 3.17 | 867 | 13.83 |
| 均质岩 | 2.74 | 32.83 | 0.18 | 30.52 | 33.67 | 3.6% | 4.67 | 4.06 | 922 | 9.94 |

为考虑工程尺度上破裂和断裂岩体刚度降低（相对于完整的实验室试样），引入了弹性模量折减系数（$k_E = 0.1$），将完整岩石模量转换为等效岩体模量。该值的选取主要基于现场工程评估和隧道现场的变形行为观测。先前数值研究中广泛采用了类似的破裂或扰动岩体模量折减处理方法（Zhao 等，2022；Zhao 等，2021），支持了本研究采用的方法论框架。类似地，为表示工程尺度上破裂岩体的等效渗透系数，引入了渗透系数放大系数（$k_k = 1 \times 10^8$），将基质渗透系数放大以表示工程尺度上破裂岩体的等效水力行为。该处理与先前破裂岩体数值研究中常用方法一致（Liu 等，2022；Qian 和 Xu，2017；Rueda Cordero 等，2019）。所选量级与隧道现场观测到的地下水涌入条件和场地尺度水力响应一致。为评估弹性模量折减系数和渗透系数放大系数对数值结果的影响，在附录 A 中进行了敏感性分析。

此外，围岩的弹性模量在研究温度范围内随温度升高呈现近似线性降低。为定量表征这种热弱化效应，建立了弹性模量的热退化函数，如公式 (12) 所示。

(12) $E(T) = E_1 \left[ 1 - \varphi_E \frac{T_2 - T_1}{(T - T_1)} \right]$

其中，$E(T)$ 为弹性模量的热退化函数，$T_1$ 和 $T_2$ 为实验室条件下的参考温度，$\varphi_E$ 为弹性模量折减率，$E_1$ 为温度 $T_1$ 下的弹性模量。上述参数可通过标准化实验室测试程序确定。

在渗流分析中，使用 COMSOL Multiphysics 材料库提供的依赖于温度的材料属性将水指定为流体。密度和动态粘度均随温度变化。动态粘度随温度升高显著降低。例如，从约 1.00 × 10?3 Pa·s（20°C）降至 5.9 × 10?? Pa·s（45°C）和 4.7 × 10?? Pa·s（60°C）。相比之下，密度在所考察温度范围内的变化相对较小，从约 998 kg/m3（20°C）降至 992 kg/m3（45°C）和 983 kg/m3（60°C）。基于这些依赖于温度的流体属性，通过流体密度和动态粘度的变化在本模拟中将温度对渗流的影响纳入考虑。

**3. 数值建模**

本节介绍深埋构造混杂岩隧道数值模型的建立过程。基于隧道现场的横截面层理结构和原位地质条件，构建了几何模型，随后对各区进行了网格划分和生成。基于现场条件，详细规定了耦合热-水-力分析所需的初始条件和边界条件。此外，设计了多个对比数值模型以研究岩体非均匀性、渗流、地温、渗流驱动力（SDF）和原始岩体温度（VRT）对围岩响应行为的影响。此外，基于现有现场渗流观测和报告涌入数据采用了半定量验证方法。将模拟渗流行为与现场观测在空间定位和涌入量级方面进行了比较，为数值模型的可靠性提供了支持证据。

**3.1 几何建模与网格划分**

对于深埋构造混杂岩隧道的几何模型，三维问题可简化为符合平面应变准则的二维公式。因此，使用 COMSOL Multiphysics 建立了深埋构造混杂岩隧道的二维有限元模型。隧道横截面呈马蹄形，尺寸为 12.82 m × 10.87 m，埋深 $H_0 = 500$ m。基于圣维南原理，计算域定义为 200 m × 200 m（$L_0 = 200$ m），以减少边界效应并容纳构造混杂岩的非均匀性。数值模型中心取为坐标原点，与隧道中心重合。

在数值模拟中，采用多级网格离散化策略，在计算效率与解的精度之间取得适当平衡。在隧道几何体附近的近场区域应用局部细化网格，而在远场区域使用逐渐粗化的网格。这种网格布置提供了近场区域（尤其是隧道边界和岩性界面处）足够的数值分辨率，同时在远场域保持计算效率。深埋构造混杂岩隧道的数值模型和网格配置，包括全域模型、近场区域和隧道几何体，如图 2 所示。

值得注意的是，数值模型中的细化网格区域对应于包含横截面层理结构的近场区域，仅由矿化云母石英岩 X、矿化云母石英岩 Y 和石榴石-矽线石-黑云母片岩组成（定义于第 2.3 节）。材料属性按照图 1 所示的地质建模程序和第 2.2 节的描述进行分配。远场区域使用粗网格离散化，建模为均质岩。数值模型在二维平面应变假设下建立，以表示隧道系统的横截面热-水-力响应。

**3.2 初始条件与边界条件**

在大尺度数值变形分析建模中，初始条件中必须考虑重力对应力分布的影响。此外，模型自重产生不可忽略的初始地应力，直接影响隧道开挖结果，需要进行应力平衡分析。因此，数值模型包含两个分析步骤。在开挖前分析步骤中，数值建立模型的初始应力场，随后将其作为开挖后分析步骤的预应力条件施加。数值模型的初始和边界条件如图 3 所示。

数值模型的初始条件与边界条件：(a) 力学条件；(b) 水力条件；(c) 热力条件。如图3(a)所示，下边界设为固定约束，隧道边界设为无牵引力边界。此外，模型的上、左、右边界施加了初始地应力条件。左右边界应力需考虑深度依赖效应，因此定义为竖直位置的线性函数。应力由隧道埋深和模型尺寸决定，公式如式(13)、(14)所示。(13) ??upper_y = ???g(H? ? L?/2) (14) ??left,right_x(y) = k????g(H? ? y) 其中，??upper_y、??left_x和??right_x分别为上、左、右边界应力，???为均质岩石密度（见第2.3节表1），g为重力加速度，H?和L?分别为隧道埋深和模型高度（均见第3.1节）。k?为水平侧压力系数（取值为1.5），代表中等构造应力环境。y为模型竖直坐标，数值模型的几何原点（y = 0 m）位于隧道几何中心。符号约定遵循标准连续介质力学记法，压应力为正。如图3(b)所示，Darcy定律控制的流体流动由重力决定的压力水头驱动。左右边界设为无流量边界。隧道边界施加零水头以表示具有自由排水条件的无衬砌隧道。此外，压力水头分布需考虑深度依赖静水压力，定义为深度的线性函数，公式如式(15)所示。(15) H_p(y) = k_SDF(H? ? y) 其中，H_p(y)为数值模型中的压力水头。k_SDF为SDF系数，作为Darcy公式中重力驱动压力水头项的无量纲缩放系数。k_SDF = 1.0代表参考水力驱动条件，k_SDF < 1.0和k_SDF > 1.0分别代表相对较弱和较强的水力驱动条件。不同模型中使用的数值见表2、表3（见第3.3节）。此外，上下边界的压力水头由其深度给出，公式如式(16)、(17)所示。(16) H_upper_p = k_SDF(H? ? L?/2) (17) H_lower_p = k_SDF(H? + L?/2) 其中，H_upper_p和H_lower_p分别为上下边界的压力水头。表2 模型的多物理场耦合条件。模型 | 围岩 | 渗流驱动力 | 围岩初始地温 | 岩壁温度 THM-TM模型 | 构造混合岩 | 1.0 | 60°C | 45°C THM-HR模型 | 均质岩石 | 1.0 | 60°C | 45°C TM-TM模型 | 构造混合岩 | 0 | 60°C | 45°C HM-TM模型 | 构造混合岩 | 1.0 | 25°C | 25°C 表3 SDF模型和VRT模型的多物理场耦合条件。模型 | 渗流驱动力 | 围岩初始地温 | 岩壁温度 SDF 0.8模型 | 0.8 | 60°C | 45°C SDF 1.0模型 | 1.0 | 60°C | 45°C SDF 1.2模型 | 1.2 | 60°C | 45°C SDF 1.4模型 | 1.4 | 60°C | 45°C VRT 50模型 | 1.0 | 50°C | 35°C VRT 60模型 | 1.0 | 60°C | 45°C VRT 70模型 | 1.0 | 70°C | 55°C VRT 80模型 | 1.0 | 80°C | 65°C 如图3(c)所示，根据隧道现场监测数据，开挖后围岩温度约为60°C。经约4小时机械通风后，监测到的岩壁温度下降并趋于稳定，约为45°C。基于此监测到的冷却过程，本研究采用的热力条件被视为准稳定运行状态，而非开挖后立即出现的瞬态热演化。因此，稳态热力条件适用于评估该准稳定条件下的后续THM响应。因此，初始岩体温度和外边界温度设为T_R = 60°C，隧道边界温度定义为T_C = 45°C。为评估T_R与T_C之间温差对数值结果的影响，在附录B中进行了敏感性分析。3.3 模型设计与监测点选取 在深埋构造混合岩隧道中，岩体非均质性、渗流和地温直接影响围岩的力学行为。为研究这些影响因素，建立了多个数值模型，包括构造混合岩的热-水-力耦合模型（THM-TM模型）、均质岩石的热-水-力耦合模型（THM-HR模型）、构造混合岩的热-力耦合模型（TM-TM模型）以及构造混合岩的水-力耦合模型（HM-TM模型），多物理场耦合条件详见表2。SDF（k_SDF）值1.0对应第3.2节定义的参考水力驱动条件。选取的0.8、1.0、1.2和1.4值分别代表相对较弱、参考、中等较强和较强的重力驱动压力水头条件。这些值用于开展水力驱动强度影响的参数分析，而非直接表示特定现场水头的测量值。各模型对比分析了岩体非均质性、渗流和地温对深埋构造混合岩隧道中应力分布、变形模式和渗流行为的影响，对比关系如图4所示。图4. 模型间对比关系（箭头标注因素）。为研究SDF和VRT对围岩响应行为的影响，建立了两组对比模型。针对不同SDF条件（系数定义见第3.2节）和不同VRT条件建立了构造混合岩的THM耦合模型。在VRT模型中，岩壁温度随围岩初始地温同步调整，使其间温差保持不变。这两组分别称为SDF模型和VRT模型。SDF模型和VRT模型的多物理场耦合条件汇总于表3。这些模型的围岩均为构造混合岩。在深埋构造混合岩隧道的数值分析中，横截面划分为拱顶、拱肩、边墙、起拱线和仰拱等区域。为分析THM耦合效应下不同区域围岩的复杂响应行为，在隧道边界选取了代表性监测点（1~8号点）进行对比分析，各点位置及对应区域如图5所示。具体而言，1号和8号点分别代表拱顶和仰拱，位于隧道中心线上。2号和3号点（左拱肩和右拱肩）、4号和5号点（左边墙和右边墙）、6号和7号点（左起拱线和右起拱线）对称分布于中心线两侧。图5. 监测点位置及对应区域。3.4 模型验证 为半定量评估数值结果与现有现场渗流观测之间的一致性，将模拟结果与现场渗流观测及报告入流量数据进行对比。现场渗流观测和Darcy速度数值结果如图6所示。图6. 现场渗流观测与模拟Darcy速度分布对比。红色虚线表示现场观测到的渗流位置。黑色虚线划定高渗流潜力区域。红色和蓝色分别代表正向（向上/向右）和负向（向下/向左）方向：(a) 现场渗流观测；(b) 竖直Darcy速度数值结果；(c) 水平Darcy速度数值结果。如图6(a)所示，附近隧道断面的现场观测表明，可见漏水主要发生在拱肩处。在图6(b)和6(c)中，数值结果表明地下水的入流具有局部化特征。高Darcy速度区域主要预测在拱肩和起拱线处。与现场观测相比，数值结果显示出部分一致的空间分布模式。现场渗流主要在隧道拱肩处观测到，与预测的拱肩处优先渗流局部化特征相符。数值模型还识别出起拱线附近的额外潜在渗流区域，这可能代表非均质渗透率条件下的保守预测。需要注意的是，数值模型基于所研究隧道的代表性地质断面建立（如图1所示），该断面的构造混合岩层理结构清晰可辨，适合进行几何和岩性建模。相比之下，用于验证的现场照片取自附近开挖面，该处局部渗流表现更为清晰暴露，因此更适合识别优先渗流位置。由于现有工程记录的局限性，模型断面与拍摄断面之间无法建立严格的逐一对应关系。因此，本验证旨在评估优先渗流局部化和入流量级顺序的一致性，而非对同一开挖面上渗流分布的严格逐点复现。渗流强度通过基于模拟Darcy速度的半定量估算进一步评估。数值结果表明，高渗流潜力区域（A~E区域）的代表性Darcy速度量级v_Darcy约为1~10 mm/s（0.001~0.010 m/s）。考虑到现场渗流以局部漏水而非均匀入流的形式出现，基于观测到的可见湿痕宽度保守地采用有效渗流宽度L_eff = 0.05 m（可见湿痕宽度的下限）。对于入流换算，有效渗流面积定义为L_eff × 10 m，其中10 m为现场报告入流量数据中采用的隧道长度。每10 m隧道长度的等效地下水入流量估算为Q_10m = v_Darcy·L_eff·10 m，预测入流量范围约为0.0005~0.0050 m3/s·10 m（30~300 L/min·10 m）。预测入流量范围高于现场报告的入流强度25~125 L/min·10 m，但量级一致，表明该模型提供了对地下水入流强度的保守半定量估算。数值模型在优先渗流位置和入流量级方面与现场观测具有半定量一致性。渗流局部化的空间一致性和入流强度的半定量一致为非均质围岩中模拟的优先渗流局部化和等效入流量级提供了支持。由于现有现场记录中断面对应的局限性，本验证旨在评估工程尺度上的一致性，而非精确的逐面对复现。因此，力学和热力响应的解释在采用的THM建模框架和假设内进行，而非直接通过现场测量进行验证。4. 结果 本节呈现THM-TM模型的模拟结果，包括von Mises应力、位移（竖直和水平）以及Darcy速度（竖直和水平）的等值线图，以及监测点（位置见图5）处的计算结果。此外，还提供了THM-HR、THM-TM、TM-TM和HM-TM模型（定义见第3.3节）的等值线图对比和监测点结果，以及SDF和VRT模型（定义见第3.3节）的监测点结果。在应力等值线图中，von Mises等效应力作为标量指标用于可视化应力集中，而塑性屈服则根据第2.1节所述的Mohr-Coulomb准则进行评估。4.1 THM-TM模型结果 THM-TM模型结果如图7所示。数据符号表示方向性，绝对值表示量级。隧道开挖通过首先建立初始地应力平衡状态，然后移除隧道区域来模拟。如图7(a)所示，开挖前初始von Mises应力分布相对均匀，范围为10~20 MPa，最大应力19.02 MPa位于不同岩体之间的岩性界面处。如图7(b)所示，开挖后拱顶（37.07 MPa）和两侧起拱线（35.60 MPa和36.98 MPa）的von Mises应力显著高于其他区域。中心线两侧的应力呈不对称分布，最大应力集中在拱顶附近的岩性界面处（40.92 MPa）。结果表明开挖在围岩中引起了显著的应力重分布。最大应力从开挖前的19.02 MPa增加到开挖后的40.92 MPa，增幅为115.14%，反映了开挖引起的隧道边界周围应力集中效应。如图7(c)所示，拱顶（?56.77 mm）、两侧拱肩（?44.02 mm和?44.22 mm）以及两侧边墙（?22.20 mm和?29.51 mm）均出现显著的竖直向下位移。两侧起拱线的竖直向下位移（?1.14 mm和?3.79 mm）相对较小，而仰拱（18.24 mm）表现出显著的竖直向上位移。最大向下位移出现在拱顶（?56.99 mm），最大向上位移出现在仰拱（18.70 mm）。如图7(d)所示，中心线左侧和右侧区域分别表现出向右和向左的水平位移，而拱顶（?3.37 mm）和仰拱（2.91 mm）处的水平位移相对较小。最大向右位移出现在左边墙（71.19 mm），最大向左位移出现在右边墙（?70.72 mm）。

在图7(e)和图7(f)中，高渗潜势区局部分布在拱顶（A区和E区）、右肩拱（B区）及两侧腰线（C区和D区）的特定岩层单元内，表现出显著的水流集中现象。A区（?6.57 mm/s和?2.83 mm/s）、B区（?4.84 mm/s和?9.48 mm/s）、C区（7.28 mm/s和?3.04 mm/s）以及D区（6.19 mm/s和5.38 mm/s）的达西速度远高于E区（?2.61 mm/s和0.90 mm/s）。下载：下载高分辨率图片（558KB）下载：下载全尺寸图片

图7. THM-TM模型的结果。黑点表示监测点；虚线勾勒高渗潜势区边界，标注数值表示各区内的最大达西速度；紫色和白色星号表示最大值。红色和蓝色分别表示正方向（向上/向右）和负方向（向下/向左）：(a) 初始地应力冯·米塞斯应力；(b) 开挖后冯·米塞斯应力；(c) 竖向位移；(d) 水平位移；(e) 竖向达西速度；(f) 水平达西速度。

**4.2 THM-HR与THM-TM模型的对比**

THM-HR和THM-TM模型的冯·米塞斯应力如图8所示。在图8(a)和(b)中，THM-HR模型在两侧腰线处均呈现最大应力（均为40.28 MPa），而THM-TM模型的最大应力出现在拱顶附近的岩性界面处（40.92 MPa）。与THM-HR模型相比，THM-TM模型的最大应力增加了1.59%（从40.28 MPa增至40.92 MPa）。在图8(c)中，拱底处应力（从28.59 MPa降至23.52 MPa）受岩性非均质性的影响更为显著（降幅达17.71%），该降幅明显大于其他区域（降幅为5.19%至8.72%）。下载：下载高分辨率图片（283KB）下载：下载全尺寸图片

图8. THM-HR和THM-TM模型的冯·米塞斯应力。黑点表示监测点，紫色星号表示最大值：(a) THM-HR模型的冯·米塞斯应力；(b) THM-TM模型的冯·米塞斯应力；(c) 监测点处冯·米塞斯应力值及百分比变化。

THM-HR和THM-TM模型的竖向位移如图9所示。数据符号表示方向，绝对值表示大小。在图9(a)和(b)中，THM-HR模型在两侧腰线处呈现中等程度的竖向向下位移（均为?8.56 mm），而THM-TM模型在相同区域的位移可忽略不计（?1.14 mm和?3.79 mm）。此外，两个模型的竖向向下最大位移均出现在拱顶（?53.28 mm）或拱顶附近（?56.99 mm），竖向向上最大位移均出现在拱底（12.60 mm）或拱底附近（18.70 mm）。与THM-HR模型相比，THM-TM模型的竖向向上最大位移和竖向向下最大位移分别增加了48.41%（从12.60 mm增至18.70 mm）和6.96%（从?53.28 mm增至?56.99 mm）。在图9(c)中，拱顶（点1）和两侧肩拱（点2和点3）表现出较大的竖向位移，但百分比变化较小（分别为6.56%、3.00%和3.45%）；而两侧腰线（点6和点7）及拱底（点8）的百分比变化则显著较大（分别为?86.70%、?55.74%和44.75%），尽管其位移量较小。下载：下载高分辨率图片（242KB）下载：下载全尺寸图片

图9. THM-HR和THM-TM模型的竖向位移。黑点表示监测点，紫色星号表示最大值。红色和蓝色分别表示正方向（向上）和负方向（向下）：(a) THM-HR模型的竖向位移；(b) THM-TM模型的竖向位移；(c) 监测点处竖向位移值及百分比变化。

THM-HR和THM-TM模型的水平位移如图10所示。数据符号表示方向，绝对值表示大小。在图10(a)和(b)中，两个模型均在左侧壁处呈现最大向右水平位移，在右侧壁处呈现最大向左水平位移。与THM-HR模型相比，THM-TM模型的最大向右水平位移和最大向左水平位移分别增加了6.96%（从66.56 mm增至71.19 mm）和6.25%（从?66.56 mm增至?70.72 mm）。在图10(c)中，两侧壁（点4和点5）的水平位移大于两侧肩拱（点2和点3）和两侧腰线（点6和点7）的水平位移。两侧壁（点4和点5）的水平位移百分比变化（分别为7.18%和6.06%）高于两侧腰线（点6和点7，分别为2.87%和?1.19%）。两侧肩拱（点2和点3）的变化相对较小（分别为?0.58%和0.87%）。对于拱顶（点1）和拱底（点8），THM-HR模型中的水平位移均为零，因此未报告这两个区域的百分比变化。虽然THM-TM模型中对应的水平位移不为零，但其幅值仍然很小（点1处为?3.37 mm，点8处为2.91 mm）。下载：下载高分辨率图片（265KB）下载：下载全尺寸图片

图10. THM-HR和THM-TM模型的水平位移。黑点表示监测点，紫色星号表示最大值。红色和蓝色分别表示正方向（向右）和负方向（向左）：(a) THM-HR模型的水平位移；(b) THM-TM模型的水平位移；(c) 监测点处水平位移值及百分比变化。

THM-HR和THM-TM模型的达西速度如图11所示。数据符号表示方向，绝对值表示大小。THM-HR模型的高渗潜势区集中在拱顶、两侧壁及两侧腰线。THM-TM模型的高渗潜势区分布在拱顶（A区和E区）、右肩拱（B区）及两侧腰线（C区和D区）的特定岩层单元内。与THM-HR模型相比，THM-TM模型在高渗潜势区内的竖向达西速度平均绝对值和水平达西速度平均绝对值分别增加了127.27%（从2.42 mm/s增至5.50 mm/s）和57.45%（从2.75 mm/s增至4.33 mm/s）。下载：下载高分辨率图片（582KB）下载：下载全尺寸图片

图11. THM-HR和THM-TM模型的达西速度。黑色虚线勾勒高渗潜势区边界，标注数值表示各区内的最大达西速度：(a) THM-HR模型的竖向达西速度；(b) THM-TM模型的竖向达西速度；(c) THM-HR模型的水平达西速度；(d) THM-TM模型的水平达西速度。

**4.3 TM-TM与THM-TM模型的对比**

TM-TM和THM-TM模型的冯·米塞斯应力如图12所示。在图12(a)和(b)中，两个模型的最大应力均出现在拱顶附近的岩性界面处。与TM-TM模型相比，THM-TM模型的最大应力降低了1.73%（从41.64 MPa降至40.92 MPa）。在图12(c)中，拱顶（点1）和两侧腰线（点6和点7）的应力大于其他区域，但所有区域的变化均较小（?5.20%至?1.80%）。下载：下载高分辨率图片（299KB）下载：下载全尺寸图片

图12. TM-TM和THM-TM模型的冯·米塞斯应力。黑点表示监测点，紫色星号表示最大值：(a) TM-TM模型的冯·米塞斯应力；(b) THM-TM模型的冯·米塞斯应力；(c) 监测点处冯·米塞斯应力值及百分比变化。

TM-TM和THM-TM模型的竖向位移如图13所示。数据符号表示方向，绝对值表示大小。在图13(a)和(b)中，与TM-TM模型相比，THM-TM模型在拱顶、两侧肩拱及两侧壁处的竖向向下位移幅值均有所增大，而拱底处的竖向向上位移则有所减小。两侧腰线从TM-TM模型中的较大竖向向上位移（17.02 mm和15.61 mm）转变为THM-TM模型中的较小竖向向下位移（?1.14 mm和?3.79 mm）。两个模型的最大竖向向上位移和最大竖向向下位移分别出现在拱底附近和拱顶附近。与TM-TM模型相比，THM-TM模型的最大竖向向上位移降低了50.78%（从37.99 mm降至18.70 mm），最大竖向向下位移增加了56.61%（从?36.39 mm增至?56.99 mm）。在图13(c)中，竖向位移的百分比变化显著（?93.31%至575.86%），尤其是两侧壁（575.86%和212.22%），这是由于TM-TM模型中的基线值较小所致。其他区域的百分比变化则处于较窄的范围内（?93.91%至87.74%）。下载：下载高分辨率图片（260KB）下载：下载全尺寸图片

图13. TM-TM和THM-TM模型的竖向位移。黑点表示监测点，紫色星号表示最大值。红色和蓝色分别表示正方向（向上）和负方向（向下）：(a) TM-TM模型的竖向位移；(b) THM-TM模型的竖向位移；(c) 监测点处竖向位移值及百分比变化。

TM-TM和THM-TM模型的水平位移如图14所示。数据符号表示方向，绝对值表示大小。在图14(a)和(b)中，两个模型的最大向右和向左水平位移分别出现在左侧壁和右侧壁。与TM-TM模型相比，THM-TM模型的最大向右水平位移降低了0.17%（从71.31 mm降至71.19 mm），最大向左水平位移增加了1.27%（从?69.83 mm增至?70.72 mm）。在图14(c)中，拱顶处水平位移的百分比变化（20.60%）大于其他区域（?0.54%至4.86%），这是由于TM-TM模型中的基线值较小所致。下载：下载高分辨率图片（279KB）下载：下载全尺寸图片

图14. TM-TM和THM-TM模型的水平位移。黑点表示监测点，紫色星号表示最大值。红色和蓝色分别表示正方向（向右）和负方向（向左）：(a) TM-TM模型的水平位移；(b) THM-TM模型的水平位移；(c) 监测点处水平位移值及百分比变化。

**4.4 HM-TM与THM-TM模型的对比**

HM-TM和THM-TM模型的冯·米塞斯应力如图15所示。在图15(a)和(b)中，两个模型的最大应力均出现在拱顶附近的岩性界面处。与HM-TM模型相比，THM-TM模型的最大应力降低了0.73%（从41.22 MPa降至40.92 MPa）。在图15(c)中，拱顶（点1）和两侧腰线（点6和点7）的应力大于其他区域，但所有区域的变化均较小（?1.62%至?0.43%）。下载：下载高分辨率图片（283KB）下载：下载全尺寸图片

图15. HM-TM和THM-TM模型的冯·米塞斯应力。黑点表示监测点，紫色星号表示最大值：(a) HM-TM模型的冯·米塞斯应力；(b) THM-TM模型的冯·米塞斯应力；(c) 监测点处冯·米塞斯应力值及百分比变化。

HM-TM和THM-TM模型的竖向位移如图16所示。数据符号表示方向，绝对值表示大小。在图16(a)和(b)中，与HM-TM模型相比，THM-TM模型的竖向位移幅值均有所增大。两个模型的最大竖向向上位移和向下位移分别出现在拱底附近和拱顶附近。与HM-TM模型相比，THM-TM模型的最大竖向向上位移和最大竖向向下位移分别增加了31.50%（从14.22 mm增至18.70 mm）和41.55%（从?40.26 mm增至?56.99 mm）。在图16(c)中，左侧腰线（2078.64%）和右侧腰线（93.52%）的竖向位移百分比变化显著大于其他区域（31.23%至45.81%），这是由于HM-TM模型中的基线值较小所致。下载：下载高分辨率图片（269KB）下载：下载全尺寸图片

图16. HM-TM和THM-TM模型的竖向位移。黑点表示监测点，紫色星号表示最大值。红色和蓝色分别表示正方向（向上）和负方向（向下）：(a) HM-TM模型的竖向位移；(b) THM-TM模型的竖向位移；(c) 监测点处竖向位移值及百分比变化。

HM-TM和THM-TM模型的水平位移如图17所示。数据符号表示方向，绝对值表示大小。在图17(a)和(b)中，与HM-TM模型相比，THM-TM模型的水平位移幅值均有所增大。两个模型的最大向右和向左水平位移分别出现在左侧壁和右侧壁。与HM-TM模型相比，THM-TM模型的最大向右水平位移和最大向左水平位移分别增加了38.85%（从51.27 mm增至71.19 mm）和38.94%（从?50.90 mm增至?70.72 mm）。

如图17(c)所示，拱顶（1号点）和拱底（8号点）的水平位移明显小于其他区域，而所有区域水平位移的变化百分比处于较窄的范围内（38.86%至40.96%）。

图17. HM-TM和THM-TM模型的水平位移。黑点表示监测点，紫色星号表示最大值。红色和蓝色分别代表正向（向右）和负向（向左）方向：(a) HM-TM模型的水平位移；(b) THM-TM模型的水平位移；(c) 监测点处的水平位移数值及变化百分比。

HM-TM和THM-TM模型的达西速度如图18所示。数据的正负号指示方向，绝对值量化大小。在图18(a)、(b)、(d)和(e)中，两个模型在特定岩性单元内的拱顶（A区和E区）、右肩（B区）以及两侧泉线（C区和D区）均表现出相同的高渗流势区域。在图18(c)和(f)中，A区至E区的垂直达西速度和水平达西速度的变化百分比分别为100.90%至113.74%和100.33%至112.01%。

图18. HM-TM和THM-TM模型的达西速度。黑色虚线勾勒高渗流势区域，标注数值代表各区内的最大达西速度：(a) HM-TM模型的垂直达西速度；(b) THM-TM模型的垂直达西速度；(c) 垂直达西速度的数值及变化百分比；(d) HM-TM模型的水平达西速度；(e) THM-TM模型的水平达西速度；(f) 水平达西速度的数值及变化百分比。

4.5 SDF和VRT模型的结果

SDF模型在监测点和高渗流势区域的结果如图19所示。模型中监测点和高渗流势区域的位置与THM-TM模型中的位置相同。图中绘制的位移和达西速度数值为幅值，用于比较SDF对响应强度的影响，而相应方向则根据THM-TM模型的原始带符号结果进行解释。在图19(a)和(c)中，随着SDF的增加，冯·米塞斯应力和水平位移基本保持不变。在图19(b)中，垂直位移对SDF表现出明显的敏感性。拱顶（1号点）、两侧肩（2号和3号点）、两侧边墙（4号和5号点）以及两侧泉线（6号和7号点）的位移幅值通常随SDF的增加呈近似线性增长。对于6号点，在SDF 0.8模型中原始带符号的垂直位移为向上，但在SDF 1.0、SDF 1.2和SDF 1.4模型中变为向下。由于图19(b)旨在比较SDF对变形幅值的影响，垂直位移以幅值形式绘制。然而，由于SDF 0.8模型中的位移方向与SDF 1.0、SDF 1.2和SDF 1.4模型相反，该点在处理中被视为方向反转响应。相比之下，拱底（8号点）显示出相反的趋势，随着SDF增加，隆起幅值从22.42 mm降至9.87 mm。这种不均匀响应表明，垂直位移不仅受SDF幅值控制，还受不同隧道位置周围的局部渗流力方向、孔隙压力重新分布和应力重新分布的控制，详见第5.4节的进一步讨论。在图19(d)和(e)中，垂直和水平达西速度均随SDF的增加呈近似线性增长，且初始达西速度较大的区域显示出更大的绝对增量。

图19. SDF模型在监测点和高渗流势区域的结果。图中绘制的位移和达西速度数值为幅值，其方向根据THM-TM模型的原始带符号结果解释：(a) 冯·米塞斯应力；(b) 垂直位移幅值；(c) 水平位移幅值；(d) 垂直达西速度幅值；(e) 水平达西速度幅值。

VRT模型在监测点和高渗流势区域的结果如图20所示。模型中监测点和高渗流势区域的位置与THM-TM模型中的位置相同。图中绘制的位移和达西速度数值为幅值，用于比较VRT对响应强度的影响，而相应方向则根据THM-TM模型的原始带符号结果进行解释。在图20(a)、(b)和(c)中，随着VRT的增加，冯·米塞斯应力、垂直位移和水平位移基本保持不变。在图20(d)和(e)中，垂直和水平达西速度均随VRT的增加呈近似线性增长，且初始达西速度较大的区域显示出更大的绝对增量。

图20. VRT模型在监测点和高渗流势区域的结果。图中绘制的位移和达西速度数值为幅值，其方向根据THM-TM模型的原始带符号结果解释：(a) 冯·米塞斯应力；(b) 垂直位移幅值；(c) 水平位移幅值；(d) 垂直达西速度幅值；(e) 水平达西速度幅值。

5 讨论

本节讨论第4节中呈现的数值模拟结果，以进一步阐明岩体非均质性、渗流、地温和SDF及VRT对深埋构造混杂岩隧道周围岩体响应行为的影响。特别关注多场耦合条件下控制应力重新分布、变形特征和渗流局部化的机制。通过解读数值结果，本节旨在明确地质非均质性和耦合THM过程如何共同影响隧道变形模式和地下水涌水行为。

5.1 岩体非均质性的影响

图8的结果表明，构造混杂岩隧道中的应力分布受岩体非均质性的影响。非均质性通常降低了应力的大小，尤其在拱底处。其对整体分布模式的影响相对较小。拱顶处的岩石因大埋深引起的较高原位应力而承受显著的垂直压缩，而泉线处的岩石经历集中的高压缩和剪切应力。因此，这些区域的应力明显大于其他区域。此外，岩体非均质性加剧了拱顶处的应力集中，导致最大应力发生在岩性界面处。这一现象源于该界面处两类岩体之间力学性质存在显著差异，导致不协调变形并增加最大应力。当应力在拱顶岩性界面等关键位置高度集中时，其他区域的应力相应降低，表明围岩内部存在应力重新分布。

图9和图10的结果表明，岩体非均质性增大了拱底处的底板隆起变形和边墙处的径向收敛变形，同时减小了泉线处的垂直位移。岩体非均质性对垂直位移的影响大于对水平位移的影响。这些不同岩性对开挖引起的应力变化表现出完全不同甚至相反的力学响应，因为不同刚度的岩石表现出不同的变形响应。这些岩石之间的相互作用和不协调导致了拱底处显著的底板隆起变形和边墙处的径向收敛变形。

图11的结果表明，岩体非均质性显著改变了渗流的分布，并在拱顶、肩部和泉线的特定区域导致明显的渗流集中。这些结果表明，岩体非均质性在控制非均质围岩中的渗流集中方面起着重要作用。

5.2 渗流的影响

图12的结果表明，渗流可轻微降低构造混杂岩隧道中的应力。这一行为可归因于有效应力原理。渗流的存在诱发了孔隙压力分布，降低了有效应力，而岩石应力的整体分布和大小仍主要受开挖引起的应力重新分布控制。

图13和图14的结果表明，渗流显著改变了垂直位移模式，增大了拱顶、肩部和边墙处的沉降，同时减小了拱底和泉线处的底板隆起。相比之下，径向收敛一般仅略有增加。这些行为主要与渗流引起的孔隙压力有关，孔隙压力降低了有效应力，从而放大了拱顶和肩部处的沉降变形。此外，虽然渗流作用下拱顶处水平位移的变化率较大，但变化前后的绝对值仍然可以忽略不计（小于1.00 mm）。其他区域的径向收敛变化均低于5.00%，表明渗流对隧道围岩径向收敛变形的总体影响不显著。

5.3 地温的影响

图15的结果表明，地温对构造混杂岩隧道中的应力仅有有限影响，应力值整体略有降低。最大百分比降幅低于2.00%，表明在所研究范围内热效应对应力幅值的影响相对较弱。这种小幅降低主要与围岩弹性模量随温度升高而降低有关，这略微削弱了岩体刚度。此外，热膨胀可能对隧道周围的局部应力重新分布有所贡献。

图16和图17的结果表明，地温增大了构造混杂岩隧道中的变形，且影响程度相似（垂直方向：31.23%至45.81%，由于泉线处变化异常大故予以排除；水平方向：38.86%至40.96%）。这一行为主要归因于围岩弹性模量随温度升高而降低，这削弱了岩体刚度。此外，围岩的热膨胀改变了隧道周围的应力状态，导致变形增加。弹性模量退化和热膨胀的综合效应导致了径向收敛变形的显著增加。

图18的结果表明，地温并未改变渗流区域的分布，但相对于基准工况，达西速度增加了不止一倍（垂直方向：100.90%至113.74%；水平方向：100.33%至112.01%）。这一行为主要归因于水动力黏度对温度的依赖性。随着温度升高，水动力黏度降低，从而提高了水力传导性，进而根据达西定律增加了达西速度。模型中使用的水动力黏度随温度的变化关系在第2.3节中有描述。此外，域内温度变化相似，导致不同区域的达西速度增加幅度相当。

5.4 SDF和VRT的影响

图19和图20的结果表明深埋构造混杂岩隧道对SDF和VRT的力学响应。SDF与拱顶、肩部、边墙和泉线处的沉降变形呈正相关，而与拱底处的隆起变形呈近似线性反比关系。两个因素均显著增加了达西速度，且呈近似线性趋势。与达西速度较低的区域相比，达西速度较高的区域表现出更大的流体动力响应。

图19(b)的结果表明，不同区域的不同垂直位移趋势可归因于渗流力方向、孔隙压力重新分布和隧道周围局部应力重新分布的综合效应。随着SDF增加，重力驱动的压力梯度增强，导致更强的渗流力效应和更显著的孔隙压力重新分布。在拱顶、肩部、边墙和大部分泉线位置，这种增强的水力作用促进了沉降变形。相比之下，拱底受不同的局部水力-力学平衡支配。拱底的底板隆起主要与开挖引起的应力释放有关，而SDF的增加改变了隧道底部附近的孔隙压力场和局部有效应力状态，从而减小了向上隆起的趋势。这解释了为什么拱底隆起从22.42 mm降至9.87 mm，随着SDF的增加。泉线处原始带符号垂直位移方向的变化进一步表明泉线附近主导变形机制发生了局部转变。在SDF 0.8的相对较低条件下，开挖引起的应力释放在局部占主导，产生向上变形。

随着渗流驱动力系数（SDF）增大至1.0–1.4，增强的渗流力效应和孔隙压力重分布促进了向下沉降，导致拱腰处发生位移方向反转。为了分析SDF对不同区域竖向变形的影响程度，建立了一个名为竖向水-力变形指数（VHMD Index）的指标。同时，还制定了四个互补指标：竖向驱动-达西响应指数（VDDR Index）、水平驱动-达西响应指数（HDDR Index）、竖向热-达西响应指数（VTDR Index）和水平热-达西响应指数（HTDR Index），以评估SDF和VRT对各区域达西速度的影响。这些指数根据SDF和VRT模型的模拟结果（见图19、图20）采用最小二乘法计算，公式如式(18)、(19)、(20)所示。

(18) $I_{VHMD} = \frac{n\sum_{i=1}^{n}(a_{SD_i} d_{V_i}) - \sum_{i=1}^{n} a_{SD_i} \sum_{i=1}^{n} d_{V_i}}{n\sum_{i=1}^{n}(a_{SD_i})^2 - (\sum_{i=1}^{n} a_{SD_i})^2}$

(19) $I_{VDDR,HDDR} = \frac{n\sum_{i=1}^{n}(a_{SD_i} v_{V,H_i}) - \sum_{i=1}^{n} a_{SD_i} \sum_{i=1}^{n} v_{V,H_i}}{n\sum_{i=1}^{n}(a_{SD_i})^2 - (\sum_{i=1}^{n} a_{SD_i})^2}$

(20) $I_{VTDR,HTDR} = \frac{n\sum_{i=1}^{n}(T_{VR_i} v_{V,H_i}) - \sum_{i=1}^{n} T_{VR_i} \sum_{i=1}^{n} v_{V,H_i}}{n\sum_{i=1}^{n}(T_{VR_i})^2 - (\sum_{i=1}^{n} T_{VR_i})^2}$

其中，$I$ 为指数，$a_{SD}$ 为SDF转换的加速度（定义见第3.3节），$T_{VR}$ 为VRT模型的围岩原始温度，$d_V$ 为竖向变形，$v_V$ 和 $v_H$ 分别为竖向和水平达西速度，$n$ 为计算数据集的基数。监测点和渗流区的指数与THM-TM模型结果对比如图21所示。

图21. 监测点和渗流区指数与THM-TM模型结果对比：(a) 竖向位移，VHMD指数；(b) 竖向达西速度，VDDR指数和VTDR指数；(c) 水平达西速度，HDDR指数和HTDR指数。

在图21(a)中，VHMD指数在大多数测点的量值保持相对稳定，表明各监测点的响应速率相对一致。VHMD指数在拱顶（测点1）、两侧拱肩（测点2和3）、两侧边墙（测点4和5）和两侧拱腰（测点6和7）处呈正值，这些位置发生沉降变形。相反，VHMD指数在底板（测点8）处呈负值，对应隆起变形。SDF与拱顶、拱肩、边墙和拱腰的沉降变形增大相关，同时抑制了底板的隆起变形。

在图21(b)和(c)中，VDDR、VTDR、HDDR和HTDR指数与达西速度的变化趋势一致，指数值与达西速度之间呈正相关关系。热指数（VTDR和HTDR）在对应区域的量值一般大于水力指数（VDDR和HDDR）。这些结果表明，SDF和VRT均与达西速度增大相关，呈近似线性趋势。对比分析表明，围岩温度的影响更为显著，而SDF在水-热-力耦合系统中作为达西速度的修正因素发挥作用。

## 6. 结论

深埋构造混杂岩隧道以强烈的岩性非均匀性和复杂的热-水-力相互作用为特征。本研究建立了结合调查隧道的地质剖面、岩性分布和现场工程条件的耦合数值模型，系统研究了岩体非均匀性、渗流、围岩温度、渗流驱动力和原始围岩温度对隧道响应的影响。研究结果揭示了工程尺度上的热-水-力响应模式和影响机制。主要结论如下：

**(1)** 岩体非均匀性对深埋构造混杂岩隧道的应力、变形和渗流行为有显著影响。岩性界面的存在引起应力重分布，最大应力出现在拱顶附近。非均匀围岩还改变了变形模式，增加了底板隆起和边墙径向收敛，同时减小了拱腰处的竖向位移，对竖向变形的影响大于水平变形。此外，岩体非均匀性改变了渗流分布，导致拱顶、拱肩和拱腰等特定区域出现明显的渗流集中现象。

**(2)** 渗流通过产生孔隙压力降低了围岩的有效应力，导致应力量值略微减小。渗流显著改变了竖向位移模式，增加了拱顶、拱肩和边墙的沉降，同时减小了底板和拱腰的隆起，而对水平位移的影响相对有限。相反，渗流对径向收敛的影响总体有限。

**(3)** 围岩温度由于温度依赖的弹性模量退化及热膨胀引起的应力重分布，略微降低了围岩的应力量值。相反，围岩温度通过刚度降低和热膨胀的综合作用显著增加了变形，同时影响竖向和水平位移。此外，围岩温度通过降低水动力黏度显著增大了达西速度，而渗流区的空间分布基本保持不变。

**(4)** 渗流驱动力和原始围岩温度对深埋构造混杂岩隧道的力学和水力响应有显著影响。渗流驱动力与拱顶、拱肩、边墙和拱腰的沉降呈正相关关系，而与底板的隆起呈反相关关系。渗流驱动力和围岩温度均有助于达西速度增大，呈近似线性趋势，在渗流强度较高的区域观察到更强的响应。原始围岩温度对渗流行为的影响总体更为显著，而渗流驱动力在水-热-力耦合系统中作为达西速度的修正因素发挥作用。

总体而言，研究结果在工程尺度上提供了深埋构造混杂岩隧道在热-水-力耦合条件下应力变形响应和渗流特性的认识。所计算的等效应力、位移和达西速度被解释为识别潜在渗流较高区域响应指标的，而非直接用于确定性隧道稳定性评估或突水预测的准则。尽管如此，本研究基于深埋构造混杂岩隧道的代表性地质剖面，采用具有空间变化材料属性的连续介质有限元框架。未直接捕获显式的不连续体控制行为，如界面滑移、张开、块体运动和断裂控制流。未来研究将采用不连续方法对界面单元进行更加详细和显式的表征。

**CRediT作者贡献声明**

Enda Zhang：撰写——初稿、验证、方法学、调查、数据整理、概念化。

Tiancai He：撰写——审阅与编辑、方法学、调查。

Hongbo Li：撰写——审阅与编辑、调查。

Feng Tian：撰写——审阅与编辑、调查。

Yubo Wang：撰写——审阅与编辑、调查。

Zaobao Liu：撰写——审阅与编辑、方法学、调查、资金获取、概念化。

Enda Zhang | Tiancai He | Hongbo Li | Feng Tian | Yubo Wang | Zaobao Liu