在裂隙岩体中，隧道开挖往往会对周围的地质环境和水环境造成显著的扰动。当原本处于平衡态的岩石地层受到隧道开挖的扰动，开挖区域附近会形成开挖扰动区或开挖损伤区^[1-2]。对于结晶岩，开挖扰动区指的是发生了可恢复变形的区域，而开挖损伤区指的是发生了不可恢复的变形区域^[3-4]。除了对地层和应力场的扰动，开挖对地下水分布以及水压力的大小也会造成影响^[5]。岩体中天然裂隙网络的赋存特征和新裂纹的出现是控制开挖损伤区力学表现和水力连通性的重要因素^[6-7]，其对裂隙扩展的研究和模拟以及许多以岩体为主体介质的地下工程，例如核废料存储、深层隧道等，具有显著影响。

目前针对完整岩体中隧道开挖损伤发展过程的研究较充分^[8-12]，但是针对含有天然裂隙网络的裂隙岩体中隧道开挖诱发损伤及其发展规律的研究较少。文献[13]使用RFPA数值方法对含有天然裂隙的岩体开挖进行了初步模拟，指出天然裂隙在损伤发展中起到控制作用。文献[14]通过数值模拟揭示了裂隙网络几何特征对隧道开挖诱发损伤特性的影响规律。除物理场的扰动外，文献[15-18]还对裂隙网络的渗流特性进行了研究。

流固耦合在裂隙岩体开挖中控制着许多重要的物理过程。首先，地下水的压力变化会影响岩体的受力和变形。Biot建立了孔隙弹性理论，对泰沙基有效应力理论进行了拓展，并建立了孔隙水压力与岩体弹性变形之间的联系，提出了比奥固结系数的概念。此外，岩体的应力应变分布也会影响地下水的流动。岩体自身和赋存裂隙在不同应力水平下具有不同的渗透率，对地下水流动有着重要影响。流固耦合可以分为直接耦合和间接耦合两种形式^[19]，其中，直接耦合指的是孔隙介质的变形可以影响压力以及渗流，同时水压力的变化会影响孔隙介质的变形；相反，间接耦合指的是应力和应变的变化会影响到岩体的渗流参数，进而对影响应力场产生影响，同时应力场的变化还会影响到应力的分布并作用于裂隙岩体，是参数到参数的耦合。间接耦合常在渗透性较低的岩层中起到重要作用。

裂隙岩体中隧道开挖会诱发强烈的流固耦合效应，从而引发显著的岩体变形和水压力扰动^[19]。开挖对岩层会造成扰动，对地下水也会产生影响，且两者之间有着密切的关联。首先，随着开挖的进行，岩层中应力的变化会导致裂隙出现开度和联通程度的变化，同时在隧道周边形成损伤，这些物理过程可能会增大隧道周边岩体的渗透性^[20]；其次，开挖形成的隧道在地层中以自由面的形式存在，引导地下水向隧道内流动从而导致隧道围岩水压力下降^[21]。水压力的变化将进一步影响岩体骨架和裂隙上的有效应力，引发孔隙弹性效应^[19]。

针对以上问题，本文研究了裂隙岩体中由于隧道开挖引发的流固耦合效应，并探究其对隧道开挖扰动区/损伤区的影响。首先，提出裂隙岩体中隧道开挖流固耦合模型及方法，通过建立固体模型和流固耦合模型两部分模型进行对比，从损伤延展、裂缝变形，以及地下水压力变化等角度论证了流固耦合模拟在富含裂隙系统的地下空间开挖时的重要性；最后，从应力分布、裂隙变形、损伤效应、地下水渗流和声发射强度等方面对模拟结果进行了分析，并探讨了裂隙网络特征参数对损伤效应的影响规律。

1 控制方程和材料本构 1.1 岩石变形与损伤延展

视岩石为固体域，采用各向同性的损伤模型描述岩石材料的峰值后特性，即

$ \boldsymbol{\sigma}^{\prime}=(1-\omega) \boldsymbol{D}: \boldsymbol{\varepsilon} $

(1)

式中：ω为损伤因子；D为线弹性刚度矩阵；ε为岩石的应变；σ ′为有效应力，σ′= σ-αpI，α为比奥固结系数，p为孔隙水压力，I为单位矩阵。

应力应变的正向定义为受压，负向为受拉，如图 1(a)所示。岩石的张拉损伤因子ω_t和受压损伤因子ω_c计算方式如下式和图 1(b)所示。

$ \omega_{\mathrm{t}}=\left\{\begin{array}{cc} 0, & \kappa_{\mathrm{t}} \geqslant \varepsilon_{10} \\ 1-\frac{f_{\mathrm{tr}}}{E \kappa_{\mathrm{t}}}, & \kappa_{\mathrm{t}}＜\varepsilon_{\mathrm{t0}} \end{array}\right. $

(2)

$ \omega_{\mathrm{c}}=\left\{\begin{array}{cl} 0, & \kappa_{\mathrm{c}} \leqslant \varepsilon_{\mathrm{c} 0} \\ 1-\frac{f_{\mathrm{cr}}}{E \kappa_{\mathrm{t}}}, & \kappa_{\mathrm{c}}>\varepsilon_{\mathrm{c} 0} \end{array}\right. $

(3)

式中：ε_t0=-f_t0/E和ε_c0 = f_c0/E分别为线弹性受拉和受压应变极限，如图 1(a)所示, f_t0和f_c0分别为岩石的抗拉和抗压强度；损伤岩石的残余强度通过残余强度比η来定义；f_tr=ηf_t0和f_cr=ηf_c0分别为残余受拉和受压强度；κ_t和κ_c分别为受拉和受压的历史最大应变。

1.2 裂隙的法向切向变形

天然裂隙的法向开度在不同的荷载条件下会呈现出较大的差异，这种开度-应力关系可以用指数函数表征(图 1(c))为

$ b_{\mathrm{n}}=b_{\mathrm{r}}+\left(b_0-b_{\mathrm{r}}\right) \exp \left(-\xi \sigma_{\mathrm{n}}^{\prime}\right) $

(4)

式中：b_n为受到法向应力情况下的开度；b₀为无应力状态下的开度；b_r为残余开度；σ′ _n为裂隙上的有效法向应力，σ′ _n=σ_n-p，其中σ_n为法向总应力，p为水压力；ξ为应力-开度相关系数，其大小与裂隙的初始法向刚度有关，ξ=1/[K_n0(b₀-b_r)]，其中K_n0为当裂隙的初始法向刚度。法向受压情况下，法向刚度表现出非线性变化特征(图 1(c)):

$ K_{\mathrm{n}}=-\frac{\partial \sigma_{\mathrm{n}}^{\prime}}{\partial b_{\mathrm{n}}}=\frac{b_0-b_{\mathrm{r}}}{b_{\mathrm{n}}-b_{\mathrm{r}}} K_{\mathrm{n} 0} $

(5)

天然裂隙的剪切变形基于库伦准则定义(图 1(d))为

$ \tau_{\mathrm{s}}=\left\{\begin{array}{cl} K_{\mathrm{s}} u_{\mathrm{s}}, & u_{\mathrm{s}}＜u_{\mathrm{p}} \\ \tau_{\mathrm{p}}, & u_{\mathrm{s}} \geqslant u_{\mathrm{p}} \end{array}\right. $

(6)

式中：τ_s为剪切应力；u_s为剪切位移；K_s为裂隙的剪切刚度；τ_p=σ′ _ntan φ_f为峰值剪切应力，φ_f为摩擦角，其对应的应变u_p=τ_p/K_s为峰值剪切应变，当裂隙的剪切位移超越这一阈值，认为裂隙开始滑移。除了切向位移，裂隙的剪切应力同样能带来法向开度的变化，其关系遵循下式(图 1(d)):

$ \mathrm{d} v_{\mathrm{s}}=\left\{\begin{array}{cc} -\tan \varphi_{\mathrm{d}} \mathrm{d} u_{\mathrm{s}}, & u_{\mathrm{p}} \leqslant u_{\mathrm{s}} \leqslant u_{\mathrm{r}} \\ 0, & \text { 其他 } \end{array}\right. $

(7)

式中：φ_d为剪胀角，u_r为残余剪切位移。结合裂隙法向与切向应力带来的开度变化，天然裂隙的总体开度可以表征为

$ b_{\mathrm{f}}=b_{\mathrm{n}}-v_{\mathrm{s}} $

(8)

1.3 裂隙多孔介质中的渗流

裂隙多孔岩体中的天然裂隙流量与裂隙开度的立方成正比，因此, 渗透率与裂隙开度平方成正比，其渗透率公式可表示为

$ k_{\mathrm{f}}=\frac{b_{\mathrm{f}}^2}{12} $

(9)

岩体材料的孔隙度与有效应力假定呈指数关系：

$ \phi_{\mathrm{m}}=\phi_{\mathrm{r}}+\left(\phi_0-\phi_{\mathrm{r}}\right) \exp \left(-\zeta \bar{\sigma}^{\prime}\right) $

(10)

式中：ϕ₀为岩体骨架的初始孔隙度，ϕ_r为残余孔隙度，σ′为有效体积应力，ζ为孔隙度-应力相关系数，其大小与岩体材料的损伤程度有关。对于未损伤的岩体，其值为c_m/(f₀-f_r)；对于已损伤的岩体单元，其值为1/[K_n(b₀-b_r)]。岩体骨架渗透率可表示为

$ k_{\mathrm{m}}=k_0\left(\frac{\phi_{\mathrm{m}}}{\phi_0}\right)^3 \exp (\zeta \omega) $

(11)

2 隧道开挖流-固耦合模型建立

全耦合的流-固有限元模型是本文基于COMSOL多物理场模拟平台完成的，该模型可以捕捉天然裂隙的刚度和变形随应力的变化、岩体孔隙的弹性变形、完整岩体的损伤延展以及岩体骨架和裂隙中的地下水渗流等。该耦合系统的控制方程通过COMSOL内置的瞬态广义-α解法实现。模型中裂隙岩体区域通过非结构化的三角形单元进行离散，其中裂隙采用节理单元进行模拟^[22-23]。

为了探究裂隙岩体中隧道开挖流固耦合效应及其对隧道周围裂隙岩体的损伤影响，本文首先针对比奥固结系数和裂缝网络的影响分别进行了全过程分析；然后，针对不同参数下的裂隙网络进行开挖模拟分析，获得了裂隙岩体中不同裂隙网络对开挖损伤演化过程的影响。

2.1 裂隙岩体属性

岩体和裂隙的材料属性基于Sellafield的BVC火成岩进行选取^[24]，其中，完整岩石的材料参数：密度ρ=2 750 kg/m³，弹性模量E=84.0 GPa，泊松比υ=0.24，抗拉强度f_t=11.6 MPa，抗压强度f_c=250 MPa，渗透率k₀=1×10^-19 m²，孔隙比φ₀=0.01，比奥固结系数α为0.5。裂隙的材料参数：初始法向刚度k_n0=50 GPa/m，初始切向刚度k_s0=50 GPa/m，摩擦角φ_f=24.9°，剪胀角φ_d=5°，初始开度h₀=3×10^-5 m，最大开度h_m=5×10^-5 m，残余开度h_r=5×10^-6 m，残余切向位移u_r=3×10^-3 m。

2.2 裂隙网络

针对隧道开挖流固耦合分析，裂隙网络基于Sellafield的勘探数据生成^[24]。据观测，此处裂隙网络主要包含4组裂隙，其倾向和倾角分别为8°/145°，88°/148°，76°/21°和69°/87°。裂隙的迹线长度介于0.5~250 m之间，迹线长度的分布服从幂律法则。根据上述数据和数值方法生成的一组离散裂隙网络(下文中简称为DFN1)，其将作为后文中隧道开挖的承载介质。

为探究裂缝网络分布特征对开挖诱发损伤的影响，需采用不同分布参数的裂隙网络。针对裂缝长度特征，其服从幂律分布：

$ n(l)=\alpha_{\mathrm{f}} l^{-\alpha}, l \subseteq\left[l_{\text {min }}, l_{\max }\right] $

(12)

式中: n(l)为长度为l的裂缝的条数，α为幂律长度指数，α_f为密度项。

裂缝的长度分布主要由长度指数α控制，其大小一般在1.5~3之间。如果长度指数较小，裂隙网络主要由长裂缝控制；若长度指数较大，裂隙网络内的裂缝长度均较小。针对裂缝数量特征，采用裂隙强度γ表征裂隙网络中的裂隙数量，其计算式^[25-26]为

$ \gamma=\frac{1}{L^2} \int n(l) l^{\prime} \mathrm{d} l $

(13)

式中：l′代表裂隙在模拟区域中的长度，L为模拟区域的长度。

由式(12)、(13)，首先生成4个裂隙网络，其长度指数分别为1.5、2.0、2.5和3.0，然后对每个长度指数均生成裂隙强度分别为1.0、2.0和3.0 m^-1的3个裂隙网络，共计12个裂隙网络。

2.3 边界条件

模拟隧道直径d为2 m，模型尺寸为20 m×20 m，基于隧道开挖影响范围的经典解析解和国内外学者的研究结果，为确保同时满足边界条件和计算效率的需要，本研究隧道边界到隧道中心的距离为隧道尺寸的5倍，可保证边界条件对开挖的影响较小。模拟区域通过三角形网格进行划分，其平均尺寸大约为10 cm。模拟区域中的裂隙视为初始存在的弹性薄层。模型底部的边界条件为滑动约束，在两侧边界和上边界分别施加均布压力，其大小为Sellafield 500 m深度处的总应力水平(竖直方向20 MPa，水平方向10 MPa)，同时在四周边界上施加静水压力(5 MPa)。模拟过程共分为4个阶段：1)将边界荷载和水压力以及边界渗流条件施加到模型边界上；2)维持当前的地应力和水压力，直至应力场和渗流场达到平衡；3)将隧道开挖区域内岩体材料的弹性模量以及开挖边界上的水压力逐步降低到0，以模拟隧道开挖扰动和排水过程；4)继续维持当前的应力和水压力边界条件，直至模型达到新的平衡。

3 结果分析 3.1 固体力学模型与流固耦合开挖模型的对比实验

为了证明流固耦合模拟在开挖损伤评估中的重要性，本文针对同一裂隙网络，在不同的比奥固结系数下分别建立的固体力学和流固耦合模型，然后分别进行了隧道开挖模拟。此节模拟了比奥固结系数分别为0.5和1.0两种情况，以研究不同程度的耦合对于开挖损伤的影响。针对固体力学模型，假设在开挖过程中水压力恒定，因此, 施加的边界条件为当前比奥固结系数下的有效应力，流固耦合模型为水压力与总应力边界条件。模拟开挖时间为0.1 h，同时模拟了开挖后20 h内的损伤延展以展示隧道开挖不同阶段的损伤变化。

图 2为开挖损伤的分布。相比于固体力学模型，流固耦合模型表现出更剧烈的损伤。对于固体力学模型，损伤主要发生在隧道周边区域，然而流固耦合模型中，可以出现在距离隧道周边2~3倍直径的区域外。当比奥固结系数较小时(α=0.50)，固体力学模型与流固耦合模型间的差异较小，然而随着比奥固结系数的增加，两类模型之间的损伤差异更加显著。当比奥固结系数为1.0时，在开挖区域右上方出现了显著且延展距离较大的损伤带，这表明随着流固耦合强度的增大(岩体比奥固结系数的增大)，开挖损伤也会增强，开挖影响区域相应也会增大。在高比奥固结系数情形下，流固耦合模型中的损伤面积能够达到固体力学模型中的3倍及以上。

在研究开挖损伤的基础上，基于所建立的模型获得了隧道开挖引发的天然裂隙变形。图 3为不同阶段的裂缝剪切变形。与固体力学模型相比，流固耦合模型模拟的裂隙滑移更剧烈，且随着α的增大，模型间的差异更加明显。对于固体力学模型，裂隙滑移主要发生在隧道周边区域，然而流固耦合模型中，滑移影响区域更大。

流固耦合模型模拟的开挖中水压力分布如图 4所示。由于开挖卸荷导致的应力突变，隧道周边区域会发生强烈的压应力变化以及强烈的孔隙弹性效应。水压力升高区出现在隧道的上方和下方，降低区出现在隧道的左侧和右侧，这与该应力边界条件下传统解析解应力压缩和拉伸应力区非常吻合。随着流固耦合程度的上升(α的增大)，水压扰动区的面积和压力峰值增大。此外，高压区的分布与剪切位移的分布模式相似(图 3)。由于裂缝的渗透性显著高于岩体骨架，在开挖阶段，裂缝相比于岩体骨架，水压下降速度快得多，t=0.1 h时，在低压区中可以明显观测到裂隙网络。在排水阶段，压力场逐渐下降，水逐渐从基质扩散到裂缝中，然后流入隧道。在1 h内，高压区面积明显缩小，而当t=20 h时，整个水压力场主要以开挖排水导致的低压为主。

3.2 裂隙网络对损伤效应的影响

为了减少单一裂隙网络的特殊性对分析结果的影响，根据式(12)、(13)生成了12个参数化裂隙网络(图 5)。如图 5所示，在隧道周围的损伤区内，由于应力和水力场的变化，岩体的损伤扩展显著。在距离隧道较远的岩体中，损伤一般以雁裂纹形式存在，这进一步增大了裂隙网络的连通性。因此，隧道开挖扰动会激活原有的裂隙，从而在远场也引发扰动。对比同一长度指数下不同裂隙密度的裂隙网络，在原有的裂隙网络中引入新的裂隙，会在裂隙延展的形式上引发变化，从而改变整体损伤的模式。而在同一裂隙密度下，随着长度指数的不断增大，裂隙网络系统从由大裂隙控制的裂隙网络转变为充满小裂隙的系统。在此过程中，随着裂隙长度减小，整体损伤面积增大，同时雁裂纹的延展范围也变广，使得距离隧道周围更远区域的岩体也受到开挖和排水的扰动而产生损伤。图 6为开挖结束时天然裂缝的剪切位移分布。对于以小裂缝为主的裂缝系统，裂缝剪切位移模式与岩石损伤模式高度一致。例如，对于α=3.0和γ=3.0 m^-1的情况，可以观测到一条宏观剪切带斜穿过模拟区域，并且形成了广泛的的岩石损伤和裂缝剪切。对于以大裂缝为主的系统(即α值较小的裂隙系统)，沿一些特定裂缝可以观测到大位移，但此类裂缝滑移几乎没有带来损坏，仅在隧道周边发生了集中开挖导致的损伤。

为了减少特定裂隙网络对于结果的影响，在同样的参数下生成了另外10组裂隙网络，对每组12个裂隙进行流固耦合效应开挖分析，以期定量分析隧道开挖造成的损伤随开挖过程的变化。图 7为在不同参数的裂隙网络下，岩体骨架的损伤和裂隙的滑移随时间的变化。其中，R_d为岩体中损伤面积占模型总面积的百分数，R_s为岩体中裂隙滑移总面积占模型总面积的百分数，折线所示大小为上述各组模型的平均值，并在图中以误差棒标记了各组模型的标准差。在各组模型中，可以发现在开挖阶段的末期(0.1 h左右)，损伤的延展和裂隙的滑移都呈现出快速增长的趋势，其主导因素是开挖引发的强烈应力场扰动和孔隙弹性效应。在排水过程中(0.1~20 h)，依旧能发现一定程度的扰动，尤其在α和γ较大的组，即小裂隙控制的高密度组其表现更为显著。对于高密度组，损伤的量级随α的变化更明显，同时各组模型间的标准差也随着α和γ的增大不断增大。

4 讨论

本文的研究结果与许多现场实测数据吻合，体现了较高的一致性。例如，本文中的隧道渗流量为10^-7 m³/s，这与瑞士Alps附近23个隧道的实测数据相符(10^-8~10^-4 m³/s)^[27]。此外，本文计算的声发射量级总体上也与加拿大Mine-by地下实验室的数据符合。根据现场实测数据，其声发射的量级为-5~-2^[28]，本文的模拟结果与其位于相同量级。

通过固体力学模型和流固耦合模型的对比分析，本文模拟结果揭示了开挖损伤预测中流固耦合分析的重要性。本文通过考虑超孔压的时空变化，可观测到在开挖阶段，较大的比奥固结系数往往会引发更大的压力变化，这是应变突变和岩体基质低渗透率共同引发的结果。这一发现与瑞士Grimsel的FEBEX隧道挖掘期间的现场监测数据以及相关的数值模拟结果^[29]高度一致。在流固耦合模型中，高压区的分布与发生裂缝滑移和基质损伤的位置非常吻合，这也支持了流固耦合模拟的重要性。在流固耦合模型中，本文注意到，在开挖过程里，排水性能差的岩石中出现的应力集中，触发了裂缝的大规模剪切位移，这与储层地质中加压引起的裂隙滑移机理一致^[30]，而裂缝滑动的出现则进一步导致了岩体损伤延展^{[6, 14]}。固体力学模型和流固耦合模型的显著差异表明，评估开挖引起的裂隙岩体压力时，孔隙压力的时空变化至关重要，这对分析系统应力、变形和损伤演化有显著影响。在实际工程中，尤其针对比奥固结系数较大的岩体，考虑到孔隙弹性带来的损伤和裂隙滑移，应适当提高支护等级，同时做好隧道防渗工作。

通过对不同参数控制的裂隙网络的模拟，本文发现在不同长度指数的裂缝网络，在固体力学场和水文力场方面表现出不同的扰动模式。对于以大裂缝为主的裂缝岩体，往往会发生天然裂缝的剪切滑移，而岩石损伤较少；应力扰动模式更多地受天然裂缝位移的影响，而不是完整岩石的破坏；此外，由于裂缝网络的高度连通性，开挖引起的超孔隙水压力变化较快消散。对于以小裂缝为主的裂缝网络，岩石基质的损伤延展和天然裂缝的摩擦滑动均会被开挖触发。因此在实际隧道开挖中，提前检测和识别裂隙的密度以及种类对于后续开挖具有重要作用，尤其是在小裂隙和大密度情况下，针对可能出现的剪切损伤带以及后续渗流区的拓展，应做好隧道排水和防渗工作。

本文在二维情形下进行了隧道开挖的模拟，尽管在实际工程中，往往会受到更复杂的应力状况，例如损伤演化和裂隙联通情况的约束，但本文模拟的渗流和声发射等结果已与现场检测数据达到了定性层面较好吻合，因此, 本文建立并开发的流固耦合模型和计算方法能够有效模拟含裂隙岩体中隧道开挖和后续排水降压过程、后续损伤发展、裂隙滑移和地下水流动过程，可再现裂隙在不同应力状况下的开度和剪切变形。与传统的隧道开挖固体力学模型相比，本文的流固耦合模型考虑了直接和间接耦合机制，并可捕获瞬态应力/压力变化，因此, 在一定程度上能提供更真实的扰动表现。本文的研究结果还强调了详细表征地质介质中裂缝网络特性对于地下工程设计和分析的重要性。

5 结论

1) 隧道开挖带来的流固耦合和孔隙弹性效应会引发局部应力、裂隙变形和水压力的出现，同时流固耦合模型与固体模型相比，损伤面积显著增大。隧道开挖中流固耦合作用引发的损伤在高比奥固结系数裂隙岩体的总体开挖损伤中占比能够达到60%或以上。

2) 不同裂隙网络的分布特征在隧道开挖引发的损伤发展和后续渗流过程中起主导作用，进而影响隧道开挖扰动区和损伤区的形成与发展。大裂缝控制的裂缝网络中以剪切滑移和局部损伤为主，而小裂缝控制时岩体损伤分布更广泛。