随着中国公路、铁路隧道的交通建设蓬勃发展，山岭隧洞穿越富水环境的工况急剧增加^[1-3]，此环境中的深埋山岭隧洞主要采取“堵水限排”支护设计准则^[4-5]，一般通过“注浆圈+衬砌”的堵水方式实现. 现阶段，地下水在围岩中的渗流特征尚不明朗，设计人员不能准确把握水压隧洞支护设计^[6]. 文献[7]指出：要在堵水限排情况下进行衬砌结构设计研究，必须从理论上研究水作用下围岩和衬砌结构的力学特征，弄清二者相互作用的机制，才能合理提供隧洞支护设计参数. 关于隧洞在渗流条件下围岩与支护结构的作用理论解析国内外诸多学者进行了研究，如文献[8]基于广义有效应力原理，提出了水压隧洞不同类型衬砌与围岩作用的解析解；文献[9]将数值试验与解析解相结合，为深埋隧道在地下水包围下初期支护与二次衬砌的初步设计提供了合理的途径；文献[10]研究了考虑衬砌和渗流场作用下海底隧道的弹塑性位移和应力解析解；文献[11]推导了围岩、注浆圈、衬砌和地下水共同作用下隧洞的弹塑性解，并提出了最优注浆圈厚度的确定方法；文献[12]基于统一强度理论，给出隧洞渗流压力、衬砌、注浆与岩体的相互作用的弹塑性解；以上支护设计研究均是基于达西定律基础上进行的，部分试验^[13-14]表明：致密砂岩、破碎岩石等介质中的渗流呈现明显的非达西流特征. 为了得到合理的支护设计，必须分析岩土材料的渗流属于线性还是非线性，部分学者对非线性渗流进行了研究，如文献[15]建立二次型高速非达西本构模型，预测了深埋隧洞的涌水量；文献[16]开展了不同颗粒粒径多孔介质在高水力梯度条件下的高速非线性渗流规律试验研究，并确定了非线性渗流模型参数与颗粒粒径之间的关系；文献[17]建立岩体破坏突水非达西渗流模型，模拟了突水瞬态流动全过程，并认为岩体破坏突水问题采用非达西流模型计算十分必要；文献[18]进行了切向位移作用下粗糙单裂隙的高速非达西渗流数值模拟分析，并给出了非线性渗流模型的经验公式；文献[19]建立了粗糙岩石裂隙低速非线性渗流模型并通过试验验证；文献[20]建立了低渗透岩石非线性渗流的运动方程，并通过实验数值验证了所建立运动方程的正确性.

综上所述，岩石渗流可能是低速非线性的、线性的，或是高速非线性的，这取决于地下水在围岩中的分布和隧洞开挖岩体破坏程度. 目前，关于隧洞围岩非线性渗流理论分析鲜少被关注，因此，亟待给出隧洞处于低速非线性渗流或是高速非线性渗流状态下的弹塑性解析解，以便精确地指导隧洞支护设计. 本文将Izbash非线性渗流模型引入隧洞渗流理论，基于统一强度理论，推导深埋隧洞在注浆、衬砌作用下的应力场和位移场解析公式，并讨论了非线性渗流对隧洞围岩应力、位移和塑性区半径的影响.

1 理论基础 1.1 力学模型及基本假定

建立如图 1所示的圆形隧洞力学模型，并作出如下假定：1)深埋隧洞处于地下水包围中；2)围岩为均质﹑各向同性的连续介质；3)水流经围岩和支护材料时流向以径向为主；4)初始地应力为p₀，r_p、r₃、r₂和r₁分别为塑性区半径、衬砌内半径、衬砌外半径和注浆圈外半径，r_w为远场水头半径，文献[21]认为隧洞远场水头半径r_w一般大于30倍r₁时，可保证工程精度，因此，可将隧洞在r_w处的径向应力视为p₀，围岩远场、弹塑性区边界、注浆圈外边界、衬砌外边界、衬砌内边界的水头分别为h_w、h_p、h₁、h₂和h₃，围岩、注浆圈和衬砌的渗透系数分别为k₁、k₂和k₃，围岩与注浆圈、注浆圈与衬砌、衬砌内边界和塑性区边界处的压力分别为p₁、p₂、p₃和p_p.

1.2 渗流场水头分析

本文考虑的非线性渗流即围岩(或注浆圈、衬砌)的水力梯度呈非线性变化，因围岩远场水压与衬砌内水压可确定，则在边界处水头^[22]有

$ {h_s} = \left\{ {\begin{array}{*{20}{l}} {{h_{\rm{w}}}, r = {r_{\rm{w}}};}\\ {{h_1}, r = {r_1};}\\ {{h_2}, r = {r_2};}\\ {{h_3}, r = {r_3}.} \end{array}} \right. $

(1)

式中：h_s为不同位置处的水头，h_s=p_s/γ，其中p_s为不同位置处的水压力，γ为地下水重度. 水头一般通过监测水压并根据水头、水压的关系式获得.

幂函数型Izbash方程因其公式简单明了被广泛应用于研究岩土材料的非线性渗流特征，Izbash方程的水力梯度公式^[23]为

$ {J_s} = {A_s}{v_s}^{{m_s}}. $

(2)

式中：s=1、2、3，分别代表围岩、注浆圈和衬砌材料；J₁、J₂、J₃分别为围岩、注浆圈和衬砌内某点水力梯度；A₁、m₁为围岩的水力梯度待定系数，A₂、m₂为注浆圈的水力梯度待定系数，A₃、m₃为衬砌的水力梯度待定系数，其中1<m₁<2时，式(2)反映了显著的惯性效应而导致的非线性渗流特征，当0<m₁(或m₂、m₃)<1时，式(2)反映了低渗透岩石介质中固液界面效应导致的非线性渗流特征，而当m₁(或m₂、m₃)=1时，式(2)服从达西定律；v₁、v₂、v₃分别为围岩、注浆圈、衬砌内某点渗流速度.

根据文献[24]，平面径向渗流速度方程有

$ \frac{\partial v_{s}}{\partial r}+\frac{v_{s}}{r}=0 . $

(3)

式中r为任意一点距离洞心的距离.

求解式(3)，代入式(2)后等式两边积分，并根据J_s=∂h_s/∂r，可得

$ h_{s}=\frac{A_{s} {c_{1}}^{m_{s}}}{1-m_{s}} r^{1-m_{s}}+c_{2}. $

(4)

式中c₁、c₂均为待定系数.

当r₁ < r < r_w时，由式(1)中边界条件h_s|_r=rw =h_w，h_s|_r=r1 =h₂，有

$ \left\{\begin{array}{l} \frac{A_{s} c_{1}{ }^{m_{s}}}{1-m_{s}} r_{\mathrm{w}}^{1-m_{s}}+c_{2}=h_{\mathrm{w}}, \\ \frac{A_{s} c_{1}{ }^{m_{s}}}{1-m_{s}} r_{1}^{1-m_{s}}+c_{2}=h_{2}. \end{array}\right. $

(5)

式(5)为二元一次方程，可解出A_sc₁^m_s和c₂，代入式(4)继而求出h_s. 同理，r₂ < r < r₁和r₃ < r < r₂时的h_s解法参照上述方法. 则有

$ h_{s}=\left\{\begin{array}{l} \frac{h_{1}\left(r_{\mathrm{w}}^{1-m_{1}}-r^{1-m_{1}}\right)+h_{\mathrm{w}}\left(r^{1-m_{1}}-r_{1}^{1-m_{1}}\right)}{r_{\mathrm{w}}^{1-m_{1}}-r_{1}^{1-m_{1}}}, r_{1}<r<r_{w}; \\ \frac{h_{2}\left(r_{1}^{1-m_{2}}-r^{1-m_{2}}\right)+h_{1}\left(r^{1-m_{2}}-r_{2}^{1-m_{1}}\right)}{r_{1}^{1-m_{2}}-r_{2}^{1-m_{2}}}, r_{2}<r<r_{1}; \\ \frac{h_{3}\left(r_{2}^{1-m_{1}}-r^{1-m_{2}}\right)+h_{2}\left(r^{1-m_{3}}-r_{3}^{1-m_{3}}\right)}{r_{2}^{1-m_{3}}-r_{3}^{1-m_{3}}}, r_{3}<r<r_{2}. \end{array}\right. $

(6)

1.3 统一强度理论

统一强度理论自俞茂宏创立以来，广泛应用于岩石、混凝土等拉压特性不同的材料，弥补了Mohr-Coulomb准则未能考虑中间主应力致使计算结果偏于保守的遗憾. 统一强度理论在平面应变状态下的表达式为

$ \left\{\begin{array}{l} \frac{\sigma_{1}-\sigma_{3}}{2}=\frac{\sigma_{1}+\sigma_{3}}{2} \sin \varphi_{t}+c_{t} \cos \varphi_{t}, \\ \sin \varphi_{t}=\frac{2(1+b) \sin \varphi}{1+b(1+\sin \varphi)}, \\ c_{t}=\frac{2(1+b) c \cos \varphi}{2+b(1+\sin \varphi)} \frac{1}{\cos \varphi_{t}}. \end{array}\right. $

(7)

式中：σ₁、σ₃为第一主应力、第三主应力；c、φ为材料粘聚力、内摩擦角，c_t、φ_t为材料统一粘聚力、内摩擦角，将围岩、注浆圈和衬砌的粘聚力分别设为c₁、c₂和c₃及内摩擦角分别设为φ₁、φ₂和φ₃，则材料统一粘聚力可对应为c_t1、c_t2和c_t3，材料统一内摩擦角对应为φ_t1、φ_t2和φ_t3；b为参数，反映中间主应力对材料强度的影响程度，取值为0≤b≤1, b的具体值可根据材料的力学试验确定.

假设σ_θ(径向有效应力) > σ_r(切向有效应力)，令σ₁=σ_θ、σ₃=σ_r，则式(7)中第1式可改写为

$ \frac{{{\sigma _\theta } - {\sigma _r}}}{2} = \frac{{{\sigma _\theta } + {\sigma _r}}}{2}\sin {\varphi _{\rm{t}}} + {c_{\rm{t}}}\cos {\varphi _{\rm{t}}}. $

(8)

2 非线性渗流作用下的隧洞围岩应力、位移分析

根据弹性力学理论，材料考虑渗透力时的平衡微分方程为

$ \frac{\mathrm{d} \sigma_{r}}{\mathrm{~d} r}+\frac{\sigma_{r}-\sigma_{\theta}}{r}+F_{\mathrm{w}}=0. $

(9)

式中：σ_r、σ_θ分别为材料径向与切向有效应力；F_w为渗透力，$F_{\mathrm{w}}=-\alpha \gamma \frac{\mathrm{d} h}{\mathrm{~d} r} $，其中α为材料等效孔隙水压力系数，由于岩石的孔隙特性不同于松散体介质，太沙基的有效孔隙应力原理需修正^[25].

几何方程为

$ \left\{\begin{array}{l} \varepsilon_{\theta}=\frac{u}{r}, \\ \varepsilon_{r}=\frac{\mathrm{d} u}{\mathrm{~d} r}. \end{array}\right. $

(10)

式中u为材料的径向位移.

弹性区服从虎克定律的平面应力应变方程为

$ \left\{\begin{array}{l} \sigma_{r}=\frac{E_{j}\left(1-\mu_{j}\right)}{\left(1+\mu_{j}\right)\left(1-2 \mu_{j}\right)}\left(\frac{\mathrm{d} u}{\mathrm{~d} r}+\frac{\mu_{j}}{1-\mu_{j}} \frac{u}{r}\right), \\ \sigma_{\theta}=\frac{E_{j}\left(1-\mu_{j}\right)}{\left(1+\mu_{j}\right)\left(1-2 \mu_{j}\right)}\left(\frac{u}{r}+\frac{\mu_{j}}{1-\mu_{j}} \frac{\mathrm{d} u}{\mathrm{~d} r}\right). \end{array}\right. $

(11)

式中：E_j、μ_j分别为弹性模量和泊松比，且围岩、注浆圈和衬砌的弹性模量分别为E₁、E₂、E₃，围岩、注浆圈和衬砌的泊松比分别为μ₁、μ₂、μ₃.

2.1 弹塑性交界处位于衬砌内

将式(9)与式(11)联立，并根据边界条件${\left. {{\sigma _r}} \right|_{r = {r_{\rm{w}}}}} = - {p_0}, {\left. {{\sigma _r}} \right|_{r = {r_1}}} = - {p_1}, $${\left. {{\sigma _r}} \right|_{r = {r_2}}} = - {p_2}, {\left. {{\sigma _r}} \right|_{r = {r_{\rm{p}}}}} = - {p_{\rm{p}}} $, 围岩、注浆圈和衬砌有弹性位移分别为

$ u_{11}=\frac{\left(1+\mu_{1}\right)\left(1-2 \mu_{1}\right)}{E_{1}}\left(a_{11} r+\frac{b_{11}}{r}+\lambda_{11} r^{2-m_{1}}\right), $

(12)

$ u_{12}=\frac{\left(1+\mu_{2}\right)\left(1-2 \mu_{2}\right)}{E_{2}}\left(a_{12} r+\frac{b_{12}}{r}+\lambda_{12} r^{2-m_{2}}\right), $

(13)

$ u_{13}=\frac{\left(1+\mu_{3}\right)\left(1-2 \mu_{3}\right)}{E_{3}}\left(a_{13} r+\frac{b_{13}}{r}+\lambda_{13} r^{2-m_{3}}\right) . $

(14)

式中：$a_{11}=-p_{0}-\lambda_{11}\left(2-m_{1}+\frac{\mu_{1}}{1-\mu_{1}}\right) r_{\mathrm{w}}^{1-m_{1}}, $

$ b_{11}=\lambda_{11}\left(2-m_{1}+\frac{\mu_{1}}{1-\mu_{1}}\right)\left(r_{1}^{1-m_{1}}-r_{\mathrm{w}}^{1-m_{1}}\right) , $

$ \lambda_{11}=\frac{\alpha r_{\mathrm{w}}\left(h_{1}-h_{\mathrm{w}}\right)}{\left(1-\mu_{1}\right)\left(m_{1}-3\right)\left[\left(r_{\mathrm{w}}^{1-m_{1}}-r_{1}^{1-m_{1}}\right)\right]}, $

$ a_{12}=-p_{1}-\lambda_{12}\left(2-m_{2}+\frac{\mu_{2}}{1-\mu_{2}}\right) r_{1}^{1-m_{2}}, $

$ b_{12}=\lambda_{12}\left(2-m_{2}+\frac{\mu_{2}}{1-\mu_{2}}\right)\left(r_{2}^{1-m_{2}}-r_{1}^{1-m_{2}}\right), $

$ \lambda_{12}=\frac{\alpha r_{\mathrm{w}}\left(h_{2}-h_{1}\right)}{\left(1-\mu_{2}\right)\left(m_{2}-3\right)\left[\left(r_{1}^{1-m_{2}}-r_{2}^{1-m_{2}}\right)\right]}, $

$ a_{13}=-p_{2}-\lambda_{13}\left(2-m_{3}+\frac{\mu_{3}}{1-\mu_{3}}\right) r_{2}^{1-m_{3}}, $

$ b_{13}=\lambda_{13}\left(2-m_{3}+\frac{\mu_{3}}{1-\mu_{3}}\right)\left(r_{\mathrm{p}}^{1-m_{3}}-r_{2}^{1-m_{3}}\right), $

$ \lambda_{13}=\frac{\alpha r_{\mathrm{w}}\left(h_{\mathrm{p}}-h_{2}\right)}{\left(1-\mu_{3}\right)\left(m_{3}-3\right)\left[\left(r_{2}^{1-m_{3}}-r_{\mathrm{p}}^{1-m_{3}}\right)\right]} . $

其中h_p在所处材料中按照差值法计算.

将式(12)~(14)代入式(11)，则围岩、注浆圈和衬砌径向、切向应力分别为

$ \left\{\begin{array}{l} \sigma_{r 11}=a_{11}+b_{11}\left(\frac{r_{1}}{r}\right)^{2}+\lambda_{11} r^{1-m_{1}}\left[\left(2-m_{1}\right)\left(1-\mu_{1}\right)+\mu_{1}\right], \\ \sigma_{\theta 11}=a_{12}-b_{12}\left(\frac{r_{1}}{r}\right)^{2}+\lambda_{11} r^{1-m_{1}}\left[1+\mu_{1}\left(1-m_{1}\right)\right] . \end{array}\right. $

(15)

$ \left\{\begin{array}{l} \sigma_{r \mathrm{12}}=a_{\mathrm{12}}+b_{\mathrm{12}}\left(\frac{r_{2}}{r}\right)^{2}+\lambda_{\mathrm{12}} r^{1-m_{2}}\left[\left(2-m_{2}\right)\left(1-\mu_{2}\right)+\mu_{2}\right], \\ \sigma_{\theta 12}=a_{12}-b_{12}\left(\frac{r_{2}}{r}\right)^{2}+\lambda_{12} r^{1-m_{2}}\left[1+\mu_{2}\left(1-m_{2}\right)\right]. \end{array}\right. $

(16)

$ \left\{\begin{array}{l} \sigma_{r 13}=a_{13}+b_{13}\left(\frac{r_{\mathrm{p}}}{r}\right)^{2}+\lambda_{13} r^{1-m_{3}}\left[\left(2-m_{3}\right)\left(1-\mu_{3}\right)+\mu_{3}\right], \\ \sigma_{\theta 13}=a_{13}-b_{13}\left(\frac{r_{\mathrm{p}}}{r}\right)^{2}+\lambda_{13} r^{1-m_{3}}\left[1+\mu_{3}\left(1-m_{3}\right)\right]. \end{array}\right. $

(17)

将式(8)与式(9)联立，并根据边界条件σ_r|_r=r3 =-p₃，则衬砌内塑性径向、切向应力分别为

$ \left\{\begin{array}{l} \sigma_{r \mathrm{p} 14}=c_{\mathrm{t} 3} \cot \varphi_{\mathrm{t} 3}+B_{1} r^{1-m_{3}}-\left(p_{3}+c_{\mathrm{t} 3} \cot \varphi_{\mathrm{t} 3}+\right. \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \left.B_{1} r_{3}^{1-m_{3}}\right)\left(\frac{r}{r_{3}}\right)^{-\frac{2 \sin \varphi_{\mathrm{t} 3}}{1-\sin \varphi_{\mathrm{t} 3}}}, \\ \sigma_{\theta \mathrm{p} 14}=-\frac{1+\sin \varphi_{\mathrm{t} 3}}{1-\sin \varphi_{\mathrm{t} 3}}\left[\left(p_{3}+c_{\mathrm{t} 3} \cot \varphi_{\mathrm{t} 3}+B_{1} r_{3}^{1-m_{3}}\right)\right. \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \left.\left(\frac{r}{r_{3}}\right)^{-\frac{2 \sin \varphi_{\mathrm{t} 3}}{1-\sin \varphi_{\mathrm{t} 3}}}-B_{1} r^{1-m_{3}}\right]+c_{\mathrm{t} 3} \cot \varphi_{\mathrm{t} 3} . \end{array}\right. $

(18)

式中

$ {B_1} = \frac{{{k_3}{r_{\rm{w}}}\left( {1 + \sin {\varphi _{{\rm{t3}}}}} \right)\left( {1 - {m_3}} \right)\left( {{h_3} - {h_2}} \right)}}{{\left[ {3\sin {\varphi _{{\rm{t3}}}} + 1 - {m_3}\left( {1 + \sin {\varphi _{{\rm{t3}}}}} \right)} \right]\left( {r_2^{1 - {m_3}} - {r_3}} \right)}}. $

假定材料塑性阶段体积应变为0，根据几何方程式(10)，则衬砌塑性区位移为

$ u_{14}=\left.\frac{r_{\mathrm{p}}}{r} u_{13}\right|_{r=r_{\mathrm{p}}} . $

(19)

联立式(12)~(14)，并根据位移在边界处连续：

$ \left.u_{11}\right|_{r=r_{1}}=\left.u_{12}\right|_{r=r_{2}}. $

(20)

可求得p₁和p₂. 将衬砌内压p₃视为0，由径向应力连续：

$ \left.\sigma_{r 3}\right|_{r=r_{\mathrm{p}}}=\left.\sigma_{r \mathrm{p} 4}\right|_{r=r_{\mathrm{p}}}. $

(21)

可求得r_p，从而获得围岩应力，下文参照此法求p₁、p₂和r_p.

2.2 弹塑性交界处位于注浆圈内

围岩仍处于弹性状态，则围岩位移u₂₁=u₁₁，围岩径向、切向应力σ_r21=σ_r11、σ_θ21=σ_θ11. 将式(9)代入式(11)，根据边界条件σ_r|_r=r1=-p₁，σ_r|_r=rp=-p_p, 注浆圈弹性部分位移为

$ u_{22}=\frac{\left(1+\mu_{2}\right)\left(1-2 \mu_{2}\right)}{E_{2}}\left(a_{22} r+\frac{b_{22}}{r}+\lambda_{22} r^{2-m_{2}}\right) . $

(22)

式中：$a_{22}=-p_{1}-\lambda_{22}\left(2-m_{2}+\frac{\mu_{2}}{1-\mu_{2}}\right) r_{1}^{1-m_{2}}, $

$ b_{22}=\lambda_{22}\left(2-m_{2}+\frac{\mu_{2}}{1-\mu_{2}}\right)\left(r_{\mathrm{p}}^{1-m_{2}}-r_{1}^{1-m_{2}}\right), $

$ \lambda_{22}=\frac{\alpha r_{\mathrm{w}}\left(h_{\mathrm{p}}-h_{1}\right)}{\left(1-\mu_{2}\right)\left(m_{2}-3\right)\left[\left(r_{1}^{1-m_{2}}-r_{\mathrm{p}}^{1-m_{2}}\right)\right]} . $

材料塑性阶段位移参照式(19), 则注浆圈、衬砌塑性部分位移分别为

$ u_{23}=\left.\frac{r_{\mathrm{p}}}{r} u_{22}\right|_{r=r_{\mathrm{p}}}, $

(23)

$ u_{24}=\left.\frac{r_{2}}{r} u_{23}\right|_{r=r_{2}}. $

(24)

将式(22)代入式(11)，将式(8)代入式(9)，并根据边界条件σ_r|_r=r2 =-p₂、σ_r|_r=r3 =-p₃，注浆圈弹性部分、注浆圈塑性部分和衬砌的径向、切向应力分别为

$ \left\{\begin{array}{l} \sigma_{r 22}=a_{22}+b_{22}\left(\frac{r_{2}}{r}\right)^{2}+\lambda_{22} r^{1-m_{2}}\left[\left(2-m_{2}\right)\left(1-\mu_{2}\right)+\mu_{2}\right], \\ \sigma_{\theta 22}=a_{22}-b_{22}\left(\frac{r_{2}}{r}\right)^{2}+\lambda_{22} r^{1-m_{2}}\left[1+\mu_{2}\left(1-m_{2}\right)\right] . \end{array}\right. $

(25)

$ \left\{ {\begin{array}{*{20}{l}} {{\sigma _{r{\rm{p}}23}} = {c_{{\rm{t2}}}}\cot {\varphi _{{\rm{t2}}}} + {B_2}{r^{1 - {m_2}}} - \left( {{p_2} + {c_{{\rm{t2}}}}\cot {\varphi _{{\rm{t2}}}} + {B_2}r_2^{1 - {m_2}}} \right){{\left( {\frac{r}{{{r_2}}}} \right)}^{ - \frac{{2\sin {\varphi _{{\rm{t2}}}}}}{{1 - \sin {\varphi _{{\rm{t2}}}}}}}}, }\\ {{\sigma _{{\theta _{{\rm{p}}23}}}} = - \frac{{1 + \sin {\varphi _{{\rm{t2}}}}}}{{1 - \sin {\varphi _{{\rm{t2}}}}}}\left[ {\left( {{p_2} + {c_{{\rm{t2}}}}\cot {\varphi _{{\rm{t2}}}} + {B_2}r_2^{1 - {m_2}}} \right){{\left( {\frac{r}{{{r_3}}}} \right)}^{ - \frac{{2\sin {\varphi _{{\rm{t2}}}}}}{{1 - \sin {\varphi _{{\rm{t2}}}}}}}} - {B_2}{r^{1 - {m_2}}}} \right] + {c_{{\rm{t2}}}}\cot {\varphi _{{\rm{t2}}}}.} \end{array}} \right. $

(26)

$ \left\{ {\begin{array}{*{20}{l}} {{\sigma _{r{\rm{p}}24}} = {c_{{\rm{t}}3}}\cot {\varphi _{{\rm{t}}3}} + {B_3}{r^{1 - {m_3}}} - \left( {{p_3} + {c_{{\rm{t}}3}}\cot {\varphi _{{\rm{t}}3}} + {B_3}r_3^{1 - {m_3}}} \right){{\left( {\frac{r}{{{r_3}}}} \right)}^{ - \frac{{2\sin {\varphi _{{\rm{t}}3}}}}{{1 - \sin {\varphi _{{\rm{t}}3}}}}}}, }\\ {{\sigma _{\theta {\rm{p}}24}} = - \frac{{1 + \sin {\varphi _{{\rm{t}}3}}}}{{1 - \sin {\varphi _{{\rm{t}}3}}}}\left[ {\left( {{p_3} + {c_{{\rm{t}}3}}\cot {\varphi _{{\rm{t}}3}} + {B_3}r_3^{1 - {m_3}}} \right){{\left( {\frac{r}{{{r_3}}}} \right)}^{ - \frac{{2\sin {\varphi _{{\rm{t}}3}}}}{{1 - \sin {\varphi _{{\rm{t}}3}}}}}} - {B_3}{r^{1 - {m_3}}}} \right] + {c_{{\rm{t}}3}}\cot {\varphi _{{\rm{t}}3}}.} \end{array}} \right. $

(27)

式中：

$ {{B_2} = \frac{{{k_2}{r_{\rm{w}}}\left( {1 + \sin {\varphi _{{\rm{t2}}}}} \right)\left( {1 - {m_2}} \right)\left( {{h_2} - {h_1}} \right)}}{{\left[ {3\sin {\varphi _{{\rm{t2}}}} + 1 - {m_2}\left( {1 + \sin {\varphi _{12}}} \right)} \right]\left( {r_1^{1 - {m_2}} - {r_2}} \right)}}, } $

$ {{B_3} = \frac{{{k_3}{r_{\rm{w}}}\left( {1 + \sin {\varphi _{{\rm{t3}}}}} \right)\left( {1 - {m_3}} \right)\left( {{h_3} - {h_2}} \right)}}{{\left[ {3\sin {\varphi _{{\rm{t3}}}} + 1 - {m_3}\left( {1 + \sin {\varphi _{{\rm{t3}}}}} \right)} \right]\left( {r_2^{1 - {m_3}} - {r_3}} \right)}}.} $

2.3 弹塑性交界处位于围岩内

将式(9)与式(11)联立，并根据边界条件σ_r|_r=rw =-p₀，σ_r|_r=rp =-p_p，围岩弹性位移有

$ {u_{31}} = \frac{{\left( {1 + {\mu _1}} \right)\left( {1 - 2{\mu _1}} \right)}}{{{E_1}}}\left( {{a_{31}}r + \frac{{{b_{31}}}}{r} + {\lambda _{31}}{r^{2 - {m_1}}}} \right). $

(28)

式中：

$ {{a_{31}} = - {p_0} - {\lambda _{31}}\left( {2 - {m_1} + \frac{{{\mu _1}}}{{1 - {\mu _1}}}} \right)r_{\rm{w}}^{1 - {m_1}}, } $

$ {{b_{31}} = {\lambda _{31}}\left( {2 - {m_1} + \frac{{{\mu _1}}}{{1 - {\mu _1}}}} \right)\left( {r_{\rm{p}}^{1 - {m_1}} - r_{\rm{w}}^{1 - {m_1}}} \right), } $

$ {{\lambda _{31}} = \frac{{\alpha {r_{\rm{w}}}\left( {{h_{\rm{p}}} - {h_{\rm{w}}}} \right)}}{{\left( {1 - {\mu _1}} \right)\left( {{m_1} - 3} \right)\left[ {\left( {r_{\rm{w}}^{1 - {m_1}} - r_{\rm{p}}^{1 - {m_1}}} \right)} \right]}}.} $

材料塑性阶段位移参照式(19), 则围岩、注浆圈和衬砌塑性部分位移分别为

$ {{u_{32}} = {{\left. {\frac{{{r_{\rm{p}}}}}{r}{u_{31}}} \right|}_{r = {r_{\rm{p}}}}}, } $

(29)

$ {{u_{33}} = {{\left. {\frac{{{r_1}}}{r}{u_{32}}} \right|}_{r = {r_1}}}, } $

(30)

$ {{u_{34}} = {{\left. {\frac{{{r_2}}}{r}{u_{33}}} \right|}_{r = {r_2}}}.} $

(31)

将式(28)代入式(11)，围岩弹性部分径向、切向应力分别为

$ \left\{ {\begin{array}{*{20}{l}} {{\sigma _{r31}} = {a_{31}} + {b_{31}}{{\left( {\frac{{{r_{\rm{p}}}}}{r}} \right)}^2} + {\lambda _{31}}{r^{1 - {m_1}}}\left[ {\left( {2 - {m_1}} \right)\left( {1 - {\mu _1}} \right) + {\mu _1}} \right], }\\ {{\sigma _{\theta 31}} = {a_{31}} - {b_{31}}{{\left( {\frac{{{r_{\rm{p}}}}}{r}} \right)}^2} + {\lambda _{31}}{r^{1 - {m_1}}}\left[ {1 + {\mu _1}\left( {1 - {m_1}} \right)} \right].} \end{array}} \right. $

(32)

将式(8)代入式(9)，并根据边界条件σ_r|_r=r2 =-p₂、σ_r|_r=r3 =-p₃，围岩塑性部分和注浆圈的径向、切向应力分别为

$ \left\{ {\begin{array}{*{20}{l}} {{\sigma _{r{\rm{p}}32}} = {c_{{\rm{t1}}}}\cot {\varphi _{{\rm{t1}}}} + {B_4}{r^{1 - {m_1}}} - \left( {{p_1} + {c_{{\rm{t1}}}}\cot {\varphi _{{\rm{t1}}}} + {B_4}r_1^{1 - {m_1}}} \right){{\left( {\frac{r}{{{r_1}}}} \right)}^{ - \frac{{2\sin {\varphi _{{\rm{t1}}}}}}{{1 - \sin {\varphi _{{\rm{t1}}}}}}}}, }\\ {{\sigma _{\theta {\rm{p}}32}} = - \frac{{1 + \sin {\varphi _{{\rm{t1}}}}}}{{1 - \sin {\varphi _{{\rm{t1}}}}}}\left[ {\left( {{p_1} + {c_{{\rm{t1}}}}\cot {\varphi _{{\rm{t1}}}} + {B_4}r_1^{1 - {m_1}}} \right){{\left( {\frac{r}{{{r_1}}}} \right)}^{ - \frac{{2\sin {\varphi _{{\rm{t1}}}}}}{{1 - \sin {\varphi _{{\rm{t1}}}}}}}} - {B_4}{r^{1 - {m_1}}}} \right] + {c_{{\rm{t1}}}}\cot {\varphi _{{\rm{t1}}}};} \end{array}} \right. $

(33)

$ \left\{ {\begin{array}{*{20}{l}} {{\sigma _{r{\rm{p}}33}} = {c_{{\rm{t}}2}}\cot {\varphi _{{\rm{t}}2}} + {B_5}{r^{1 - {m_2}}} - \left( {{p_2} + {c_{{\rm{t}}2}}\cot {\varphi _{{\rm{t}}2}} + {B_5}r_2^{1 - {m_2}}} \right){{\left( {\frac{r}{{{r_2}}}} \right)}^{ - \frac{{2\sin {\varphi _{{\rm{t}}2}}}}{{1 - \sin {\varphi _{{\rm{t}}2}}}}}}, }\\ {{\sigma _{\theta {\rm{p}}33}} = - \frac{{1 + \sin {\varphi _{{\rm{t}}2}}}}{{1 - \sin {\varphi _{{\rm{t}}2}}}}\left[ {\left( {{p_2} + {c_{{\rm{t}}2}}\cot {\varphi _{{\rm{t}}2}} + {B_5}r_2^{1 - {m_2}}} \right){{\left( {\frac{r}{{{r_2}}}} \right)}^{ - \frac{{2\sin {\varphi _{{\rm{t}}2}}}}{{1 - \sin {\varphi _{{\rm{t}}2}}}}}} - {B_5}{r^{1 - {m_2}}}} \right] + {c_{{\rm{t}}2}}\cot {\varphi _{{\rm{t}}2}}.} \end{array}} \right. $

(34)

式中：

$ {{B_4} = \frac{{{k_1}{r_{\rm{w}}}\left( {1 + \sin {\varphi _{{\rm{t}}1}}} \right)\left( {1 - {m_1}} \right)\left( {{h_1} - {h_{\rm{p}}}} \right)}}{{\left[ {3\sin {\varphi _{{\rm{t}}1}} + 1 - {m_1}\left( {1 + \sin {\varphi _{{\rm{t}}1}}} \right)} \right]\left( {r_{\rm{p}}^{1 - {m_1}} - {r_1}} \right)}}, } $

$ {{B_5} = \frac{{{k_2}{r_{\rm{w}}}\left( {1 + \sin {\varphi _{{\rm{t2}}}}} \right)\left( {1 - {m_2}} \right)\left( {{h_2} - {h_1}} \right)}}{{\left[ {3\sin {\varphi _{{\rm{t2}}}} + 1 - {m_2}\left( {1 + \sin {\varphi _{{\rm{t2}}}}} \right)} \right]\left( {r_1^{1 - {m_2}} - {r_2}} \right)}}.} $

衬砌的径向、切向应力与式(27)相同，即σ_rp34= σ_rp24、σ_θp34=σ_θp24.

3 实例验证

为了验证本文的理论解答，借助文献[15]实例参数, 隧洞衬砌内半径r₁=4 m、外半径r₂=5 m，注浆圈外半径r₃=5.2 m, 初始地应力p₀=10 MPa, 远场水头半径r_w=200 m，α=1，远场水压力p_w=1 MPa, 地下水重度γ=10 kN/m³，h₁=50 m，h₂=h₃=0 m，围岩、注浆圈和衬砌参数见表 1.

理论计算时，将隧道已知参数代入上述公式，当求得塑性区半径r_p与某一种弹塑性边界情况一致时，则位移场与应力场按照该情况计算.

利用ABAQUS有限元软件模拟文献[15]工况，为了减少模型边界效应，建立模型如图 3所示，模型为210 m×210 m监测断面为隧洞纵向中心断面. 为模拟隧洞的非线性渗流，根据地下水流速公式：

$ {k_s}{J_s} = {v_s}, $

(35)

式中k_s为材料的渗透系数，并结合式(2)可得

$ {k_s}{A_s} = v_s^{1 - {m_s}}. $

(36)

假设在隧洞远端的渗流速度v_s和水力梯度参数A_s一定，且取v_s为1×10^-4 m/s，根据式(36)，可得m₁=0.5、m₁=1.5的围岩等效渗透系数分别为5×10^-6 m/s、0.05 m/s，则隧洞监测断面洞壁的位移、切向应力的数值解与本文理论解对比见表 2、3.

从表 2可以看出，本文理论解较数值解略大，随着m₁增大，洞壁位移误差在8.8%以内，隧洞切向应力误差最大为14.1%，理论解与数值解吻合性较好，验证了理论解的可靠性.

4 非线性渗流对隧洞塑性区半径、围岩应力和位移的影响

一般来说，材料的渗透系数越小，水流穿过材料裂隙的难度越大，水力梯度待定系数越小，由于表 1中注浆圈和衬砌的渗透系数远小于围岩，不难判断，m₁值大于m₂和m₃，假设m₂=m₃，分析非线性渗流对隧洞塑性区半径、围岩应力和位移的影响，如图 3、4、5所示.

从图 3可以看出，随着b值增大，塑性区半径r_p逐渐减小且最大减小幅度为13.2%，而工程实践中r_p的精确性直接影响支护结构设计，说明考虑中间主应力是必要的，否则误差较大；m₁对r_p值的影响远大于m₂和m₃，且m₁越大r_p越大. 根据上述分析，图 4仅探讨m₁对围岩应力的影响，因此假定m₂=m₃=0.4. 当b=0.4时，随着m₁越大，围岩径向、切向应力越小. 结合图 3和图 4，当m₁=1.0、m₂=m₃=0.4时，围岩渗流服从达西定律，塑性区半径、围岩应力与文献[15]较为接近，从而进一步验证了本文理论方法的可行性. 综上所述，当围岩渗流为高渗透性非线性渗流时，若采用达西定律则低估塑性区半径、高估实际围岩应力，不利于支护设计；围岩渗流为低渗透性非线性渗流时，达西定律计算得到塑性区半径偏大、围岩应力偏保守.

图 5显示了b=0.4时隧洞位移与支护厚度的关系，D₁、D₂分别为注浆圈、衬砌厚度，当支护厚度一定时，围岩渗流从低渗透性非线性渗流向高渗透性非线性渗流转变时，隧洞位移逐渐增大，因此，当隧洞围岩处于高水压、突水严重地层时，传统的达西定律计算的洞壁位移偏于保守，不利于隧洞施工安全. 另外，随着注浆圈和衬砌厚度的增大，洞壁位移能够得到有效控制，但衬砌厚度超过1 m时，洞壁位移的控制并不显著，则支护设计时应合理设定支护材料厚度，以免造成较大浪费.

5 结论

1) 在注浆、衬砌支护条件下，基于Izbash非线性渗流模型，求得围岩及支护材料的水头分布，并运用统一强度理论，推导了塑性区在不同材料位置时的应力和位移弹塑性解.

2) 考虑非线性渗流的理论位移解较等效后的ABAQUS解略大，切向应力误差在允许范围之内，吻合性较好；中间主应力和围岩的水力梯度系数m₁对隧洞塑性区半径的影响显著，而对围岩应力的影响较塑性区半径略小.

3) 围岩水力梯度系数m₁对隧洞位移的影响不容忽视，当隧洞围岩处于高速非线性渗流时，达西定律计算的隧洞位移偏保守，体现了研究非线性渗流作用下隧洞弹塑性分析的重要意义；另外，不能盲目追求安全性而增加支护材料厚度，否则导致支护不经济.