周 建，蔡 键，杨 帆，杨新安

(1.道路与交通工程教育部重点实验室(同济大学)，上海 201804； 2.浙江温州沈海高速公路有限公司，浙江 温州 325000)

随着中国公路、铁路隧道的交通建设蓬勃发展，山岭隧洞穿越富水环境的工况急剧增加[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)初始地应力为p0，rp、r3、r2和r1分别为塑性区半径、衬砌内半径、衬砌外半径和注浆圈外半径，rw为远场水头半径，文献[21]认为隧洞远场水头半径rw一般大于30倍r1时，可保证工程精度，因此，可将隧洞在rw处的径向应力视为p0，围岩远场、弹塑性区边界、注浆圈外边界、衬砌外边界、衬砌内边界的水头分别为hw、hp、h1、h2和h3，围岩、注浆圈和衬砌的渗透系数分别为k1、k2和k3，围岩与注浆圈、注浆圈与衬砌、衬砌内边界和塑性区边界处的压力分别为p1、p2、p3和pp.

图1 力学模型

1.2 渗流场水头分析

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

(1)

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

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

Js=Asvsms.

(2)

式中：s=1、2、3，分别代表围岩、注浆圈和衬砌材料；J1、J2、J3分别为围岩、注浆圈和衬砌内某点水力梯度；A1、m1为围岩的水力梯度待定系数，A2、m2为注浆圈的水力梯度待定系数，A3、m3为衬砌的水力梯度待定系数，其中1﹤m1﹤2时，式(2)反映了显著的惯性效应而导致的非线性渗流特征，当0﹤m1(或m2、m3)﹤1时，式(2)反映了低渗透岩石介质中固液界面效应导致的非线性渗流特征，而当m1(或m2、m3)=1时，式(2)服从达西定律；v1、v2、v3分别为围岩、注浆圈、衬砌内某点渗流速度.

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

(3)

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

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

(4)

式中c1、c2均为待定系数.

当r1

(5)

式(5)为二元一次方程，可解出Asc1ms和c2，代入式(4)继而求出hs. 同理，r2

(6)

1.3 统一强度理论

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

(7)

式中：σ1、σ3为第一主应力、第三主应力；c、φ为材料粘聚力、内摩擦角，ct、φt为材料统一粘聚力、内摩擦角，将围岩、注浆圈和衬砌的粘聚力分别设为c1、c2和c3及内摩擦角分别设为φ1、φ2和φ3，则材料统一粘聚力可对应为ct1、ct2和ct3，材料统一内摩擦角对应为φt1、φt2和φt3；b为参数，反映中间主应力对材料强度的影响程度，取值为0≤b≤1,b的具体值可根据材料的力学试验确定.

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

(8)

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

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

(9)

几何方程为

(10)

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

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

(11)

式中：Ej、μj分别为弹性模量和泊松比，且围岩、注浆圈和衬砌的弹性模量分别为E1、E2、E3，围岩、注浆圈和衬砌的泊松比分别为μ1、μ2、μ3.

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

将式(9)与式(11)联立，并根据边界条件σr|r=rw=-p0，σr|r=r1=-p1，σr|r=r2=-p2，σr|r=rp=-pp,围岩、注浆圈和衬砌有弹性位移分别为

(12)

(13)

(14)

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

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

(15)

(16)

(17)

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

(18)

式中

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

(19)

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

u11|r=r1=u12|r=r2.

(20)

可求得p1和p2. 将衬砌内压p3视为0，由径向应力连续：

σr3|r=rp=σrp4|r=rp.

(21)

可求得rp，从而获得围岩应力，下文参照此法求p1、p2和rp.

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

围岩仍处于弹性状态，则围岩位移u21=u11，围岩径向、切向应力σr21=σr11、σθ21=σθ11. 将式(9)代入式(11)，根据边界条件σr|r=r1=-p1，σr|r=rp=-pp,注浆圈弹性部分位移为

(22)

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

(23)

(24)

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

(25)

(26)

(27)

式中：

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

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

(28)

式中：

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

(29)

(30)

(31)

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

(32)

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

(33)

(34)

式中：

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

3 实例验证

为了验证本文的理论解答，借助文献[15]实例参数,隧洞衬砌内半径r1=4 m、外半径r2=5 m，注浆圈外半径r3=5.2 m,初始地应力p0=10 MPa,远场水头半径rw=200 m，α=1，远场水压力pw=1 MPa,地下水重度γ=10 kN/m3，h1=50 m，h2=h3=0 m，围岩、注浆圈和衬砌参数见表1.

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

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

ksJs=vs，

(35)

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

表1 围岩、注浆圈和衬砌参数

图2 隧道模型图

(36)

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

表2 隧洞洞壁位移对比

表3 隧洞切向应力对比

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

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

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

图3 本文方法得到的塑性区半径与Li解对比图

图4 考虑非线性渗流的隧洞应力分布

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

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

图5 隧洞位移与支护厚度曲线

5 结 论

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

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

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