通过分析围岩特征曲线（GRC）与隧道纵剖面变形曲线（LDP）确定支护结构的架设参数，是收敛约束法的重要组成部分^[1–2]。其中，GRC描述围岩横断面变形行为，LDP反映隧道开挖面空间效应。学者们普遍认为GRC主要是指黏性地层特征线，且隧道围岩稳定性与成洞时间具有复杂关联^[1]。为客观地确定岩质隧道施工过程支护参数及评价隧道稳定性，成洞时间效应与开挖面空间效应的影响（时空效应）必须同时考虑，缺一不可。

关于围岩时空变形问题的研究相继得到开展^[3–7]。卞跃威等^[3]基于Kelvin体和Mohr–Coulomb组合模型与非关联流动法则，推导了应力释放影响下隧道时效变形解。蔡燕燕等^[4]基于非线性西原体和Hoek–Brown组合模型与非关联流动法则，推导了开挖面约束作用下围岩蠕变位移计算式。Manh等^[5]基于Burgers体和Mohr–Coulomb组合模型（CVISC模型）与非关联流动法则，计算了大挤压地层中隧道开挖的黏塑性变形解。温森等^[6]利用西原体和广义Hoek–Brown组合模型与非关联流动法则推导了隧道机械开挖的流变变形解。Birchall等^[7]以广义变分原理为基础采用3维能量法研究了应力释放对围岩时效变形的影响。以上研究在描述围岩流变力学行为时，采用了经典的蠕变模型（如西原体与Burgers体）与屈服准则（如Hoek–Brown准则和Mohr–Coulomb准则）。这些经典模型成功描述了围岩流变力学行为，但上述研究在描述开挖面空间效应时均采用应力释放系数法。由于应力释放率与开挖距离、围岩性态等的函数关系难以确定^[8–9]，应力释放系数法的应用存在诸多限制。为克服这一困难，不少学者假定在隧道开挖面推进过程中，应力释放率近似等效于位移释放率^[10–11]。在弹性变形阶段，这个假定尚可接受。但在塑性阶段，由于围岩变形不可逆且与应力加卸载历史相关，做此假定还需考究其合理有效性。

目前，隧道开挖过程中位移释放问题的研究渐趋成熟^[8,12–14]，位移释放系数法已广泛应用于该问题的数值研究中^[15]。经证实，位移释放系数法能够完备地描述开挖面空间效应。“应力释放率近似等效于位移释放率”这一假定正是利用了位移释放系数法的完备性。作者尝试利用位移释放系数法解析岩质隧道施工过程的时空变形问题，在计算过程中，采用Burgers体和Drucker–Prager组合模型与非关联流动法则，最终建立了岩质隧道施工过程中时空变形封闭解。详细论证了变形封闭解的正确有效性，并分析了主要参数的影响规律。

1 岩质隧道计算模型 1.1 计算模型与基本假设

如图1所示，深埋隧道半径为R₀。假设地应力场为静水压力状态，围岩为各向同性均质体。当地应力达到一定水平时，围岩形成塑性圈，其半径假设为R_p。随着时间推移，围岩蠕变行为逐渐显现。如图1所示，围岩由黏弹性区与塑性区构成。黏弹性区总应变包括弹性应变与黏弹性应变，塑性区总应变包括弹性应变与塑性应变。假设围岩瞬时力学性态服从Drucker–Prager 准则，蠕变力学性态服从Burgers模型，由于隧道开挖过程开挖面具有空间效应，其前后方一定范围内围岩力学响应均受到约束作用。如图2所示，当时间由t₁推移至t₂时，断面AA与无穷远处断面的蠕变变形量并不一致。可以推断，由于断面AA受到开挖面的约束作用，其蠕变量偏小，位移释放极不充分。作者采用位移释放系数法描述岩质隧道施工过程变形的时空问题。

1.2 位移释放系数法

位移释放系数 ${u^ * }$ 可由Vlachopoulos和Diederichs 以围岩塑性区最大半径R_max为基础的公式（简称V–D（09）方程）^[8,12]得出：

${u^ * } = \frac{{{u_{\rm{r}}}}}{{{u_{{\rm{r}},\max }}}} = \left\{ \begin{aligned}& u_0^ * {{\rm{e}}^{{X^ * }}},\;\;{X^ * } \le 0{\text{；}}\\ & 1 - \left( {1 - u_0^ * } \right){{\rm{e}}^{\textstyle{\frac{{ - 3{X^ * }}}{{2{R^ * }}}}}},\;\;{X^ * } \ge 0 \end{aligned} \right.$

(1)

$u_0^ * = \frac{1}{3}{{\rm{e}}^{ - 0.15{R^ * }}}\mathop ,\nolimits^{} \mathop {}\nolimits^{} x = 0$

(2)

式中： $X^ *$ 为相对距离， $X^ *$ =x/R₀，x为断面AA至开挖面的纵向距离；u_r为断面AA径向位移；u_r,max为距开挖面无穷远处不受开挖面约束作用的2维平面应变无支护断面产生的最大径向位移； $R^ *$ 为相对半径， $R^ *$ =R_max/R₀，R_max为2维平面应变无支护围岩的最大塑性区半径； $u_0^ *$ 为开挖面处位移释放率。

式（1）描述了隧道纵剖面变形曲线。当已知2维平面应变无支护隧道最大径向位移u_r,max时，通过式（1）可预测开挖面约束作用下各断面的变形情况。V–D（09）方程普遍适用于弹性围岩（ $R^ *$ =1）和各种弹塑性围岩（ $R^ * $ >1），且反映了隧道埋深（地应力）、岩体特性（强度准则参数与中间主应力效应）和施工方法（围岩扰动程度）等因素的综合影响，可在复杂工程地质环境中推广应用^[8,12]。图2中，卸荷段围岩蠕变历时t₁与t₂后，洞壁位移值增大。假设不同断面洞壁位移演化满足V–D（09）方程，则2维平面应变无支护隧道最大径向位移u_r,max应当视为时间的函数，即蠕变变形。

1.3 Burgers体和Drucker–Prager组合模型

Burgers体和Drucker–Prager组合模型为Burgers模型与Drucker–Prager准则组合的黏弹塑性模型。如图3所示，组合模型能模拟弹塑性体积行为和黏弹塑性偏量特性。Burgers模型本构关系为^[5]：

$\begin{aligned}[b] & \qquad 2{\eta ^{\rm{K}}}{{\ddot e}_{ij}} + 2{G^{\rm{K}}}{{\dot e}_{ij}} = \\ & \displaystyle\frac{{{\eta ^{\rm{K}}}}}{{{G^{\rm{M}}}}}{{\ddot S}_{\!ij}} + \left( {1 + \displaystyle\frac{{{G^{\rm{K}}}}}{{{G^{\rm{M}}}}} + \frac{{{\eta ^{\rm{K}}}}}{{{\eta ^{\rm{M}}}}}} \right){{\dot S}_{\!ij}} + \frac{{{G^{\rm{K}}}}}{{{\eta ^{\rm{M}}}}}{S_{\!ij}} \end{aligned}$

(3)

式中： $G^{\rm K}$ 为Kelvin剪切模量； $ \eta^{\rm K}$ 为Kelvin黏滞系数； $G^{\rm M}$ 为Maxwell剪切模量； $ \eta^{\rm M}$ 为Maxwell黏滞系数；S_ij和e_ij分别为偏应力张量与偏应变张量，有：

${S_{\!ij}} = {\sigma _{ij}} - \frac{1}{3}{\delta _{ij}}{\sigma _{kk}}, \; {e_{ij}} = {\varepsilon _{ij}} - \frac{1}{3}{\delta _{ij}}{\varepsilon _{kk}}$

(4)

式中： $\delta_{ij}$ 为Kronecker符号； $ \sigma _{ij}$ 和 $\varepsilon_{ij}$ 分别为应力张量与应变张量； $ \sigma _{kk}$ 和 $\varepsilon _{kk}$ 分别为体积应力与体积应变，存在

${\sigma _{kk}} = 3K{\varepsilon _{kk}}$

(5)

式中，K为体积模量， $K = \displaystyle\frac{{2{G^{\rm{M}}}\left( {1 + \upsilon } \right)}}{{3\left( {1 - 2\upsilon } \right)}}$ ， $ \upsilon$ 为泊松比。

Drucker–Prager 准则为^[16]：

$F = \sqrt {{J_2}} - \alpha {I_1} -\kappa= 0$

(6)

式中：I₁ 为应力张量第一不变量；J₂为偏应力张量第二不变量； $\alpha$ 和 $\kappa$ 均为与岩石的内摩擦角和黏聚力有关的试验常数，计算式为：

$\alpha = \frac{{2\sin \;\varphi }}{{\sqrt 3 \left( {3 - \sin\; \varphi } \right)}}$

(7)

$\kappa= \frac{{6c\cos \;\varphi }}{{\sqrt 3 \left( {3 - \sin\; \varphi } \right)}}$

(8)

式中，c为黏聚力， $ \varphi$ 为内摩擦角。

由非关联流动法则可确定塑性势函数为：

$Q = \sqrt {{J_2}} - {\alpha _{ψ} }{I_1} - \kappa = 0$

(9)

式中， ${\alpha _{ψ} }$ 为扩容参数，且

${\alpha _{ψ} } = \frac{{\sin \;{ψ} }}{{\sqrt 3 \sqrt {3 + {{\sin }^2}{ψ} } }}$

(10)

式中， ${ψ}$ 为扩容角。

塑性应变可通过塑性势函数求解^[17]：

$\dot \varepsilon _{ij}^p = \dot \lambda \frac{{\partial Q}}{{\partial {\sigma _{ij}}}} = \dot \lambda \left( { - {\alpha _{ψ} }{\delta _{ij}} + \frac{{{S_{\!ij}}}}{{2\sqrt {{J_2}} }}} \right)$

(11)

式中， $\dot \varepsilon _{ij}^p$ 为塑性应变率， $\dot \lambda $ 为塑性乘子。

2 围岩变形时空效应解析

岩质隧道开挖过程围岩变形的时空效应分析的关键在于2维平面应变无支护隧道最大径向位移u_r,max的求解。以下为计算过程：

2.1 黏弹性区计算

基于Drucker–Prager准则的围岩开始屈服时的原岩应力p₀为^[18]：

${p_0} = \frac{k}{{1 - 3\alpha }}$

(12)

当原岩应力大于应力阈值式（12）时，围岩力学行为表现为塑性。式（13）给出了塑性圈存在下弹性区位移的计算公式^[19]：

${u_{\rm{r}}} = \frac{{R_{\rm{p}}^2}}{{2Gr}}\left( {{p_0}\sin \;\varphi + c\cos\; \varphi } \right)$

(13)

式中，R_p为围岩塑性区半径，r为围岩径向深度。

根据弹性–黏弹性对应原理^[20]，有：

${u_{\rm{r}}}\left( t \right) = \frac{{R_{\rm{p}}^2\left( {1 + \upsilon } \right)J\left( t \right)}}{r}\left( {{p_0}\sin\; \varphi + c\cos\; \varphi } \right)$

(14)

式中，J(t)为Burgers体的蠕变柔量，可对式（3）进行拉普拉斯变换及其逆变换得到，计算公式为：

$J\left( t \right) = \frac{1}{{2{G^{\rm{M}}}}} + \frac{t}{{2{\eta ^{\rm{M}}}}} + \frac{1}{{2{G^{\rm{K}}}}}\left( {1 - {{\rm{e}}^{ - \frac{{{G^{\rm{K}}}}}{{{\eta ^{\rm{K}}}}}t}}} \right)$

(15)

将式（15）代入式（14）可得岩质隧道黏弹性变形解。进而可得围岩弹塑性交界面上的位移，即

$\begin{aligned}[b] {u_{\rm{r}}}\left( {{R_{\rm{p}}},t} \right) = {R_{\rm{p}}}\left( {1 + \upsilon } \right)\left( {{p_0}\sin\; \varphi + c\cos\; \varphi } \right) \cdot \\ \left[ {\displaystyle\frac{t}{{2{\eta ^{\rm{M}}}}} +\displaystyle\frac{1}{{2{G^{\rm{M}}}}} + \displaystyle\frac{1}{{2{G^{\rm{K}}}}} - \displaystyle\frac{1}{{2{G^{\rm{K}}}}}\exp \left( { - \frac{{{G^{\rm{K}}}}}{{{\eta ^{\rm{K}}}}}t} \right)} \right] \end{aligned}$

(16)

围岩弹塑性交界面变形协调，式（16）为计算塑性区变形的边界条件。

2.2 塑性区计算

针对平面应变问题，轴向塑性应变率 $\dot \varepsilon _{\textit{z}}^{\rm{p}} = 0$ 。将该等式与式（11）联立，可得围岩塑性区轴向偏应力 $S_{\textit{z}}^{\rm{p}}$ 和轴向应力 $\sigma _{\textit{z}}^{\rm{p}}$ ：

$S_{\textit{z}}^{\rm{p}} = 2{\alpha _{ψ} }\sqrt {{J_2}} $

(17)

$\sigma _{\textit{z}}^{\rm{p}} = S_{\textit{z}}^{\rm{p}} + \frac{{{I_1}}}{3} = 2{\alpha _{ψ} }\sqrt {{J_2}} + \frac{{{I_1}}}{3}$

(18)

由此可得：

$\sigma _{\textit{z}}^{\rm{p}} = \frac{1}{2}\left( {\sigma _\theta ^{\rm{p}} + \sigma _{\rm{r}}^{\rm{p}}} \right) + 3{\alpha _{ψ} }\sqrt {{J_2}} $

(19)

式中， $\sigma _\theta ^{\rm{p}}$ 和 $\sigma _{\rm{r}}^{\rm{p}}$ 分别为围岩塑性区环向应力和径向应力。

因而偏应力张量第二不变量J₂可表示为：

$\begin{aligned}[b] &{J_2} = \displaystyle\frac{1}{6}\left[ {{{\left( {\sigma _\theta ^{\rm{p}} - \sigma _{\textit{z}}^{\rm{p}}} \right)}^2} + {{\left( {\sigma _{\textit{z}}^{\rm{p}} - \sigma _{\rm{r}}^{\rm{p}}} \right)}^2} + {{\left( {\sigma _{\rm{r}}^{\rm{p}} - \sigma _\theta ^{\rm{p}}} \right)}^2}} \right]= \\ &\qquad \displaystyle\frac{1}{6}\left[ {\frac{3}{2}{{\left( {\sigma _\theta ^{\rm{p}} - \sigma _{\rm{r}}^{\rm{p}}} \right)}^2} + 18\alpha _{ψ} ^2{J_2}} \right] \end{aligned}$

(20)

式（20）可化简为：

$\sqrt {{J_2}} = \sqrt {\frac{{3 + {{\sin }^2}{ψ} }}{{12}}} \left( {\sigma _\theta ^{\rm{p}}- \sigma _{\rm{r}}^{\rm{p}}} \right)$

(21)

联立式（18）与（19），可得：

${I_1} = 3{\alpha _{ψ} }\sqrt {{J_2}} + \frac{3}{2}\left( {\sigma _{\rm{r}}^{\rm{p}} + \sigma _\theta^{\rm{p}}} \right)$

(22)

将式（21）代入式（22），可得：

${I_1} = \frac{{\left( {3 + \sin\; {ψ} } \right)\sigma _\theta ^{\rm{p}} + \left( {3 - \sin\; {ψ} } \right)\sigma _{\rm{r}}^{\rm{p}}}}{2}$

(23)

联立Drucker–Prager准则式（5）与式（22）～（23），可得：

$\sigma _\theta ^{\rm{p}} - \sigma _{\rm{r}}^{\rm{p}} = \frac{{2\left( {\sigma _{\rm{r}}^{\rm{p}} + c\cot\; \varphi } \right)}}{{\csc\; \varphi - 1}}$

(24)

该力学问题的平衡方程为^[21–22]：

$\sigma _\theta ^{\rm{p}} - \sigma _{\rm{r}}^{\rm{p}} = r \displaystyle\frac{{{\rm{d}}\sigma _{\rm{r}}^{\rm{p}}}}{{{\rm{d}}r}}$

(25)

针对平面应变条件，当 $r = {R_0}$ 时，有：

$\sigma _{\rm{r}}^{\rm{p}} = 0$

(26)

联立式（22）～（26），可得围岩应力场：

$\left\{ \begin{array}{l} {\sigma^{\rm{p}} _{\rm{r}}} = c\cot\; \varphi \left[ {{{\left( {\displaystyle\frac{r}{{{R_0}}}} \right)}^{\frac{2}{{\csc\; \varphi - 1}}}} - 1} \right]{\text{，}}\\ {\sigma^{\rm{p}} _\theta } = c\cot\; \varphi \left[ {\displaystyle\frac{{\csc\; \varphi + 1}}{{\csc\; \varphi - 1}}{{\left( {\displaystyle\frac{r}{{{R_0}}}} \right)}^{\frac{2}{{\csc \;\varphi - 1}}}} - 1} \right] \end{array} \right.$

(27)

围岩弹塑性交界面的应力边界条件为：

$ \sigma _{\rm{r}}^{\rm{p}} + \sigma _\theta ^{\rm{p}} = 2{p_0} $

(28)

将式（28）代入式（27），可得R_p的计算公式：

${R_{\rm{p}}} = {R_0}{\left[ {\frac{{\left( {{p_0} + c\cot\; \varphi } \right)\left( {1 - \sin\; \varphi } \right)}}{{c\cot\; \varphi }}} \right]^{\frac{{1 - \sin\; \varphi }}{{2\sin\; \varphi }}}}$

(29)

由边界条件（26）可知，式（29）塑性区半径即为2维平面应变无支护围岩的最大塑性区半径R_max。Park和Kim^[21]给出了2维平面应变无支护围岩径向位移u_r,max的微分方程：

$\displaystyle\frac{{{\rm{d}}{u_{{\rm{r}},\max }}}}{{{\rm{d}}r}} + {\kappa _{ψ} }\displaystyle\frac{{{u_{{\rm{r}},\max }}}}{r} = f\left( r \right)$

(30)

$f\left( r \right) = {\varepsilon _{\rm{r}}} + {\kappa _{ψ} }{\varepsilon _\theta } = \left( {\varepsilon _{\rm{r}}^{\rm{e}} + \varepsilon _{\rm{r}}^{\rm{p}}} \right) + {\kappa _{ψ} }\left( {\varepsilon _\theta ^{\rm{e}} + \varepsilon _\theta ^{\rm{p}}} \right)$

(31)

其中，塑性应变分量满足^[21]：

$\varepsilon _{\rm{r}}^{\rm{p}} + {\kappa _{ψ} }\varepsilon _\theta ^{\rm{p}} = 0$

(32)

式中， $ {\kappa _{ψ} }$ 为扩容系数，计算公式为：

${\kappa _{ψ} } = \frac{{1 + \sin\; {ψ} }}{{1 - \sin \;{ψ} }}$

(33)

式（30）的位移通解为：

${u_{{\rm{r}},\max }} = \frac{1}{{{r^{{\kappa _{ψ} }}}}}\int_{{R_{\rm p}}}^r {{r^{{\kappa _{ψ} }}}f\left( r \right)} {\rm{d}}r + {u_{\rm{r}}}\left( {{R_{\rm{p}}},t} \right) \cdot {\left( {\frac{{{R_{\rm{p}}}}}{r}} \right)^{{\kappa _{ψ} }}}$

(34)

式中，u_r(R_p，t)为围岩弹塑性交界面上的位移。

围岩塑性区弹性应变 $\varepsilon _{\rm{r}}^{\rm{e}}$ 、 $\varepsilon _\theta ^{\rm{e}}$ 可表示为^[21]：

$\left\{ \begin{aligned} &\varepsilon _{\rm{r}}^{\rm{e}} = \frac{1}{{2{G^{\rm{M}}}}}\left[ {\left( {1 - \upsilon } \right)\left( {\sigma _{\rm{r}}^{\rm{p}} - {p_0}} \right) - \upsilon \left( {\sigma _\theta ^{\rm{p}} - {p_0}} \right)} \right]{\text{，}}\\ &\varepsilon _\theta ^{\rm{e}} = \frac{1}{{2{G^{\rm{M}}}}}\left[ {\left( {1 - \upsilon } \right)\left( {\sigma _\theta ^{\rm{p}} - {p_0}} \right) - \upsilon \left( {\sigma _{\rm{r}}^{\rm{p}} - {p_0}} \right)} \right] \end{aligned} \right.$

(35)

将式（35）代入式（31），可得：

$\begin{aligned}[b] &{f\left( r \right) = \dfrac{1}{{2{G^{\rm{M}}}}}\left[ {\left( {1 - \upsilon - {\kappa _{ψ} }\upsilon } \right)\sigma _{\rm{r}}^{\rm{p}} + \left( {{\kappa _{ψ} } - {\kappa _{ψ} }\upsilon - \upsilon } \right)} \right.}\cdot\\ & {\quad \quad \;\;\;\left. {\sigma _\theta ^{\rm{p}} - \left( {1 - 2\upsilon } \right)\left( {1 + {\kappa _{ψ} }} \right){p_0}} \right]} \end{aligned}$

(36)

将式（36）代入式（34），可得：

$\begin{aligned}[b] {u_{r,\max }}= & {C_1}\left[ {{{\left( {\displaystyle\frac{r}{{{R_0}}}} \right)}^{\frac{2}{{\csc\; \varphi - 1}}}}r - {{\left( {\displaystyle\frac{{{R_{\rm{p}}}}}{{{R_0}}}} \right)}^{\frac{2}{{\csc\; \varphi - 1}}}}\displaystyle\frac{\mathop R\nolimits_p^{{\kappa _{ψ} } + 1}}{{{r^{{\kappa _{ψ} }}}}}} \right]- \\ &{C_2}\left[ {r - \displaystyle\frac{\mathop R\nolimits_{\rm{p}}^{{\kappa _{ψ} } + 1}}{{{r^{{\kappa _{ψ} }}}}}} \right] + {u_{\rm{r}}}\left( {{R_{\rm{p}}},t} \right) \cdot {\left( {\displaystyle\frac{{{R_{\rm{p}}}}}{r}} \right)^{{\kappa _{ψ} }}} \end{aligned}$

(37)

式中，与材料力学参数相关的子函数为：

$\begin{aligned}[b] {C_1} =& \left[ {\displaystyle\frac{{1 - \upsilon - {\kappa _{ψ} }\upsilon }}{{2{G^{\rm{M}}}}} + \displaystyle\frac{{{\kappa _{ψ} } - {\kappa _{ψ} }\upsilon - \upsilon }}{{2{G^{\rm{M}}}}} \cdot \displaystyle\frac{{\csc\; \varphi + 1}}{{\csc\; \varphi - 1}}} \right] \cdot \\ & \quad\quad\quad {\displaystyle\frac{{c\cot \;\varphi }}{{\displaystyle\frac{2}{{\csc\; \varphi - 1}} + {\kappa _{ψ} } + 1}}} \end{aligned}$

(38)

$\begin{aligned}[b] {C_2} =& \left[ {\displaystyle\frac{{1 - \upsilon - {\kappa _{ψ} }\upsilon }}{{2{G^{\rm{M}}}}}c\cot\; \varphi + \frac{{{\kappa _{ψ} } - {\kappa _{ψ} }\upsilon - \upsilon }}{{2{G^{\rm{M}}}}}} \right.\cdot \\ & c\cot\; \varphi +\left. { \displaystyle\frac{{\left( {1 - 2\upsilon } \right)\left( {1 + {\kappa _{ψ} }} \right){p_0}}}{{2{G^{\rm{M}}}}}} \right] \cdot \displaystyle\frac{1}{{{\kappa _{ψ} } + 1}} \end{aligned}$

(39)

将式（16）代入式（37），可得2维平面应变无支护隧道最大径向位移u_r,max。进而联立式（1），可获得开挖面约束作用下各断面的变形情况。式（37）的数学函数形式与Manh解^[5]一致，但由于所考虑问题不一致，两者存在系数项上的差异。当不考虑时间因素时，围岩表现为弹塑性；或当 $\eta^{\rm M}$ →∞，且 $\eta^{\rm K}$ →∞，材料黏弹塑性退化为弹塑性，本文解的数学函数形式与Park解^[21]一致。

3 论证与分析

通过对比已有的数值解和解析解，进而分析隧道时空变形曲线（隧道蠕变特征曲线、围岩变形径向分布曲线与隧道纵剖面变形曲线）关于准则参数c和 $\varphi$ 、扩容角 ${ψ}$ 、延迟时间 $T^{\rm K}$ 等参数的敏感性规律。

3.1 对比分析 3.1.1 与数值解对比

陈炳瑞等^[23]以锦屏Ⅱ级水电站辅助洞断面BK14+190为例，利用FLAC自带的Burgers体和Mohr–Coulomb 组合模型（BMC模型^[23]）计算了该断面时效变形。辅助洞为直壁半圆拱形，埋深1 735 m，地层岩性为大理岩。由于BMC模型与Burgers体和Drucker–Prager组合模型具有一致的力学参数，且Drucker–Prager准则正是Mohr–Coulomb准则的修正模型，二者具有可比性。参数取值^[23]：p₀=36.36 MPa，c=9.61 MPa， $\varphi$ =28.97°， $E^{\rm M}$ =131.19 GPa， $\eta^{\rm M}$ =21 079 GPa·d， $E^{\rm K}$ =51.97 GPa， $\eta^{\rm K}$ =217.21 GPa·d。由于BMC模型不计扩容影响，令 ${ψ}$ =0°， $ \upsilon$ =0.5，进而有G^M=43.73 GPa，G^K=17.32 GPa。直壁半圆拱形隧道近似用等面积圆形隧道等效，其等效半径为3.68 m。图4给出了本文解与BMC模型数值解（平均收敛值）的对比。由图4可知，隧道变形均随时间增长，趋势较为一致，但相比BMC模型数值解，本文计算结果偏小，这可能与直壁半圆拱形的近似等效处理有关。表1给出了围岩变形径向分布和隧道纵剖面变形，其中围岩塑性区半径R_p为8.38 m。由表1可知，围岩变形随着围岩深度的增大而减小，应力释放越充分的远处断面变形较大，与隧道实际变形规律一致。另外，围岩塑性区半径8.38 m与文献[23]计算所得结果较为接近。

3.1.2 与解析解对比

余东明等^[19]以四川省某拟建隧道为例，假设围岩为Burgers体与Drucker–Prager 组合模型，利用平面应变理论计算了圆形隧道黏弹塑性蠕变解。参数取值^[19]：p₀=50 MPa，R₀=6 m，c=6 MPa， $\varphi $ =40°， ${ψ}$ =20°， $E^{\rm M}$ =50 GPa， $\eta^{\rm M}$ = 4 000 GPa·d， $E^{\rm K}$ = 20 GPa， $\eta^{\rm K}$ = 40 GPa·d， $\upsilon$ =0.3。计算可得： $G^{\rm M}$ =19.23 GPa， $G^{\rm K}$ =7.69 GPa。图5给出了2种不同解法结果的对比。表2给出了围岩变形径向分布值和隧道纵剖面变形值，其中围岩塑性区半径R_p为8.03 m。由图5可知：2种解法得到的隧道变形量与变形速率较为一致，但由于余东明等^[19]考虑了支护力，其计算得到的变形值相对于本文解偏小；围岩变形为其深度的递减函数，亦为纵向距离的递增函数，与隧道的实际变形规律是一致的。

通过2个算例的对比分析，可知位移释放系数法成功地描述了岩质隧道施工过程变形的时空效应。基于此，对参数敏感性，包括黏聚力c、内摩擦角 $\varphi$ 、扩容角 ${ψ}$ 和延迟时间 $T^{\rm K}$ 进行研究。

3.2 参数敏感性分析

以四川省某拟建隧道为分析对象^[19]，隧道半径R₀=6 m，c=6 MPa， $\varphi$ =40°， ${ψ}$ =20°， $E^{\rm M}$ =50 GPa， $\eta^{\rm M}$ =4 000 GPa·d， $E^{\rm K}$ =20 GPa， $\eta^{\rm K}$ =40 GPa·d， $\upsilon$ =0.3。隧道最大埋深1 800 m，岩石容重27.5 kN/m³，因而可取地应力p₀=50 MPa。由计算式 $ G^{\rm M}$ = $E^{\rm M}$ /2（1+ $\upsilon $ ）、 $G^{\rm K}$ = $E^{\rm K}$ /2（1+ $\upsilon$ ）可得， $G^{\rm M}$ =19.23 GPa， $G^{\rm K}$ =7.69 GPa。

3.2.1 黏聚力c

图6～8分别为不同黏聚力c时隧道蠕变特征曲线、围岩变形径向分布曲线与隧道纵剖面变形曲线；图7～8中，取蠕变时间为30 d。

总体上，隧道变形为时间t与纵向距离x的递增函数，为围岩深度r的递减函数。由图6～8可知，当黏聚力c线性增大时，隧道时空变形均呈非线性特性。以蠕变时间30 d为例，当c=1 MPa时，洞壁变形量高达273.6 mm，为c=2 MPa时的1.74倍，为c=3 MPa时的2.36倍，为c=4 MPa时的2.92倍，为c=5 MPa时的3.41倍，为c=2 MPa时的3.85倍。这个规律与围岩塑性区半径紧密相关。如表3所示，当黏聚力c线性增大时，塑性区半径非线性递减。黏聚力c的影响规律与其物理意义一致。

隧道开挖面空间效应主要是指隧道纵向的“半圆穹”约束^[8]，也即隧道掌子面前方未开挖岩体对后方卸荷围岩变形发展的约束作用。“半圆穹”是指隧道纵向洞壁变形的形状呈半圆穹形，其约束程度取决于隧道断面形状、围岩特性、地应力、隧道埋深及施工对岩体的扰动程度等因素。由图8可知：当纵向距离小于5～6倍隧道半径时，掌子面约束作用较强；当纵向距离为隧道半径的5～6倍时，断面变形受开挖面的约束作用基本消失。

3.2.2 内摩擦角 $\varphi$

图9～11分别为不同内摩擦角 $ \varphi$ 下隧道蠕变特征曲线、围岩变形径向分布曲线与隧道纵剖面变形曲线；图10～11中，取蠕变时间为30 d。

由图9～11可知，隧道变形为时间t与纵向距离x的递增函数，为围岩深度r的递减函数。当内摩擦角 $\varphi$ 线性增大时，隧道时空变形均呈非线性特性。以蠕变时间30 d为例，当 $\varphi $ =15°时，洞壁变形量高达888.4 mm，相比 $\varphi $ =20°时递增了136%，相比 $\varphi$ =25°时递增了335%，相比 $ \varphi$ =30°时递增了583%，相比 $ \varphi$ =35°时递增了863%，相比 $\varphi$ =40°时递增了1 151%。

表4给出了不同 $\varphi$ 值下围岩塑性区半径，可知当内摩擦角 $\varphi$ 线性增大时，塑性区半径非线性递减。塑性范围的大小在一定程度上决定了隧道变形量，其间存在正相关关系。岩质隧道变形关于内摩擦角 $\varphi$ 的敏感性规律反映了该参数的物理意义。由图11可知：在掌子面附近，“半圆穹”约束作用明显；距开挖面5～6倍隧道半径处的断面，其变形基本不受开挖面的约束作用。结合黏聚力c的敏感性分析，可知这2个材料屈服参数的影响规律极为相似。

3.2.3 扩容角ψ

图12～14分别为不同扩容角ψ下隧道蠕变特征曲线、围岩变形径向分布曲线与隧道纵剖面变形曲线；图13～14中，取蠕变时间为30 d。同样地，隧道变形为时间t与纵向距离x的递增函数，为围岩深度r的递减函数。由图12～14可知：隧道时空变形关于扩容角 ${ψ}$ 的敏感性几乎呈线性特性，随着扩容角 ${ψ}$ 增大，隧道变形几乎线性增加；围岩塑性区半径与扩容角无关。这些规律与扩容角的实际作用一致^[19]。另外，图14表明在距掌子面5～6倍隧道半径范围内，掌子面“半圆穹”约束作用较强。

3.2.4 延迟时间T^K

图15～17分别为不同延迟时间 $ T^{\rm K}(\eta^{\rm K}/G^{\rm K})$ 下隧道蠕变特征曲线、围岩变形径向分布曲线与隧道纵剖面变形曲线；图16～17中，取蠕变时间为30 d。相似地，隧道变形为时间t与纵向距离x的递增函数，为围岩深度r的递减函数。

由图15可知：延迟时间 $T^{\rm K}$ 描述了围岩变形的时程规律，尤其决定了围岩蠕变衰减阶段的时间尺度，但对成洞瞬时变形与最终变形没有影响。当延迟时间 $T^{\rm K}$ 增大时，围岩将在较长时间内处于蠕变衰减阶段。塑性区半径与延迟时间无关，这是因为延迟时间为时效参数。由图17可知，距开挖面5～6倍隧道半径处的断面变形基本不受开挖面的约束作用。

针对四川省该拟建隧道，根据隧道蠕变特征曲线和隧道纵剖面变形曲线，基本可认为成洞20 d以内或距开挖面5～6倍隧道半径以内的卸荷段围岩不宜施作二次支护结构。这是因为成洞20 d以内卸荷段围岩蠕变仍处于衰减阶段，其变形迅速发展；距开挖面5～6倍隧道半径以内的卸荷段围岩变形尚不稳定，相比无穷远处其断面变形仍持续发生。如在施工过程中盲目施作衬砌，有可能导致衬砌结构受力过大，变形突出，最终支护失效酿成工程灾害与经济损失。当然，若在施工过程中预留一定的变形量，推延围岩与支护结构相互作用，便可减小支护结构受力。这种情况下，即使围岩与支护结构共同承载，其间变形也能够相互协调。预留变形量可通过比较支护断面与无穷远断面的时效位移量差值判定。另外，根据围岩变形径向分布曲线，可大致确定诸如锚杆支护深度等支护参数。因此，位移释放系数法可对岩质隧道施工过程提供理论指导。

4 结　论

采用位移释放系数法解析了岩质隧道施工过程的时空变形，并论证了解的正确有效性。最后探讨了隧道时空变形曲线，包括隧道蠕变特征曲线、围岩变形径向分布曲线和隧道纵剖面变形曲线关于黏聚力、内摩擦角、扩容角与延迟时间等参数的敏感性。结果表明，隧道时空变形对介入参数的敏感性与其物理意义一致，位移释放系数法成功描述了岩质隧道施工过程变形的时空效应。

从流变力学角度解释围岩变形的时效性是目前常用的一种方式。作者在解析围岩变形时对隧道力学模型进行了均匀化处理，忽视了岩体结构的非线性特性。在解析围岩变形的时效性问题时，若考虑岩体结构非线性特性，不可避免地要重新评估基本方程（包括平衡方程、变形协调方程）的成立条件。另外，需要一个有充分试验数据支撑的能用于解析推导的变形时效模型。目前，考虑岩体结构非线性特性的流变试验较少开展，建立简单完备的变形时效模型的理论基础也不甚清楚。今后，将开展相应的流变试验研究，并结合试验数据开发完备的变形时效模型，以用于考虑岩体结构非线性特性影响下围岩时效变形问题。