0.   引言

开挖过程中深部巷道经历径向卸荷，在无支护的情况下，巷道内壁围岩由三向受力变为二向受力，原有的应力和能量平衡被打破，此时围岩远离平衡状态并发生破坏。塑性区围岩在变形破坏过程中不断接受内部围岩传递的能量，进行自我组织、形成耗散结构，并逐步向围岩内部扩展与传播，直至重新达到平衡状态。此时，围岩中将形成了一个松弛的碎裂带，即围岩松动圈^[1-4]，而后方靠近松动圈的围岩则处于塑性状态，松动圈和松动圈后处于塑性状态的围岩构成了围岩塑性区。塑性区的厚度深刻地影响巷道的支护安全与稳定性，为了确定围岩塑性区的厚度，国内外众多学者开展了相关研究，在原位测试、模型试验、数值模拟和基础理论等方面取得了许多成果。由于之前许多学者没有充分明确松动圈和塑性区的概念，研究时经常不区分松动圈和塑性区。高树棠^[5]提出围岩松动圈内的岩体为多裂隙体，具有一定的残余强度是塑性区的一部分，且松动圈半径R_t与塑性区半径R存在一定关系，即R_t=α R（0 < α < 1）。需要指出的是，本文引言中文献所说的松动圈大多是指围岩的塑性区，但是为了尊重作者，笔者将沿用原文中的叫法。

1974年池田和彦等人采用声波检测技术观察到了巷道围岩松动圈，并给出松动圈厚度的计算公式，首次证明了松动圈的存在。宋宏伟等^[6]在实测数据的基础上，运用数理统计和计算机模拟的方法得到了松动圈厚度的数学表达式。阎凤山^[7]对岩石类相似材料进行了真三轴试验，得到了不同地应力作用下松动圈的厚度。试验表明：随着原岩应力的增大，松动圈厚度越大；支护对松动圈的形成和松动圈厚度影响不大；原岩应力越大，松动圈越厚，支护越困难。王汉鹏等^[8]、李利平等^[9]研制了低强度和弱黏结特性的软弱破碎围岩的相似材料，在大型三维地质力学模型试验台架上模拟分析了不同施工过程中隧道围岩受力和变形的三维空间演化规律。试验发现：围岩变形具有很强的时空效应；围岩变形的三维扰动深度一般在3倍洞径内。牛双健等^[10]利用大尺度真三轴试验系统分析了开挖过程对围岩松动圈的影响，并着重分析了松动圈内剪应力的演化规律。黄锋等^[11]采用声波探测的方法获取了不同工况条件下巷道的松动圈范围，并结合弹塑性及损伤理论对巷道围岩松动圈半径进行了预测。Park^[12]分别利用Mohr-Coulomb和Drucker-Prager屈服准则，理论求解了围岩软化条件下圆形硐室的塑性区范围。Zareifard等^[13]结合Mohr-Coulomb破坏准则计算了理想弹–脆–塑性圆形巷道围岩的塑性半径。平雯等^[14]运用数值计算软件研究了不同应力条件下巷道围岩松动圈的厚度，并讨论了侧压力系数λ变化对松动圈的影响。当λ < 1时，洞壁松动圈厚度随λ增大而减小；当λ > 1时，洞壁松动圈厚度随λ增大而增大。张頔等^[15]利用3DEC模拟手段，结合实测数据提出了洞室围岩稳定性的评价方法。

岩石作为流变材料其自身的应力状态和变形过程具有很强的时间依赖性，对于深部巷道围岩变形破坏问题，类似的时间依赖特性更加显著。董方庭等^[3]认为围岩松动圈的发展形成伴随着应力调整与重分布，是一个时间过程。根据其现场视察发现（如图 1所示），松动圈的形成通常需要1~3个月甚至更长时间，且巷道收敛变形时间与松动圈的发展时间一致。高树棠^[5]同样发现围岩松动圈的形成具有一定的时间效应，且随时间增长松动圈逐步向围岩内部扩张，如图 2所示。由上述研究可知，深部巷道围岩塑性区的发展和传播需要历经几个月甚至更长时间，在如此长的时间跨度内由开挖引起的动力效应已然消失。而且，从变形波传播的角度分析，根据图 2中数据估算可得塑性区的传播速度约为0.01 m/d，远小于围岩的纵波波速（约2500~4000 m/s）。因此，采用经典弹塑性动力学分析围岩塑性区的演化过程缺乏必要的物理合理性。流变理论的持续发展为描述围岩塑性区的演化规律提供了一个可能的解决方法，但由于三维黏弹塑性流变模型数学表达式复杂且缺乏必要的初边值条件，导致相应控制方程的求解（解析解和数值解）异常困难。另外，黏弹塑性流变模型的建立还严重依赖围岩未知的峰后物理力学性质，涉及复杂的本构关系。因此，通过流变力学等观点定量研究深部围岩塑性区的演化过程理论上可行，但实际执行却困难异常。

图  1  超声波测孔曲线^[3]

Figure  1.  Curves of ultrasonic hole logging^[3]

下载: 全尺寸图片 幻灯片

图  2  围岩塑性区随时间的变化^[5]

Figure  2.  Evolution of plastic zone of rock^[5]

下载: 全尺寸图片 幻灯片

围岩塑性区的演化和传播过程对于巷道支护方案的设计以及围岩的稳定性具有重要意义，而现有的理论研究尚不能很好地定量描述围岩塑性区的演化和传播规律。本文从能量和变形的角度出发，结合统计物理学中的二级相变理论，分析了深部围岩塑性区的演化规律：采用有限差分法求解了围岩塑性区演化的kink波解，分析了各参数对塑性区演化过程的影响规律；通过岩石密度–波速关系曲线以及质量守恒定量，建立了围岩波速与塑性应变之间的定量关系；通过与现场实测数据进行对比，发现kink波解可以很好地描述围岩塑性区的分布及演化规律。

1.   深部围岩塑性区演化的相变模型与kink波解

如图 3所示，岩石变形破坏的典型模式主要有拉伸破坏、剪切破坏、多重剪切破坏和劈裂破坏^[17]。在高围压作用下，深部围岩通常发生多重剪切破坏，其剪应变γ与剪应力τ和体应变${\varepsilon _{\text{v}}}$的关系如图 4所示。

图  3  岩石破坏的基本模式^[16]

Figure  3.  Basic fracturing model for rock samples^[16]

下载: 全尺寸图片 幻灯片

图  4  深部围岩应力–应变曲线

Figure  4.  Stress-strain curves of deep surrounding rock mass

下载: 全尺寸图片 幻灯片

此时，围岩内部将产生多个剪切面并将岩体切割成大小各异的块体，由于块体间摩擦力的作用，围岩并不会完全丧失剪切强度，仍具有一定的承载能力（即残余应力），但块体之间却可发生相对滑动与转动。

众所周知，固体材料与流体材料力学性质的最大区别在于能否保持自身形状，即抵抗剪切变形的能力。在变形破坏过程中，固体材料的抗剪强度降低，抵抗剪切变形的能力减弱，其力学性质逐步向流体过度。观察图 4中剪应力–剪应变（τ – γ）曲线可以发现，深部围岩的变形破坏过程可看作是其力学性质由固体向流体转换的过程。尤其当围岩进入残余强度阶段，围岩中块体在外力作用下发生滑动与转动，表现出了较为显著的流体性质，可类比与散体介质，兼备了固体和流体两种力学性质^[17]。此时，围岩力学性质的转变是一个连续的过程，并未发生突变，具有连续渐进相变的特点，因此可以将深部围岩的变形破坏过程看作是多阶段渐进相变，即二级相变。需要说明的是，本文涉及的连续相变并非通常意义上的固体-液体-气体3种相态之间的转换，而是围岩物理力学性质和状态的变化，类似于物理学中的超导体和磁性转变。

为了方便描述围岩变形破坏过程中应变能H的变化规律，引入如下无量纲塑性剪切应变$\psi$作为序参量

式中${\gamma _{\text{e}}}$为剪应变弹性极限；${\gamma _{\text{c}}}$为剪应变破坏极限；$\gamma$为任意时刻的剪应变，且围岩塑性体积应变$\varepsilon _{\text{v}}^{\text{p}} = \kappa \psi$，$\kappa$为扩容系数，通常由试验确定。

此时，围岩应变能密度H($\psi$)按无量纲塑性剪切应变$\psi$展开^[18]

在高地应力作用下深部围岩内部存在着较大的应力梯度和变形梯度，因此围岩中某一物质点的应力和变形状态受相邻物质点的影响很大。分析开挖过程中围岩能量的变化规律可知，能量由外部围岩以径向压力做功的形式向内部围岩传递^[19]。因此，围岩应变能密度H($\psi$)的表达式应考虑相邻物质点之间的相互作用，可以类比于Ginzburg-Landau展开式^[20]，即

式中，C为系数，且C > 0。

在整个体积Ω中对应变能密度H($\psi$)积分，可得总能量

对式（4）求变分并进行高斯变化，可得

式（5）对时间t求导

对于耗散系统，能量随时间变化应满足

为了保证不等式（7）恒成立，序参量$\psi$随时间t的变化应满足

对于圆形轴对称问题，忽略高阶项，式（8）可化为

式中，$C'$，a，b，c为系数，且$C' = - \mathit{\Gamma} C$，$a = - \mathit{\Gamma} {V_2}$，$b = - \mathit{\Gamma} {V_4}$，$c = - \mathit{\Gamma} {V_6}$。

经观察分析可知，式（9）（或式（8））为对流扩散型偏微分方程，描述了序参量的演化特征。其中，${\nabla ^2}\psi$为“对流”项，描述了物质点之间的长程相互作用；$a\psi + b{\psi ^3} + c{\psi ^5}$为“源”项，表征了物质点自身的变形和含能特征。

式（9）为强非线性偏微分方程，虽然其显式解析解尚未求得，但却可以通过相平面分析法对其解的形式进行预估。利用相平面分析法，可知：当$a > 0$，$b < 0$，$c > 0$时，式（9）对应的解的形式如图 5所示。关于该方程解的其他形式（如周期解和孤立子解）所对应的相图，读者可参看文献[21]。

图  5  相曲线及其对应的解的形式

Figure  5.  Phase curve and form of solution

下载: 全尺寸图片 幻灯片

2.   相变模型求解及参数敏感性分析

本章借助Matlab数学软件对式（9）进行数值求解，采用有限差分法，选取空间步长h=0.25 m，时间步长t =0.01 d，对半径为r₀=5 m的圆形巷道塑性区的演化规律进行分析。此时，方程中各系数取值如下：$C'$=4，a =1，b =$-$15，c=18。

围岩的变形破坏需要经历一个时间过程，因此当t =0时巷道围岩的初始塑性剪切变形为零，即初始条件为$\psi (r,0) = 0$，其物理意义为巷道开挖后t =0时刻塑性剪切应变的分布情况。而边界条件的选取参照文献[22]，可得$\psi (0,t) = 0.2705 \times (1 - {{\text{e}}^{ - {\text{2}}t}})$，其物理意义为巷道壁围岩剪切变形的收敛曲线。

图 6所示为巷道壁围岩序参量的演化曲线，反映了围岩的变形收敛规律，而序参量的幅值则对应于巷道壁围岩的收敛变形。图 7为序参量ψ随时间的变化规律，图中曲线清晰地描绘出序参量ψ以扭结波（kink波）波阵面的形式向围岩内部传播的过程。观察图 8可以发现，序参量ψ的分布可分为3个阶段：I阶段围岩进入残余应力状态，此区域对应于巷道围岩的松动圈。此时，序参量ψ的幅值不随时间变化，这是因为序参量ψ是描述围岩塑性剪切变形的无量纲参数，塑性变形一旦发生便不会自行恢复，如不出现新的能量输入，塑性变形的幅值亦不会增大；Ⅱ阶段内围岩处于塑性状态，相较于Ⅰ阶段围岩的塑性变形较小，且随着距离的增加逐步向弹性阶段过渡。Ⅰ、Ⅱ阶段共同构成了巷道围岩的塑性区；Ⅲ阶段内围岩仍处于弹性阶段，对应于巷道围岩的弹性区。kink波的形状和传播特性与文献[3，5，11]测得的塑性区分布形状类似，为描述围岩塑性区的演化规律提供了可能。

图  6  巷道壁围岩变形收敛曲线

Figure  6.  Convergence of circular excavation

下载: 全尺寸图片 幻灯片

图  7  巷道围岩塑性区演化过程

Figure  7.  Evolution process of plastic zone

下载: 全尺寸图片 幻灯片

图  8  序参量ψ随距离的分布规律

Figure  8.  Distribution of order parameter ψ with distance

下载: 全尺寸图片 幻灯片

图 9为不同时刻下，围岩塑性剪切变形的空间分布情况，从靠近巷道壁开始依次为松动圈、塑性区和弹性区。观察各个图形可以发现，松动圈与塑性区半径随时间逐渐增大，可以很好地描述塑性变形的传播以及松动圈和塑性区的发展演化过程。

图  9  不同时刻下巷道围岩塑性区分布

Figure  9.  Distribution of plastic zone of surrounding rock mass at different time

下载: 全尺寸图片 幻灯片

式（9）中各系数的取值对于kink波传播规律影响显著。为了能够直观的得到各个系数对方程解的影响规律，下面对系数$C'$，a，b，c进行敏感性分析。此处采用控制变量法，每次仅考虑一个系数变化。

（1）系数$C'$敏感性分析

在讨论系数$C'$对方程解的影响时，系数a=1，b=-15，c=18保持不变，分别讨论当$C'$=1，2，3，4时方程解的性质。初边值条件选择如下：

初始条件：$\psi (r,0) = 0$；

边界条件：$\psi (0,t) = 0.2705 \times (1 - {{\text{e}}^{ - 2t}})$；

具体计算结果如图 10所示。

图  10  不同C '值时序参量ψ的演化过程

Figure  10.  Evolution of order parameter ψ with different C '

下载: 全尺寸图片 幻灯片

图 10描绘了不同$C'$值对序参量ψ演化过程的影响规律。系数$C'$取不同数值时，序参量ψ的幅值并没发生变化，这表明系数$C'$对序参量的幅值没有影响；相同时刻下，扭结波（kink波）的波阵面距洞室的距离随着系数$C'$的增加而增大，说明在相同时间内，系数$C'$越大，波阵面传播速度越快、传播距离越远。从系数$C'$对序参量的影响规律可以看出，随着系数$C'$的增大松动圈和塑性区的扩展速度加快。由对流扩散型方程的性质可知，系数$C'$主要用于描述“对流”程度（即物质点之间长程作用）的强弱，因此系数$C'$的取值除了受围岩自身物理力学性质的影响外，还取决于围岩的初始应力状态和应力与变形梯度。从物理学角度，系数$C'$为波阵面的“传播”速度，由量纲分析可知系数$C'$的量纲为m²/s。数值计算过程中速度$C'$的取值并不影响塑性变形的幅值，侧面印证了有限差分法求解的正确性。

（2）系数a敏感性分析

在分析一次项系数a的敏感性时，选取了4组不同的a值：a=1，a=1.2，a=1.5，a=2。此时，保持速度$C'$=1，系数b=-15，c=18不变。初始条件和边界条件如下：

当a=1.0时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.2705 \times (1 - {{\text{e}}^{ - 2t}})$；

当a=1.2时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.2995 \times (1 - {{\text{e}}^{ - 2t}})$；

当a=1.5时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.3408 \times (1 - {{\text{e}}^{ - 2t}})$；

当a=2.0时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.4082 \times (1 - {{\text{e}}^{ - 2t}})$。

由图 11可知，随着系数a的增加，kink波的传播速度增加，序参量的幅值增大。上述现象表明一次项系数a的取值不仅影响序参量幅值，还影响波阵面的传播速度，且波阵面的传播速度和幅值与一次项系数a的取值正相关。

图  11  不同a值时序参量ψ的演化过程

Figure  11.  Evolution of order parameter ψ with different a

下载: 全尺寸图片 幻灯片

（3）系数b敏感性分析

三次项系数b的敏感性分析，选取了4组不同的b值：b=-10，b=-12，b=-15，b=-18。在分析系数b时，保持速度$C'$=1，系数a=1.0，c=18不变。初始条件和边界条件如下：

当b=-10时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.3616 \times (1 - {{\text{e}}^{ - 2t}})$；

当b=-12时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.3125 \times (1 - {{\text{e}}^{ - 2t}})$；

当b=-15时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.2705 \times (1 - {{\text{e}}^{ - 2t}})$；

当b=-18时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.2430 \times (1 - {{\text{e}}^{ - 2t}})$。

计算结果见图 12。

图  12  不同b值时序参量ψ的演化过程

Figure  12.  Evolution of order parameter ψ with different b

下载: 全尺寸图片 幻灯片

分析图 12可知，系数b对于序参量的幅值有较大影响，对于波阵面的传播速度影响很小可忽略。随着系数b代数值的增加，序参量的幅值增大，波阵面的传播速度几乎不变。

（4）系数c敏感性分析

五次项系数c的敏感性，选取了4组不同的c值：c=10，c=18，c=25，c=30。在分析系数c时，保持速度$C'$=1，系数a=1.0，b=-15不变。初始条件和边界条件如下：

当c=10时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.2644 \times (1 - {{\text{e}}^{ - {\text{2}}t}})$；

当c=18时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.2705 \times (1 - {{\text{e}}^{ - {\text{2}}t}})$；

当c=25时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.2764 \times (1 - {{\text{e}}^{ - 2t}})$；

当c=30时，初始条件：$\psi (r,0) = 0$，边界条件：$\psi (0,t) = 0.2815 \times (1 - {{\text{e}}^{ - 2t}})$。

计算结果见图 13。

图  13  不同c值时序参量ψ的演化过程

Figure  13.  Evolution of order parameter ψ with different c

下载: 全尺寸图片 幻灯片

系数c对于序参量的幅值有较明显的影响，但是对于波阵面的传播速度几乎没有影响。且随着系数c代数值的增加，序参量的幅值增大。系数c对方程解的影响规律与系数b相似，但程度较弱。

3.   围岩塑性区力学状态与纵波波速的关系

岩石密度是围岩基本物理力学性质参数之一，其与围岩纵波波速的定量关系为围岩塑性区的探测提供了重要的依据。通过岩石密度–波速关系曲线，不仅可以确定围岩塑性区的位置和厚度，还可以反演出围岩的变形与应力状态。本节主要就如何通过纵波波速反演围岩变形进行介绍。

Birch^[23]对深部火成岩进行了试验研究，发现岩石的密度和纵波波速呈现出线性关系，即C_p=αρ+β。其中，C_p为岩石纵波波速，ρ为岩石密度，α，β为常数，主要由岩石自身性质决定。于是，可得围岩破坏前后纵波波速与对应密度之间的关系

式中$C_{\text{p}}^0$，$C_{\text{p}}^*$分别为围岩破坏前后的纵波波速；${\rho _0}$，${\rho _*}$则为围岩破坏前后对应的密度。此时

根据质量守恒定律，可知

式（14）化简可得

式中，${V_0}$，${V_*}$分别为围岩破坏前后体积，$\Delta V$为围岩破坏前后的体积增量，${\varepsilon _{\text{V}}}$为围岩的体积应变。

将式（13）代入式（11）并化简可得围岩纵波波速与体积变形之间的关系

考虑${\varepsilon _{\text{V}}} = \kappa \psi$，式（14）可化为

4.   围岩塑性区演化理论与实测数据对比分析

为了验证本文提出理论的准确性，需要与实测结果进行对比。限于围岩塑性区演化数据，本文仅选取了文献[5]中的实测数据进行对比分析，实测数据详见表 1，表中纵波波速单位为km/s。

表  1  围岩纵波波速实测数据^[5]

Table  1.  In-situ ultrasonic data of plastic zone in surrounding rock mass^[5]

时间	测点1	测点2	测点3	测点4	测点5	测点6
4月25日	3.17	3.55	4.08	3.91	3.94	3.96
4月29日	2.50	3.25	3.68	3.80	3.83	3.88
5月12日	2.51	2.50	2.79	4.04	3.75	3.76

下载: 导出CSV 

| 显示表格

金解放等^[24]通过室内测试发现岩石密度与纵波波速存在如下统计学关系：

故本文在理论计算时，式（10）取$\beta = - 2.33$。

以文献[5]中实测数据（4月29日）为参考，确定式（9）中系数及其对应的初边值条件，具体如下：a =1，b =-15，c=18，kink波传播速度$C'$=0.3，扩容系数$\kappa = 1.15$；初始条件：$\psi (r,0) = 0$；边界条件：$\psi (0,t) =$ $0.2705 \times (1 - {e^{ - 0.6t}})$。采用上述相同参数，对剩余两组实测数据进行理论预测与对比，具体结果如图 14所示。

图  14  围岩塑性区实测数据与理论对比

Figure  14.  Comparison between theoretical prediction and in-situ data

下载: 全尺寸图片 幻灯片

比较图 14中结果可以发现，理论曲线与现场实测数据拟合程度很好。此时，kink波可以很好地描述围岩塑性区的分布及演化规律，验证了理论模型的准确性。但需要指出的是，理论计算结果整体上较实测数据偏小，由对流扩散方程的性质可知，造成该现象的原因可能有二：①对流效果较弱，此时理论上能量的对流速度小于实际传递速度（即$C'$取值小于真实值）；②“源”项描述围岩所包含的能量小于围岩实际存在的能量，导致理论扩散速度小于实际值。

5.   结论

深部围岩的变形破坏可以看作是其力学性质由固体向流体转变的过程，该转变过程是一个连续的过程，具有连续渐进相变的特点。从能量与变形的角度分析了深部围岩的演化规律，具体结论如下：

（1）采用相平面分析法预测了演化方程的kink波解，并分析了其特点与规律，结果表明kink波的形状和传播特性与围岩塑性区分布和演化规律类似，为描述围岩塑性区的演化规律提供了可能。

（2）采用有限差分法求解了围岩塑性变形区演化的数值解，分析了各参数对塑性区演化的影响规律：系数$C'$仅影响围岩塑性区的扩展速度，并不影响围岩塑性变形幅值的大小，其物理意义可理解为能量的“对流”速度；系数a对塑性变形幅值及波阵面的传播速度都有很大影响，且随着a值的增大，围岩塑性变形增大、塑性区传播速度加快；系数b仅影响塑性变形的幅值，对塑性区传播速度基本无影响；系数c对塑性区变形和传播速度的影响规律与系数b相似，但程度较弱。

（3）通过岩石密度–波速关系曲线和质量守恒定律，建立了围岩波速与塑性变形之间的定量关系，可反演出围岩的变形状态。

（4）通过比较理论计算与现场实测数据可知，kink波可以很好地描述围岩塑性区的分布及演化规律，验证了理论模型的准确性。