黏弹性围岩隧道纵向变形曲线与虚拟支护力解析

Analysis of Longitudinal Displacement Profile and Fictitious Support Pressure of Tunnels in Viscoelastic Rocks

摘要

黏弹性围岩隧道支护设计需同时考虑开挖面空间效应与围岩自身流变效应复杂叠加的影响。为此，本文从隧道施工的3维时空效应出发，通过大量数值模拟揭示黏弹性围岩隧道纵向变形时空演化特性，基于响应面法拟合位移释放系数回归模型。以此为基础，进一步探究黏弹性围岩的虚拟支护力时空演化规律，推导黏弹性围岩的虚拟支护力表达式。研究表明：1）在隧道纵向上，黏弹性围岩位移释放率沿隧道纵向发展可划分为未扰动状态、时空效应共同影响和时间效应控制阶段，与弹性围岩相比在空间效应影响基础上耦合了时间效应，其演化更为复杂；2）开挖面前方位移释放率波动较小，可采用统一表达式进行描述，而后方随开挖速度变化较大，为此将其分为时间控制型、空间控制型和时空效应共同影响型3种，并分别给出相应的拟合公式；3）黏弹性围岩虚拟支护力受时空效应共同影响，不能参照弹性围岩通过位移释放系数翻转直接获得，而需通过围岩特性曲线耦合求解；4）不同围岩参数对应的虚拟支护力分布较为集中，且主要受开挖速度影响，据此给出虚拟支护力的简化拟合公式以供初步设计时使用。研究结果可为流变软岩隧道的支护设计提供理论基础。

Abstract

Objective The design of support structures for tunnels in viscoelastic rocks requires consideration of the complex, superimposed spatial influence of the tunnel face and the rheological properties of the surrounding rock. This study investigates the spatiotemporal evolution of the longitudinal displacement of tunnels in viscoelastic surrounding rocks through extensive numerical simulations and establishes a regression model for the displacement release coefficient using the response surface method. This research provides a theoretical foundation for the design of tunnel support systems in rheological soft rock. Methods It examines the space-time evolution behavior of the fictitious support pressure in viscoelastic surrounding rock and derives an analytical expression for this pressure. A numerical simulation model of the longitudinal deformation curve of viscoelastic surrounding rock is developed, revealing the longitudinal displacement profile (LDP) characteristics of the viscoelastic constitutive model. An empirical formula for the displacement release coefficient is derived based on the spatiotemporal dynamic evolution mechanism. Results and Discussions The study analyzes the effects of varying surrounding rock parameters on LDP, noting that the LDP under viscoelastic conditions demonstrates similarities to, but apparent differences from, those under elastic and elastoplastic conditions. Specifically, the undisturbed region of surrounding rock beyond a certain distance from the unexcavated face shows a displacement release rate of zero. Near the excavation face, the displacement release rate of the surrounding rock gradually increases and then rises sharply. This rate is affected by both temporal and spatial factors, which dynamically evolve depending on the interaction between the excavation speed and the time-dependent characteristics of the surrounding rock. Beyond a certain distance after tunnel face excavation, the displacement release rate is controlled solely by the time effect, making the spatial effect negligible. Viscoelastic surrounding rocks exhibit a superimposition of time effects on the spatial effect range, defined as space-time effects. The region where the spatial effect disappears is influenced only by the time effect, corresponding to the traditional rheological stage. The displacement release rate of viscoelastic surrounding rock within the space-time effect zone correlates with the excavation speed. An empirical formula applicable to viscoelastic surrounding rock for determining the displacement release coefficient is derived utilizing the dynamic interaction between time and space effects. This coefficient, when reduced to a purely spatial effect, aligns with the existing elastic displacement release coefficient and extends its applicability to unexcavated surrounding rock. This study also investigates the space-time evolution behavior of the fictitious support force in viscoelastic surrounding rock and develops an empirical formula for this force. The fictitious support force for the viscoelastic surrounding rock is determined based on the elastic surrounding rock solution using the Laplace transformation. The longitudinal deformation curve of the surrounding rock, as established in this study, facilitates the calculation of the fictitious support force. Conclusions Analysis of various parameters indicates that the variation in the fictitious support force is related to the tunnel excavation speed. When the excavation speed is sufficiently slow, the viscoelastic surrounding rock behaves similarly to elastic rock, and the influence of the time effect on the fictitious support force becomes minimal, being primarily governed by the spatial effect. In contrast, at higher excavation speeds where the time effect becomes significant, the fictitious support force and the displacement release coefficient for viscoelastic surrounding rock cannot be derived directly through inversion; they must be coupled with the surrounding rock characteristic curve, thus being influenced by both time and spatial effects. The empirical formula for the fictitious support force, considering space-time effects, demonstrates that for different surrounding rock parameters, the fictitious support force in viscoelastic surrounding rock is more sensitive to excavation speed.

Graphical abstract

关键词

隧道工程 / 黏弹性围岩 / 围岩变形 / 纵向变形曲线 / 虚拟支护力 / 时空效应

Key words

tunnel engineering / viscoelastic rock / surrounding rock deformation / longitudinal displacement profile / fictitious support pressure / space-time effect

本刊网刊

随着国家战略向中西部地区调整，中国西部交通基础设施得到了大规模发展。而中西部地质条件复杂多变，隧道工程多处于软弱破碎围岩中，表现为变形大且流变性显著，给支护设计带来极大挑战。在进行隧道支护初步设计时，收敛约束法因其概念清晰且能描述支护‒围岩关系而广为流行^[1‒3]。作为收敛约束法的重要组成部分，纵向变形曲线是定量描述隧道开挖时空效应的重要手段，也是支护与围岩相互作用全过程模拟的依据^[4‒6]。因此，准确评估隧道围岩的3维时空效应演化机制，对于收敛约束法的可靠使用至关重要。

现有研究表明，开挖面后方围岩应力转移和位移释放并非瞬时完成，而是随隧道施工推进逐渐达到稳定，即表现为开挖面空间效应^[7‒9]。国内外众多学者对空间效应及其定量表征方法进行了研究，并得到了诸多经典的纵向变形曲线拟合公式。针对弹性围岩，Panet等^[10]基于3维有限元数值模拟软件拟合了位移释放系数公式并进行了修正，Liang等^[11]综合考虑多因素影响，提出了扩展收敛函数进行修正；针对弹塑性围岩，Chern等^[12]根据地下水电站现场实测数据，采用最优拟合方式得到了位移释放系数拟合函数。Vlachopoulos等^[13]研究了开挖速度对塑性区发展的影响，并得到了基于最大归一化塑性半径的纵向变形曲线V‒D公式，Alejano等^[14]将其进一步推广应用至应变软化围岩。Wang等^[15]在V‒D基础上考虑开挖面加固措施对位移释放系数的影响，提出了基于开挖面挤压变形的纵向变形曲线。Luo等^[16]采用黄金分割法对围岩参数进行反分析，并结合室内试验对位移释放系数进行了拟合。上述研究对于开挖面空间效应的认识具有积极的推进作用，但均未考虑岩体本身的流变性。事实上，实际工程中软岩隧道即使支护过后仍需要较长时间达到稳定，当施工停滞时围岩变形仍在发生，无法用空间效应解释。因此，软岩隧道施工过程受时间效应与空间效应的叠加作用，无法直接使用既有纵向变形曲线进行支护设计。

对于流变性软岩的变形特性，现有研究主要采用理论^[17‒21]或现场实测方法^[22]进行分析。理论分析时，隧道3维开挖效应的描述实则来源于弹性或弹塑性空间效应研究结果，显然未考虑流变性软岩的时效性。由于隧道施工过程中必然要施加支护，现场实测对应的分析结果实则为支护与围岩共同作用下的隧道变形特征。因此，对流变性软岩在无支护条件下因隧道开挖而导致岩体自身时空效应的认知不够深入，仍需进一步探索。

此外，为了定量描述开挖面空间约束效应，在进行收敛约束分析时常引入虚拟支护力作为3维问题与2维平面应变问题的桥梁，即：将约束效应等效为作用于隧道边界的支护力，并随着开挖距离而逐渐衰减^[23‒25]。虚拟支护力一般可由纵向变形曲线和围岩特性曲线耦合求解^[26]，而对于流变性软岩，因缺乏准确的纵向变形曲线表征方法，虚拟支护力演化及其影响因素尚不明确，这也使得流变性软岩支护的定量化设计步履维艰。

鉴于此，本文针对软岩隧道施工的3维时空效应，建立黏弹性围岩纵向变形曲线数据库，基于广义Kelvin本构模型揭示黏弹性围岩隧道施工的3维时空效应，拟合黏弹性围岩隧道纵向变形曲线经验公式。在此基础上，通过围岩纵向变形曲线与围岩特性曲线的耦合分析，获得黏弹性围岩虚拟支护力求解方法，并就其时空演化特性和影响因素进行系统分析，最后提出虚拟支护力简化拟合公式，以便于工程人员使用。

1 黏弹性围岩纵向变形时空效应演化特征

1.1 数值模型

数值模拟是计算围岩纵向变形曲线的主要方法。本文采用FLAC3D软件对黏弹性围岩纵向变形曲线进行计算，本构模型采用广义Kelvin模型，其参数包括体积模量K、Maxwell剪切模量G_{\mathrm{0}}、Kelvin剪切模量G_{\mathrm{1}}、Kelvin黏滞系数{\eta }_{\mathrm{1}}。在地层反演中通常使用1维广义Kelvin模型来描述反演结果，其中，模型参数包括串联弹簧的弹性模量E_{\mathrm{0}}、Kelvin弹簧的弹性模量E_{\mathrm{1}}及Kelvin体中黏壶的黏滞系数{\eta }_{{}_{\mathrm{1}}}^{\mathrm{E}}，注意到常体积模量3维假设条件下Kelvin弹簧原件对应泊松比为0.5，因此，上述反演参数与FLAC3D模型参数换算关系为^[27]：

\left\{\begin{array}{l}\begin{array}{l}K=\frac{E_{\mathrm{0}}}{\mathrm{3}(\mathrm{1}-\mathrm{2}\mu )},\\ G_{\mathrm{0}}=\frac{E_{\mathrm{0}}}{\mathrm{2}(\mathrm{1}+\mu )},\end{array}\\ \begin{array}{l}{G}_{\mathrm{1}}=\frac{E_{\mathrm{1}}}{\mathrm{3}},\\ {\eta }_{\mathrm{1}}=\frac{{\eta }_{{}_{\mathrm{1}}}^{\mathrm{E}}}{\mathrm{3}}\end{array}\end{array}\right.（1）

式中，\mu为泊松比。

假设深埋隧道受静水应力作用，初始地应力选取3、5、10、15、20、30 MPa这6个梯度。隧道开挖半径为6 m，模型的X、Z方向取20倍洞径（120 m），Y方向取200 m，研究断面对应为100 m处。模型上侧和右侧设为自由边界，不添加约束，其余面均设置为法向约束边界，地应力加载在自由边界处。数值模型及其边界条件如图1所示。

本文选择分步开挖，并记录每次开挖循环后的研究断面径向位移。随着开挖面不断推进，可获得一系列位移变形与空间位置的对应关系，即具备时空效应的纵向变形曲线。由于目前尚无Kelvin模型取值范围及建议值，根据既有研究^[28]，选取14种不同岩土Kelvin反演结果作为计算参数。

1.2 典型因素对围岩纵向变形曲线影响

1.2.1 地应力对围岩纵向变形曲线影响分析

通过式（1）的换算关系，得到14组实际工程的Kelvin反演参数，如表1所示。

表1中，同一工况的岩体参数一致，各工况根据地应力分为6种情况，如A‒1、A‒2、A‒3、A‒4、A‒5、A‒6的地应力分别为3、5、10、15、20、30 MPa。以工况A为例，不同地应力条件下纵向变形曲线变化如图2所示。图2中，u_{\mathrm{r}}为围岩位移，R为隧道半径，x为离开挖面的距离。由图2可见，在同一Kelvin参数下，地应力提升会导致对应的围岩位移增加。

图2中，3维条件下广义Kelvin模型围岩特性曲线只与时间有关，围岩位移u_{\mathrm{r}}(t)为：

u_{\mathrm{r}}(t)=\frac{{\sigma }^{\mathrm{0}}R}{\mathrm{2}}\left(\frac{\mathrm{1}}{G_{\mathrm{0}}}+\frac{\mathrm{1}}{G_{\mathrm{1}}}\left(\mathrm{1}-{\mathrm{e}}^{-\frac{G_{\mathrm{1}}}{{\eta }_{\mathrm{1}}}t}\right)\right)（2）

式中，{\sigma }^{\mathrm{0}}为原岩应力，t为时间。

当t\to \mathrm{\infty }时，得：

\frac{u_{\mathrm{m}\mathrm{a}\mathrm{x}}}{{\sigma }^{\mathrm{0}}}=\frac{R}{\mathrm{2}}\left(\frac{\mathrm{1}}{G_{\mathrm{0}}}+\frac{\mathrm{1}}{G_{\mathrm{1}}}\right)=\frac{R}{\mathrm{2}}\left(\frac{G_{\mathrm{0}}+G_{\mathrm{1}}}{G_{\mathrm{0}}{G}_{\mathrm{1}}}\right)（3）

可见，黏弹性围岩最大位移与地应力比值只与隧道半径、Maxwell剪切模量、Kelvin剪切模量3个参数有关，这与本文数值模拟研究结论一致。

为了进一步分析地应力对围岩纵向变形曲线影响，定义位移释放系数u^{\mathrm{*}}(x)=u(x)/u_{\mathrm{m}\mathrm{a}\mathrm{x}}，其中，u(x)为x处的围岩位移，u_{\mathrm{m}\mathrm{a}\mathrm{x}}为无支护条件下围岩最大位移。图3为随机选取4组工况所得结果。由图3可见，地应力对位移释放系数无影响。

1.2.2 Kelvin参数对围岩纵向变形曲线影响分析

为进一步探究Kelvin参数对位移释放系数的影响，选取工况A的参数作为标准组，开挖速度设置为2 m/d，围岩体积模量、Kelvin剪切模量、Kelvin黏滞系数、Maxwell剪切模量均额外设计4个水平进行计算，其中，2、3水平对应参数与标准组具备同等差值，1水平对应参数为标准组的1/10，4水平对应参数为标准组的10倍，具体参数见表2。

将改变不同Kelvin参数对应的围岩位移曲线与位移释放系数曲线分别绘制，首先考虑体积模量变化对位移释放系数的影响，如图4所示。对于工况A而言，体积模量的变化对两者的影响都不大，考虑FLAC3D数值模拟计算过程为有限差分计算，可以认为体积模量对位移变形和位移释放系数无影响。

Maxwell剪切模量变化对位移释放的影响如图5所示。由图5可见，Maxwell剪切模量直接影响位移，而对位移释放系数几乎无影响，5种工况所得位移释放系数基本重合，研究断面的变化趋势皆从开挖面前2倍隧道半径R开始变化，到开挖面后6倍隧道半径抵达最大值。

Kelvin剪切模量变化和黏滞系数变化对位移释放系数的影响分别如图6和7所示。由图6和7可见，当Kelvin剪切模量减小，对应的洞周位移增加，位移释放系数几乎不变。然而，Kelvin黏滞系数对两者均可被认为无影响。

1.2.3 开挖速度对围岩纵向变形曲线影响分析

由式（2）可知，3维条件下广义Kelvin模型与时间相关的围岩特性可用延滞时间t_{\mathrm{d}}（t_{\mathrm{d}}={\eta }_{\mathrm{1}}/G_{\mathrm{1}}）表征，将14种不同工况的延滞时间汇总于表3。

在对单一变量Kelvin参数的影响研究中，工况N相较工况A对应的时间效应更为显著，同时在计算模型范围内得出最大收敛位移所需的速度也越来越苛刻。根据延滞时间对不同的Kelvin参数重新排序，并以延滞时间1 d为界，将延滞时间小于1 d的工况列入表4。表4的各工况又根据开挖速度被分为9种情况，如1‒1、1‒2、1‒3、1‒4、1‒5、1‒6、1‒7、1‒8、1‒9的开挖速度分别为16.0、8.0、4.0、2.0、1.3、1.0、0.8、0.4、0.2 m/d。

当延滞时间较小时（<588 s），位移释放系数存在较为相似的性质。以工况1为例，低延滞时间（<588 s）下开挖速度对位移释放系数的影响如图8所示。当开挖速度大于4 m/d时，距开挖面相同空间位置的位移释放系数变化随开挖速度增加而减小；当开挖速度小于4 m/d时其变化可忽略不计。

对于延滞系数大于等于588 s的Kelvin参数，以工况7为例，中延滞时间（0.22～0.93 d）下开挖速度对位移释放系数的影响如图9所示。由图9可以看出，当x/R时，位移释放系数初始值相较于图8中掌子面（0点）前方2倍隧道半径内提前到了3倍隧道半径，而在开挖速度为2 m/d的条件下，已无法使时间效应抵达最大值。在工况7～9中，开挖速度需达1 m/d方能使时间效应达到最大，对应位移释放系数趋于同一值。

在工况7～9中，开挖断面处（x/R=\mathrm{0}）的位移释放系数随时间波动较大。工况7‒1下，开挖面位移释放系数仅为0.11，而工况7‒10下，开挖面的位移释放系数为0.31。当开挖速度较快时，距开挖面相同空间位置处的位移释放系数随着开挖速度的减小而增加，对应隧道达到最终收敛的空间位置也会后延，工况7中，所有工况在模型计算范围内均可得到最大位移值。

为了进一步讨论开挖速度与高延滞时间的相互作用机制，将延滞时间高于1 d的Kelvin参数及对应开挖速度工况分别列入表5、6，高延滞时间（>1.56 d）下开挖速度对位移释放影响如图10所示。由图10可见：在开挖面前方（x/R），位移释放初始值局限于3倍隧道半径内，与低延滞时间围岩性质相似，时间效应释放完毕时，所要求的开挖速度更慢；当x/R>\mathrm{0}时，随着延滞时间的增加，获得正常“S”形曲线需要的开挖速度越来越慢，对于工况14，仅有工况14‒7、14‒8、14‒9可得到完整的“S”形曲线，对应的最快开挖速度为0.4 m/d，隧道每次推进的开挖时长为5 d，低于5 d的单次开挖时长无法得到完整的“S”形曲线。

2 黏弹性围岩位移释放系数的函数表达

2.1 拟合函数基本形式

黏弹性围岩位移释放系数需在空间效应的基础上进一步耦合时间效应，相较于弹性或弹塑性围岩，除与开挖面空间距离有关外，还与自身黏滞系数相关，这导致在模型尺寸内会存在无法达到最大位移的情况。以工况14为例，当速度小于0.8 m/d时才能在100 m计算范围内出现最大位移。基于现有研究^[9]，空间效应于10倍隧道半径后基本消失，故认为10倍隧道半径内存在时空共同影响，而此范围外仅受时间效应影响。

由式（2）、（3）可得，广义Kelvin本构模型中时间效应导致的位移释放系数公式为：

\frac{u_{\mathrm{r}}\left(t\right)}{u_{\mathrm{m}\mathrm{a}\mathrm{x}}}=\mathrm{1}-\frac{G_{\mathrm{0}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}{\mathrm{e}}^{-\frac{G_{\mathrm{1}}}{{\eta }_{\mathrm{1}}}t}（4）

将对应数值模拟结果与解析结果对比分析，10倍隧道半径外时空效应分析如图11所示。由图11可见，在10倍隧道半径外，空间效应的影响可忽略不计，此时可直接采用不计空间效应的时变解析解。需要注意的是，由于前期的时空效应共同影响，时间效应的起始点需通过10倍隧道半径处位移修正。

考虑10倍隧道半径内位移变化，此范围位移释放受空间效应和时间效应共同影响。参考Panet等^[10]研究成果，拟合函数f(x)为：

f(x)=a{\mathrm{e}}^{bx}+c（5）

式中，a、b、c为与时空效应相关的待定系数。

结合具体工况，10倍隧道半径内位移释放系数的拟合结果汇总于表7。由表7可知，拟合优度r²均在0.97以上，式（5）可较好地说明变化趋势，拟合结果有效。

2.2 基于时空效应的位移释放系数经验公式

由于拟合函数式（5）中待定参数与时空效应有关，因此需要分类拟合。分析可知，位移释放系数在研究断面开挖前波动性较小，可直接以计算结果均值进行拟合：

（6）

对于开挖面后方（x>0），在隧道开挖10倍隧道半径内时间效应与空间效应共同影响隧道变形，导致位移释放系数曲线的形式复杂多样。为了更好地分析位移释放系数曲线特征，将时空效应影响效果分为空间效应控制型、时空效应共同影响型和时间效应控制型3种类型。

1）空间效应控制型

将计算获得的所有工况位移释放系数汇总于图12。由图12发现，在工况8‒2上下出现了较为明显的分界线，工况8‒2以上的位移释放系数曲线密集，工况8‒2以下则较为分散。密集区域内的位移释放系数对应不同围岩和开挖速度工况，与空间效应控制型结果相接近，此状态下时间效应影响不显著，对应的位移释放系数有较好的“集中”态势。

由Kelvin参数、开挖速度对位移释放系数影响的研究结果可知：不同Kelvin参数的位移释放系数一定程度上随着开挖速度的降低而变大；由于空间效应影响，位移释放系数存在上限。对应不同的Kelvin参数，存在相应临界开挖速度使位移释放系数靠近工况8‒2所得结果，开挖速度较该速度大则对应的位移释放系数低于工况8‒2，反之则高于工况8‒2。空间效应控制型临界开挖速度见表8。

以表8中临界开挖速度为下界，考虑当黏滞系数为0时，对应最大开挖速度16 m/d已足够满足隧道实际施工速度判别需求，可得位移释放系数是否为空间效应控制型的判据如下：

（7）

式中，v为隧道掘进速度，t_{\mathrm{d}}为延滞时间。

开挖速度满足该判据则对应位移释放系数呈现较好的“集中”趋势，可直接借助包络曲线及对应平均值直接进行预测，空间效应控制型位移释放系数拟合结果如图13所示。

2）时空效应共同影响型

空间效应和时间效应共同点表现为，都会阻止隧道洞周位移；不同点在于，空间效应影响下，位移释放系数存在上限包络线，而时间效应的核心机理为围岩延滞变形在隧道开挖过程中的3维体现，在相同开挖速度下，时间效应随延滞时间增加而效应显著。时空效应的实际效果由时间、空间效应的相对显著关系决定，如隧道开挖初期，由于时间效应对应变形较小，随着隧道掘进，较快的空间效应衰减对应较快的变形速率。依据该机理分析，发现工况10和11出现明显分界。时空效应共同影响对应的位移释放系数如图14所示。由图14可见，工况9‒1、9‒2、10‒1、10‒2、10‒3与工况11‒1、11‒2、11‒3彼此都存在着交叉，呈现先期变形小而后变形速率大的特点。

当延滞时间为0.43 d（工况8）且开挖速度在8.0 m/d内，可认为主要受空间效应影响，已足够满足隧道施工需求。而当延滞系数大于3.36 d时，时间效应为主要影响因素。时空效应共同影响可在空间效应判别基础上进一步借助延滞时间进行判别：

（8）

时空效应共同影响对应的延滞时间分布较为集中，由图14可知，时空效应共同影响对应的计算结果具有较好的空间分布，依据本文计算公式对其余位移释放系数有较好的预测结果。时空效应共同影响型纵向变形曲线待定系数如表9所示，使用插值法确定系数a、b、c的值。

为了验证预测的可靠性，以工况10‒2为例，取工况10‒1、10‒3对应的系数a、b、c进行线性插值，计算可得预测工况10‒2对应的系数a=-0.79、b=-0.35、c=0.94。将预测值与实测结果进行拟合对比，所得r²=0.928 3，说明所提出的拟合公式可以较好地预测位移释放系数，预测结果与实际值对比结果如图15所示。

3）时间效应控制型

分析可知，延滞时间在3.36 d以上且开挖速度较快条件下，时间效应为主要影响因素，延滞时间高于3.36 d的工况下，时间效应控制型纵向变形曲线如图16所示。由图16发现，时间效应控制型位移释放系数与延滞时间有较为明显的相关性。

取延滞时间3.36 d作为下限，判断围岩是否为主要受时间效应影响，具体判据如下：

\left\{\begin{array}{l}v>\mathrm{16}{\mathrm{e}}^{-\mathrm{2.5}{t}_{\mathrm{d}}},\\ t_{\mathrm{d}}>\mathrm{3.36}\end{array}\right.（9）

考虑位移释放系数曲线在时间效应控制型条件下的变化趋势，在该效应影响下曲线发展趋势应主要与围岩自身的延滞时间参数相关。选取16.0 m/d开挖速度为代表，系数b与延滞时间的函数关系，r²=0.970 0，拟合函数如下：

b=-\mathrm{0.439}\mathrm{4}{\mathrm{e}}^{-\mathrm{0.156}\mathrm{5}{t}_{\mathrm{d}}}+\mathrm{0.678}\mathrm{6}（10）

系数a与c不仅受延滞时间的影响，还同时受开挖速度的影响，以该二元变量为自变量，考虑交互项作用，基于二次曲面对待定系数a、c进行拟合，拟合优度分别为0.928 3、0.953 5，具体如下：

a=-\mathrm{0.858}+\mathrm{0.009}{t}_{\mathrm{d}}+\mathrm{0.044}v-\mathrm{0.000}\mathrm{27}v{t}_{\mathrm{d}}-\mathrm{0.001}\mathrm{5}{v}^{\mathrm{2}}（11）

c=\mathrm{1.119}-\mathrm{0.012}{t}_{\mathrm{d}}-\mathrm{0.046}v+\mathrm{0.000}\mathrm{34}v{t}_{\mathrm{d}}+\mathrm{0.001}\mathrm{33}{v}^{\mathrm{2}}（12）

通过式（10）～（12）可求解时间效应控制型纵向变形曲线。需要指出，本文计算过程中延滞时间最大为27.67 d，当实际围岩延滞时间超过此值时，仍需数值模拟手段另行计算。

综上，在开挖面前方，位移释放系数经验公式形式统一，在开挖面后方，依据黏弹性围岩时空协同效果，位移释放系数，采用不同的求解方法。黏弹性围岩位移释放系数经验公式为式（5），弹性围岩位移释放系数经验取值见表10。

2.3 验证与分析

1）与实测数据对比

向家坝水电站位于金沙江水系，是金沙江水电基地的最末一个梯级电站^[29‒30]。水电枢纽的水平埋深为126～371 m，垂直埋深为105～225 m。中导洞开挖等效半径为3.75 m，施工进度为1.7 m/d，采用台阶法施工。向家坝水电站地质的广义Kelvin参数具体为：K=\mathrm{40.243}\mathrm{2}\mathrm{G}\mathrm{P}\mathrm{a}、G_{\mathrm{0}}=\mathrm{8.623}\mathrm{5}\mathrm{G}\mathrm{P}\mathrm{a}、G_{\mathrm{1}}=\mathrm{0.946}\mathrm{5}\mathrm{G}\mathrm{P}\mathrm{a}、{\eta }_{\mathrm{1}}=\mathrm{0.136}\mathrm{8}\mathrm{M}\mathrm{P}\mathrm{a}·\mathrm{h}。由广义Kelvin参数可获得延滞时间t_{\mathrm{d}}=\mathrm{0.144}\mathrm{5}\mathrm{d}。用判据式（7）计算出临界开挖速度为11.15 m/d，即在开挖速度低于11.15 m/d条件下，洞周位移主要受空间效应影响，可直接采用图13的经验表达式，即：

u^{\mathrm{*}}\left(x\right)=\frac{u\left(x\right)}{u_{\mathrm{m}\mathrm{a}\mathrm{x}}}=-\mathrm{0.729}{\mathrm{e}}^{-\mathrm{0.963}\frac{x}{R}} +\mathrm{1}（13）

向家坝中导洞洞周位移预测值与实测值对比情况如图17所示。由图17可见，本文函数的拟合优度为0.924 5，可较好地描述隧道开挖过程中位移释放随开挖面推进的关系。

研究断面在开挖瞬间预测值与实测值相差较大，这是由于实际监测仪器布置的滞后性，只能捕捉到掌子面（0点）后方一定距离的洞周位移，因此，掌子面位移释放27%更符合工程实际。

对比距掌子面0.5 m后，预测值与实际值可知，排除开挖瞬间的位移释放差异，本文预测值略大于实测值，这是由于实际施工采用台阶法开挖，工法本身对位移释放具备一定约束作用，且支护作用下变形必然低于无支护条件。

2）与既有公式对比

分析可知，当开挖速度足够慢时，对应位移释放系数会退化为仅存弹性解特征。将本文计算结果与既有位移释放系数曲线对比，结果如图18所示。由图18可见，本文弹性解与Corbetta等^[31]所得弹性结果更相似，表明本文拟合公式的相对误差较小。既有公式在黏弹性围岩应用时可能导致支护设计偏于危险。因此，本文方法获得的围岩位移释放率更能反映黏弹性围岩隧道的变形规律，可使支护设计更具科学性和可靠性。

3 黏弹性围岩虚拟支护力理论解析

3.1 弹性围岩虚拟支护力解答

弹性围岩隧道围岩位移u_{\mathrm{r}}表达式为^[32]：

u_{\mathrm{r}}=\frac{R}{\mathrm{2}G}\left({\sigma }^{\mathrm{0}}-p_{\mathrm{x}}\right)（14）

式中，G为围岩剪切模量，p_{\mathrm{x}}为虚拟支护力。

则虚拟支护力p_{\mathrm{x}}可表示为：

p_{\mathrm{x}}={\sigma }_{\mathrm{0}}-\frac{\mathrm{2}G{u}_{\mathrm{r}}^{p_{\mathrm{x}}}}{R}（15）

当虚拟支护力完全消失，则对应隧道在无支护条件下的洞周最大位移u_{\mathrm{r},\mathrm{m}\mathrm{a}\mathrm{x}}为：

u_{\mathrm{r},\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{R}{\mathrm{2}G}{\sigma }^{\mathrm{0}}（16）

对于弹性围岩，可根据Panet等^[10]的公式求得对应弹性围岩虚拟支护力为：

p_{\mathrm{x}}={\sigma }^{\mathrm{0}}\left(\mathrm{0.75}{\left(\frac{\mathrm{0.75}}{\mathrm{0.75}+x/R}\right)}^{\mathrm{2}}\right)（17）

由于Panet等^[10]公式仅适用于开挖面后方的位移释放表达，为此结合李煜舲等^[33]实测获得的位移释放系数计算对应虚拟支护力：

p_{\mathrm{x}}={\sigma }^{\mathrm{0}}\left(\frac{\mathrm{1}}{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\frac{\mathrm{1}}{\mathrm{3}}-\frac{x}{D}\right)\right)（18）

将由Panet^[10]和Lee^[33]等的公式拟合数据获得的位移释放系数和虚拟支护力进行对比，弹性围岩归一化虚拟支护力p(t)/{\sigma }^{\mathrm{0}}如图19所示。图19中，p(t)为虚拟支护力与时间的函数。由图19可见，弹性围岩虚拟支护力和位移释放系数数值加和总为1，两者在其交点对称，通过交点翻转后完全统一。

3.2 黏弹性围岩虚拟支护力解答

基于弹性围岩求解结果，对u_{\mathrm{r}}取Laplace变换，得{\tilde{u}}_{(\mathrm{r},s)}为：

{\tilde{u}}_{(\mathrm{r},s)}={\int }_{\mathrm{0}}^{\mathrm{\infty }}\frac{R}{\mathrm{2}G}({\sigma }^{\mathrm{0}}-p(t)){\mathrm{e}}^{-st}\mathrm{d}t（19）

式中，s为拉式空间。

式（19）可改写为：

{\tilde{u}}_{(\mathrm{r},s)}=\frac{{\sigma }^{\mathrm{0}}R}{\mathrm{2}s\tilde{G}\left(s\right)}-\frac{R}{\mathrm{2}\tilde{G}\left(s\right)}{\int }_{\mathrm{0}}^{\infty }p(t){\mathrm{e}}^{-st}\mathrm{d}t（20）

式中，\tilde{G}\left(s\right)为G\left(s\right)的Laplace变换，其中，G\left(s\right)为中间变量，G\left(s\right)=\frac{G_{\mathrm{0}}+G_{\mathrm{1}}+{\eta }_{\mathrm{1}}s}{G_{\mathrm{0}}{G}_{\mathrm{1}}+G_{\mathrm{0}}{\eta }_{\mathrm{1}}s}。

1维广义Kelvin本构方程为：

\sigma +\frac{{\eta }_{\mathrm{1}}^{\mathrm{E}}}{E_{\mathrm{0}}+E_{\mathrm{1}}}·\frac{\mathrm{d}\sigma }{\mathrm{d}t}=\frac{E_{\mathrm{0}}{E}_{\mathrm{1}}}{E_{\mathrm{0}}+E_{\mathrm{1}}}\epsilon +\frac{E_{\mathrm{0}}{\eta }_{\mathrm{1}}^{\mathrm{E}}}{E_{\mathrm{0}}+E_{\mathrm{1}}}·\frac{\mathrm{d}\epsilon }{\mathrm{d}t}（21）

式中，\sigma为应力，\epsilon为应变。

考虑均匀且各向同性岩体，可将\sigma、\epsilon、弹性模量E、{\eta }_{\mathrm{1}}^{\mathrm{E}}经Laplace变换为s_{ij}（i、j分别为方向和作用平面）、e_{ij}、\mathrm{2}G（G为剪切模量）、\mathrm{2}{\eta }_{\mathrm{1}}，对应变换获得3维条件下的广义Kelvin本构方程：

s_{ij}+\frac{{\eta }_{\mathrm{1}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}·\frac{\mathrm{d}{s}_{ij}}{\mathrm{d}t}=\frac{\mathrm{2}{G}_{\mathrm{0}}{G}_{\mathrm{1}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}{e}_{ij}+\frac{\mathrm{2}{\eta }_{\mathrm{1}}{G}_{\mathrm{0}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}·\frac{\mathrm{d}{e}_{ij}}{\mathrm{d}t}（22）

则可得微分算子Q_{\mathrm{d}}\left(D\right)和P_{\mathrm{d}}\left(D\right)分别为：

Q_{\mathrm{d}}\left(D\right)=\frac{G_{\mathrm{0}}{G}_{\mathrm{1}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}+\frac{{\eta }_{\mathrm{1}}{G}_{\mathrm{0}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}·\frac{\mathrm{d}}{\mathrm{d}t}（23）

P_{\mathrm{d}}\left(D\right)=\mathrm{1}+\frac{{\eta }_{\mathrm{1}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}·\frac{\mathrm{d}}{\mathrm{d}t}（24）

式（23）、（24）中，D为实空间。

对Q_{\mathrm{d}}\left(D\right)和P_{\mathrm{d}}\left(D\right)进行Laplace变换，分别得{\tilde{Q}}_{\mathrm{d}}\left(s\right)、{\tilde{P}}_{\mathrm{d}}\left(s\right)：

{\tilde{Q}}_{\mathrm{d}}\left(s\right)=\frac{G_{\mathrm{0}}{G}_{\mathrm{1}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}+\frac{{\eta }_{\mathrm{1}}{G}_{\mathrm{0}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}s（25）

{\tilde{P}}_{\mathrm{d}}\left(s\right)=\mathrm{1}+\frac{{\eta }_{\mathrm{1}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}s（26）

则可得\tilde{G}\left(s\right)为：

\tilde{G}\left(s\right)=\frac{Q{}_{\mathrm{d}}{}^{}(s)}{P{}_{\mathrm{d}}{}^{}(s)}=\frac{G_{\mathrm{0}}{G}_{\mathrm{1}}+G_{\mathrm{0}}{\eta }_{\mathrm{1}}s}{G_{\mathrm{0}}+G_{\mathrm{1}}+{\eta }_{\mathrm{1}}s}（27）

将式（27）代入式（20）可得：

\begin{array}{l}{\tilde{u}}_{(\mathrm{r},s)}=\frac{{\sigma }^{\mathrm{0}}R}{\mathrm{2}s}·\frac{G_{\mathrm{0}}+G_{\mathrm{1}}+{\eta }_{\mathrm{1}}s}{G_{\mathrm{0}}{G}_{\mathrm{1}}+G_{\mathrm{0}}{\eta }_{\mathrm{1}}s}-\frac{R}{\mathrm{2}}·\\ \frac{G_{\mathrm{0}}+G_{\mathrm{1}}+{\eta }_{\mathrm{1}}s}{G_{\mathrm{0}}{G}_{\mathrm{1}}+G_{\mathrm{0}}{\eta }_{\mathrm{1}}s}{\int }_{\mathrm{0}}^{\infty }p(t){\mathrm{e}}^{-st}\mathrm{d}t\end{array}（28）

此即为围岩特性曲线与虚拟支护力在Laplace象空间的对应关系。该围岩特性曲线与时间有关，而前文所得位移释放系数为空间函数，无法耦合求解，因此，需将位移释放系数转换为时间函数。

在开挖面后方，可直接将开挖时刻定义为0时刻起点，利用x=vt将空间距离转换为以时间为自变量的函数。在开挖面前方，位移释放系数的空间函数起始范围为x\to -\mathrm{\infty }。基于起始点为无穷远处的函数无法将时间与空间建立联系，x=\mathrm{0}对应的时间t无法定义。根据前文结果，当x/R=-\mathrm{5}时，隧道位移释放系数基本为0。为减少误差，本文将开挖面距研究断面10倍隧道半径处定为t=\mathrm{0}，在此空间位置处研究断面未受开挖扰动，此处t无实际物理意义，之后的t正常流动，则可将位移释放系数f(t)表达式转换为：

f(t)=a{\mathrm{e}}^{b\frac{\left(vt-\mathrm{10}R\right)}{R}}+c=a{\mathrm{e}}^{\frac{vb}{R}t-\mathrm{10}b}+c,\mathrm{0}\le t\le \frac{\mathrm{10}R}{v}（29）

对式（29）进行Laplace变换：

\frac{{\tilde{u}}_{(\mathrm{r},s)}}{u_{\mathrm{m}\mathrm{a}\mathrm{x}}}=\frac{c}{s}+\frac{a{\mathrm{e}}^{-\mathrm{10}b}}{-\frac{vb}{R}+s}（30）

式中，u_{\mathrm{m}\mathrm{a}\mathrm{x}}=\frac{{\sigma }^{\mathrm{0}}R}{\mathrm{2}}\left(\frac{\mathrm{1}}{G_{\mathrm{0}}}+\frac{\mathrm{1}}{G_{\mathrm{1}}}\right)=\frac{{\sigma }^{\mathrm{0}}R}{\mathrm{2}}·\frac{G_{\mathrm{0}}+G_{\mathrm{1}}}{G_{\mathrm{0}}{G}_{\mathrm{1}}}。

联立式（28）和（30）可得：

\begin{array}{l}\frac{{\sigma }^{\mathrm{0}}R}{\mathrm{2}s}·\frac{\mathrm{1}}{G\left(s\right)}-\frac{R}{\mathrm{2}}·\frac{\mathrm{1}}{G\left(s\right)}{\int }_{\mathrm{0}}^{\infty }p(t){\mathrm{e}}^{-st}\mathrm{d}t=\\ \left(\frac{c}{s}+\frac{a{\mathrm{e}}^{-\mathrm{10}b}}{-\frac{vb}{R}+s}\right)u_{\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}（31）

\begin{array}{l}{\int }_{\mathrm{0}}^{\infty }p(t){\mathrm{e}}^{-st}\mathrm{d}t=\frac{{\sigma }^{\mathrm{0}}}{s}-\frac{\mathrm{2}c}{s}·\frac{G\left(s\right)}{R}{u}_{\mathrm{m}\mathrm{a}\mathrm{x}}-\\ \frac{\mathrm{2}a{\mathrm{e}}^{-\mathrm{10}b}}{\left(-\frac{vb}{R}+s\right)}·\frac{G\left(s\right)}{R}{u}_{\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}（32）

对式（32）进行Laplace逆变换，得到黏弹性围岩中虚拟支护力表达式为：

\begin{array}{l}p(t)={\sigma }^{\mathrm{0}}-\mathrm{2}c\frac{u_{\mathrm{m}\mathrm{a}\mathrm{x}}}{R}\left(\frac{G_{\mathrm{0}}{G}_{\mathrm{1}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}+\frac{G_{\mathrm{0}}^{\mathrm{2}}}{G_{\mathrm{0}}+G_{\mathrm{1}}}{\mathrm{e}}^{-\frac{G_{\mathrm{0}}+G_{\mathrm{1}}}{{\eta }_{\mathrm{1}}}t}\right)-\\ \frac{\mathrm{2}a{\mathrm{e}}^{-\mathrm{10}b}}{R}{u}_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(\frac{G_{\mathrm{0}}^{\mathrm{2}}}{G_{\mathrm{0}}+G_{\mathrm{1}}+{\eta }_{\mathrm{1}}\frac{vb}{R}}{\mathrm{e}}^{-\frac{G_{\mathrm{0}}+G_{\mathrm{1}}}{{\eta }_{\mathrm{1}}}t}+\right.\\ \frac{G_{\mathrm{0}}{G}_{\mathrm{1}}+G_{\mathrm{0}}{\eta }_{\mathrm{1}}\frac{vb}{R}}{G_{\mathrm{0}}+G_{\mathrm{1}}+{\eta }_{\mathrm{1}}\frac{vb}{R}}{\mathrm{e}}^{\frac{vb}{R}t}\left.\begin{array}{c}\\ \end{array}\right)\end{array}（33）

4 黏弹性围岩虚拟支护力影响因素分析

4.1 Kelvin参数对虚拟支护力影响分析

基于式（33）进行分析，将工况N对应的归一化虚拟支护力与位移释放系数曲线绘制成图20。图20中，归一化虚拟支护力随着x/R的增加而减小，位移释放系数随着x/R的增加而变大。

由图20可见，虚拟支护力分布较为集中，待定参数具有较小的波动性，因此，认为Kelvin参数对虚拟支护力无显著影响。

此外，从图20可以看出，对于同一组围岩，位移释放率分布较为离散，而虚拟支护力曲线分布较为集中。因此，对于不同的黏弹性围岩，虚拟支护力表达式的同一组拟合系数预测相较于围岩纵向变形曲线，其波动性更小、鲁棒性更高。

4.2 开挖速度对虚拟支护力影响分析

为了探究不同开挖速度对虚拟支护力的影响，以表4中工况1为例，开挖速度对归一化虚拟支护力的影响如图21所示。图21中，编码越大则开挖速度越慢。由图21可知，开挖速度越快，在距开挖面相同位置处虚拟支护力更大。当开挖速度低于一定值后，对应虚拟支护力趋于稳定，此时速度对虚拟支护力的影响不大。

将表4、5中130组工况的虚拟支护力变化汇总成图22。由图22可见，开挖面前方虚拟支护力更为集中，在隧道开挖面处的归一化虚拟支护力处于0.664～0.841区间。开挖面后方虚拟支护力与开挖速度具有强相关性，曲线由上包络线至下包络线之间由疏至密。不同开挖速度对应的虚拟支护力待定参数见图23。图23中，a-、b-、c-为开挖面前方（x\le0）的虚拟支护力待定参数，a+、b+为开挖面后方（x>0）的虚拟支护力待定参数；横坐标的系列1～9对应工况1‒1～1‒9，系列10～18对应工况2‒1～2‒9，以此类推，系列118～126对应工况14‒1～14‒9。由图23可见，随着开挖速度减小，各拟合系数整体趋势呈现周期性波动，但波动范围较小。

5 虚拟支护力简化拟合公式

由于式（33）涉及参数过多，难以满足实际工程使用的便捷性要求，为此，本节结合前文影响因素研究规律，提出虚拟支护力的简化拟合公式。虚拟支护力有着更高的鲁棒性，与开挖速度有较为明显的相关关系，可由上下包络线及平均值曲线进行预测，隧道开挖速度越快，虚拟支护力越靠近上包络线。为此，提出虚拟支护力简化拟合公式如下：

p(t)=A{\mathrm{e}}^{\frac{t/R}{B}}+C（34）

式中，A、B、C为待定参数。

由于时空效应的演化规律，开挖面前方（x\le0）和后方（x>0）所对应的拟合参数不同，本文算例中虚拟支护力变化及包络线如图24所示。

由于开挖速度的影响，上述待定参数在开挖速度变化时可通过线性差值得到。在传统的弹塑性研究中，学者广泛认可的虚拟支护力公式主要包括两种形式^[34‒35]。

第1类负指数形式为：

p_{\mathrm{x}}/{\sigma }^{\mathrm{0}}=\alpha {\mathrm{e}}^{-x/R_{\mathrm{L}}},x\ge \mathrm{0}（35）

式中：\alpha与围岩参数相关，一般取0.55～0.75；R_{\mathrm{L}}为开挖影响半径，取决于围岩性质和施工方法，可取为0.7～1.0倍隧道开挖半径。

第2类反比例函数形式为：

p_{\mathrm{x}}/{\sigma }^{\mathrm{0}}=\alpha {\left(\frac{\mathrm{1}}{x/R_{\mathrm{L}}+\mathrm{1}}\right)}^{\mathrm{2}},x\ge \mathrm{0}（36）

将本文求解结果与既有经验公式^[34‒36]对比，\alpha取值分别为0.55、0.6、0.65、0.70、0.75，R_{\mathrm{L}}同样取5个梯度，分别为0.700、0.775、0.850、0.925、1.000。简化拟合公式与既有虚拟支护力表达式对比，结果如图25所示。

由图25可见，既有研究^[34‒36]的虚拟支护力计算结果位于本文计算所得包络线范围内且靠近下限，这是由于在开挖速度较慢时，黏弹性围岩的虚拟支护力趋近于弹性状态下，而既有研究^[34‒36]的解答无法适用较快的开挖速度工况。当开挖速度较大时，现有方法计算的虚拟支护力偏小，将导致支护设计偏强。较既有研究的经验公式，本文提出的虚拟支护力简化拟合公式适用范围更广，能更好地模拟不同开挖速度下黏弹性围岩的3维时空效应。

6 结 论

本文从黏弹性围岩纵向变形的时空演化特征出发，通过大量数值模拟和回归分析获得了纵向变形曲线的经验公式，进一步对黏弹性围岩的虚拟支护力进行了解析，分析其主要影响因素，并提出了虚拟支护力的简化拟合方法。主要结论如下：

1）揭示了黏弹性围岩纵向变形的时空演化规律，分析表明地应力对位移释放系数无影响。当开挖速度较快时，位移释放系数随开挖速度增加而减小；当开挖速度足够慢时，黏弹性围岩位移释放系数只存在空间效应影响。

2）黏弹性围岩位移释放系数在隧道纵向上可分为3部分：第1部分为，此阶段由于无开挖效应而分布较为集中；第2部分为，此阶段时空效应共同作用，受诸多因素影响而较为离散；第3部分为x\ge \mathrm{10}R，此阶段空间效应影响可忽略，时间效应影响显著。

3）开挖面前方位移释放系数曲线簇较为集中，而后方随时空效应变化较大。开挖面后方的时空效应可分为空间效应控制型、时间效应控制型和时空效应共同影响型等3种，并给出了各类型的判据及对应经验公式求解方法。

4）基于黏弹性围岩位移释放系数获得了虚拟支护力解答，分析表明，当位移释放系数为空间效应控制型时，虚拟支护力与位移释放系数满足翻转对称。总体而言，虚拟支护力曲线簇鲁棒性更高，呈现集中趋势，这也使得其简化拟合成为可能。

5）虚拟支护力随开挖速度减慢而降低，最终趋近于只存在空间效应解，而当速度加快时则应力释放减弱。为此提出了其简化拟合公式以便于工程初步设计时使用，较既有经验公式适用范围更广，能更好地反映黏弹性围岩的3维时空效应。

参考文献

[1] Zhang Dingli.Essential issues and their research progress in tunnel and underground engineering[J].Chinese Journal of Theoretical and Applied Mechanics,2017,49(1):3‒21.

[2] 张顶立.隧道及地下工程的基本问题及其研究进展[J].力学学报,2017,49(1):3‒21.

[3] Sun Zhenyu, Zhang Dingli, Fang Qian,et al.Analysis of the interaction between tunnel support and surrounding rock considering pre-reinforcement[J].Tunnelling and Underground Space Technology,2021,115:104074. doi:10.1016/j.tust.2021.104074

[4] Sun Zhenyu, Zhang Dingli, Fang Qian,et al.Displacement process analysis of deep tunnels with grouted rockbolts considering bolt installation time and bolt length[J].Computers and Geotechnics,2021,140:104437. doi:10.1016/j.compgeo.2021.104437

[5] Sun Zhenyu, Zhang Dingli, Fang Qian,et al.Convergence-confinement analysis for tunnels with combined bolt-cable system considering the effects of intermediate principal stress[J].Acta Geotechnica,2023,18(6):3323‒3348. doi:10.1007/s11440-022-01758-4

[6] Cui Lan, Zhen Junjie, Miao Chenxi,et al.Coupling analysis of longitudinal deformation profile and ground reaction curve[J].Chinese Journal of Geotechnical Engineering,2014,36(4):707‒715.

[7] 崔岚,郑俊杰,苗晨曦,等.隧道纵向变形曲线与围岩特征曲线耦合分析[J].岩土工程学报,2014,36(4):707‒715.

[8] Sun Zhenyu, Zhang Dingli, Liu Daoping,et al.Analysis of the whole-process mechanical response of tunnel surrou-nding rock under the effect of anchorage system[J].Engin-eering Mechanics,2022,39(7):170‒182.

[9] 孙振宇,张顶立,刘道平,等.锚固体系作用下隧道围岩力学响应的全过程解析[J].工程力学,2022,39(7):170‒182.

[10] Zhang Changguang, Zeng Kaihua.Comparisons of spatial-effect approaches for tunnel excavation using convergence-confinement method[J].Rock and Soil Mechanics,2016,37(5):1417‒1424.

[11] 张常光,曾开华.收敛约束法中隧道开挖面空间效应方法比较[J].岩土力学,2016,37(5):1417‒1424.

[12] Wu Shunchuan, Geng Xiaojie, Gao Rongtao,et al.A study of the longitudinal deformation of tunnels based on the ge-neralized Hoek‒Brown failure criterion[J].Rock and Soil Mechanics,2015,36(4):946‒953.

[13] 吴顺川,耿晓杰,高永涛,等.基于广义Hoek‒Brown准则的隧道纵向变形曲线研究[J].岩土力学,2015,36(4):946‒953.

[14] Sun Zhenyu, Zhang Dingli, Fang Qian,et al.Spatial and temporal evolution characteristics of interaction between primary support and tunnel surrounding rock[J].Chinese Journal of Rock Mechanics and Engineering,2017,36(S2):3943‒3956.

[15] 孙振宇,张顶立,房倩,等.隧道初期支护与围岩相互作用的时空演化特性[J].岩石力学与工程学报,2017,36(S2):3943‒3956.

[16] Panet M, Guenot A.Analysis of convergence behind the face of a tunnel[C]//In:Proceedings of the International Symposium Tunnelling,London:IMM,1982:197‒204. doi:10.1016/0148-9062(83)91744-8

[17] Liang Yiwen, Zha Wenhua, Liu Zaobao,et al.Research on the prediction of the time of initial support of soft rock tunnel based on the optimized longitudinal deformation profile[J].Applied Sciences,2022,12(14):7221. doi:10.3390/app12147221

[18] Chern J C, Shiao F Y, Yu C W.An empirical safety criterion for tunnel construction[C]//Proceedings of the regional sy-mposium on sedimentary rock engineering,Taiwan:Balke-ma,1998:222‒227.

[19] Vlachopoulos N, Diederichs M S.Improved longitudinal displacement profiles for convergence confinement analysis of deep tunnels[J].Rock Mechanics and Rock Engineering,2009,42(2):131‒146. doi:10.1007/s00603-009-0176-4

[20] Alejano L R, Rodríguez‒Dono A, Veiga M.Plastic radii and longitudinal deformation profiles of tunnels excavated in strain-softening rock masses[J].Tunnelling and Underground Space Technology,2012,30:169‒182. doi:10.1016/j.tust.2012.02.017

[21] Wang Zhilong, Wang Mingnian, Liu Dagang.Analysis and modeling of longitudinal deformation profiles considering the reinforcement effect of tunnel face in weak rocks[J].KSCE Journal of Civil Engineering,2022,26(12):5372‒5383. doi:10.1007/s12205-022-1664-2

[22] Luo Yanbin, Chen Jianxun, Chen Yi,et al.Longitudinal deformation profile of a tunnel in weak rock mass by using the back analysis method[J].Tunnelling and Underground Space Technology,2018,71:478‒493. doi:10.1016/j.tust.2017.10.003

[23] Chu Zhaofei, Liu Baoguo, Liu Kaiyun,et al.Analytical viscoelastic solutions for lined circular tunnels under two co-ntact conditions in a non-hydrostatic stress[J].Rock and Soil Mechanics,2017,38(11):3215‒3224.

[24] 储昭飞,刘保国,刘开云,等.非静水应力场中圆形隧道衬砌与围岩间两种接触的黏弹性解析[J].岩土力学,2017,38(11):3215‒3224.

[25] Zhao Nannan, Shao Zhushan, Wu Kui,et al.Time-dependent solutions for lined circular tunnels considering rockbolts reinforcement and face advancement effects[J].Internatio-nal Journal of Geomechanics,2021,21(10):04021179. doi:10.1061/(asce)gm.1943-5622.0002130

[26] Kabwe E, Karakus M, Chanda E K.Time-dependent solution for non-circular tunnels considering the elasto-viscoplastic rockmass[J].International Journal of Rock Mechanics and Mining Sciences,2020,133:104395. doi:10.1016/j.ijrmms.2020.104395

[27] Do D P, Tran N T, Mai V T,et al.Time-dependent reliability analysis of deep tunnel in the viscoelastic burger rock with sequential installation of liners[J].Rock Mechanics and Rock Engineering,2020,53(3):1259‒1285. doi:10.1007/s00603-019-01975-6

[28] Kargar A R, Haghgouei H, Babanouri N.Time-dependent analysis of stress components around lined tunnels with circular configuration considering tunnel advancing rate effects[J].International Journal of Rock Mechanics and Mining Sciences,2020,133:104422. doi:10.1016/j.ijrmms.2020.104422

[29] Zeng G S, Wang H N, Jiang M J,et al.Analytical solution of displacement and stress induced by the sequential excavation of noncircular tunnels in viscoelastic rock[J].International Journal of Rock Mechanics and Mining Sciences,2020,134:104429. doi:10.1016/j.ijrmms.2020.104429

[30] Zhou Zihan, He Chuan, Chen Ziquan,et al.Virtual supporting force method for tunnel support timing based on the concept of active deformation control of surrounding rock[J].China Railway Science,2022,43(6):76‒88.

[31] 周子寒,何川,陈子全,等.基于围岩变形主动控制理念的隧道支护时机虚拟支撑力法[J].中国铁道科学,2022,43(6):76‒88.

[32] Zhang Dingli, Sun Zhenyu, Hou Yanjuan.Tunnel supporting structure system and its synergy[J].Chinese Journal of Theoretical and Applied Mechanics,2019,51(2):577‒593.

[33] 张顶立,孙振宇,侯艳娟.隧道支护结构体系及其协同作用.力学学报,2019,51(2):577‒593.

[34] Sun Zhenyu.The synergetic principle and design method of tunnel support system[D].Beijing:Beijing Jiaotong University,2020.

[35] 孙振宇.隧道支护体系协同作用原理与设计方法[D].北京:北京交通大学,2020.

[36] Cui Lan, Zheng Junjie, Zhang Rongjun,et al.A numerical procedure for the fictitious support pressure in the application of the convergence-confinement method for circular tunnel design[J].International Journal of Rock Mechanics and Mining Sciences,2015,78:336‒349. doi:10.1016/j.ijrmms.2015.07.001

[37] Huang Xiaohua, Jin Yanli, Huang Shuling,et al.Conversion between two types of parameters in the viscoelastic combination model for common hypotheses and its applicability[J].Chinese Journal of Solid Mechanics,2019,41(1):59‒70.

[38] 黄小华,金艳丽,黄书岭,等.常用假设下黏弹组合模型两类表述参数间转换及适用性[J].固体力学学报,2020,41(1):59‒70.

[39] Chu Zhaofei, Wu Zhijun, Liu Quansheng,et al.Viscos-elastic-plastic solution for deep buried tunnels considering tunnel face effect and sequential installation of double linings[J].Computers and Geotechnics,2024,165:105930. doi:10.1016/j.compgeo.2023.105930

[40] Zhang Zhiliang, Xu Weiya, Wang Wei.Study of triaxial creep tests and its nonlinear visco-elastoplastic creep model of rock from compressive zone of dam foundation in Xiajiaba Hydropower Station[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(1):132‒140.

[41] 张治亮,徐卫亚,王伟.向家坝水电站坝基挤压带岩石三轴蠕变试验及非线性黏弹塑性蠕变模型研究[J].岩石力学与工程学报,2011,30(1):132‒140.

[42] Sun Zhenyu, Zhang Dingli, Hou Yanjuan,et al.Whole-process deformation laws and determination of stability criterion of surrounding rock of tunnels based on statistics of field measured data[J].Chinese Journal of Geotechnical Engineering,2021,43(7):1261‒1270.

[43] 孙振宇,张顶立,侯艳娟,等.基于现场实测数据统计的隧道围岩全过程变形规律及稳定性判据确定[J].岩土工程学报,2021,43(7):1261‒1270.

[44] Corbetta F, Bernaud D, Nguyen Minh D.Contribution à la méthode convergence-confinement par le principe de la similitude[J].Revue Française de Géotechnique,1991(54):5‒11. doi:10.1051/geotech/1991054005

[45] Sun Zhenyu, Zhang Dingli, Fang Qian,et al.Analytical solutions for deep tunnels in strain-softening rocks modeled by different elastic strain definitions with the unified strength theory[J].Science China Technological Sciences,2022,65(10):2503‒2519. doi:10.1007/s11431-022-2158-9

[46] Lee Yulin, Xu Wengui, Lin Mingyi.Analysis of advancing effect and confinement loss by using deformation measurements in tunnel[J].Chinese Journal of Rock Mechanics and Engineering,2009,28(1):39‒46.

[47] 李煜舲,许文贵,林铭益.以隧道变形量测资料分析掘进效应与约束损失[J].岩石力学与工程学报,2009,28(1):39‒46.

[48] Paraskevopoulou C, Diederichs M.Analysis of time-de-pendent deformation in tunnels using the Convergence‒Confinement Method[J].Tunnelling and Underground Space Technology,2018,71:62‒80. doi:10.1016/j.tust.2017.07.001

[49] Liu Chang, Zhang Dingli, Zhang Sulei,et al.Long-term mechanical analysis of tunnel structures in rheological rock considering the degradation of primary lining[J].Undergr-ound Space,2023,10:217‒232. doi:10.1016/j.undsp.2022.10.005

基金资助

国家自然科学基金项目(52578456)

国家自然科学基金项目(52278387)

国家自然科学基金项目(52178379)

国家自然科学基金项目(52208382)