考虑土体强度非均质和各向异性的隧道洞口仰坡地震稳定性上限解

张京伍, 荣耀, 李栋伟, 孙洋, 何昌春, 陈欣

(1.江西省交通科学研究院有限公司， 南昌 330200；

2. 东华理工大学土木与建筑工程学院， 南昌 330013；
3.苏州科技大学土木工程学院， 苏州 215011 )

山区交通工程建设中往往形成隧道洞口仰坡，隧道开挖施工会造成仰坡土体产生巨大扰动，破坏边坡土体原有的应力平衡状态，导致坡体抗剪强度下降，极易诱发滑坡、崩塌等工程地质问题[1-4]。隧道洞口区域岩土体风化程度较大，风化裂隙发育且受地表水侵蚀严重，该区域的衬砌结构是失效破坏的多发地段，特别是当隧道洞口支护不及时或衬砌支护能力不足的情况下，容易发生衬砌变形破坏及支护拱脱离等现象。隧道洞口段稳定与边坡工程密切相关[5-6]，仰坡稳定对隧道结构的安全起到至关重要的作用。地震诱发的边坡滑动是主要地震地质灾害类型之一[7]，其对边坡稳定性的影响机制一直都是土力学研究的热点和难点[8-11]。因此，如何分析地震力作用下隧道洞口仰坡稳定性是研究隧道工程建设及运营期安全性能的关键问题。

大量试验研究表明，受开挖、天然沉积及后期加载等因素影响，岩土体黏聚力往往表现出显著的非均质性和各向异性[12-15]，且边坡表面分布着大量裂缝[16]。忽略岩土材料的各向异性和非均质性会导致边坡稳定性分析产生误差，从而影响含裂缝边坡稳定性评价的准确性[17]。在中国边坡工程相关设计规范的稳定性计算中多采用极限平衡理论，针对边坡失稳的圆弧滑动面多采用圆弧滑动法、Bishop法及Morgenstern-Price 法，针对非圆弧滑面多采用不平衡推力法。然而，极限平衡法所得到的解既不是真实解的上限，也不是真实解的下限[18-19]。基于土体塑性力学理论，极限分析法首次由Drucker等[20]于1952年引入边坡稳定性分析中，可得到边坡稳定性分析的精确上下限解。Chen[18]最早将极限分析上限法应用于土体非均质性和各向异性下边坡稳定性研究。之后，众多学者考虑土体强度非均质性和各向异性影响，将极限分析上限法推广至不同工况下边坡的稳定性分析中[21-24]。Terzaghi[25]最早提出拟静力法用于边坡地震稳定性分析，本质上是由静力稳定分析方法扩展而来，机理是将地震的瞬时作用等效为作用于潜在滑体上的水平和竖直方向的加速度作用。由于使用简便，因而积累了大量的工程实际经验，便于实际应用推广。Chen等[26]在平动破坏机制和圆弧转动破坏机制基础上，提出一种更符合实际的对数螺旋破坏机制，得到地震力作用下非均质及各向异性边坡稳定性临界高度最优上限解。

已有地震力作用下考虑岩土体非均质及各向异性影响下的边坡稳定性研究多关注于加筋土边坡稳定性问题[27-28]，抗滑桩加固边坡稳定性问题[29]，预应力锚索加固边坡稳定性问题[30]，孔隙水压力作用下边坡稳定性问题[31]，而考虑土体强度非均质及各向异性影响下隧道洞口段仰坡地震稳定性研究相对较少。基于此，将土体黏聚力强度非均质和各向异性与极限分析上限理论相结合，同时引入土体强度折减理论，推导出水平和竖向地震力作用下隧道洞口段仰坡稳定性极限分析上限解计算公式，对比验证所提计算方法应用的合理性，得到地震力作用下非均质及各向异性系数对隧道洞口仰坡临界坡高、隧道拱顶破坏位置及仰坡安全系数的影响规律。以期为山区隧道洞口仰坡抗震稳定性设计提供理论参考。

1.1 计算模型

计算模型如图1所示，依据相关联流动法则，隧道洞口段仰坡滑体边界为CD，滑体区域ABCD以角速度ω绕旋转中心O作整体转动破坏。隧道洞口段仰坡坡高为Hcr，坡面AB水平向倾角为β，仰坡滑移线起始点位置与坡肩距离为L，rh为潜在滑移面剪出口A′到图1中滑体旋转中心O的半径，仰坡滑移线穿过隧道顶部位置与隧道洞口距离为Lt，α1为隧道洞口段仰坡高度相关系数、α2为隧道拱顶与潜在滑移面剪出口之间距离的相关系数。假定隧道洞口段仰坡为理想弹塑性模型，遵循Mohr-Coulomb强度准则，基于Chen等[18,26]提出的边坡极限状态下的对数螺旋转动破坏机制，结合图1可得仰坡潜在滑移面方程，可表示为

r=r0exp[(θ-θ0)tanφ]

(1)

式(1)中：r为对数螺旋线上任一点到图1中滑体旋转中心O的半径；
r0为对数螺旋线起点C到图1中滑体旋转中心O的半径；
θ、θ0分别为对数螺旋线上任一点与滑体旋转中心O的连线与水平参考线的角度及对数螺旋线起点C与滑体旋转中心O的连线水平参考线的角度；
φ为土体内摩擦角。

由图1可得到对数螺旋转动破坏机制下的滑动面高度H、隧道洞口段仰坡高度Hcr、潜在滑移面分布L及隧道拱顶破坏位置Lt几何关系无量纲表达式为

Hr=H/r0=exp[(θh-θ0)ψm]sinθh-sinθ0

(2)

H1r=Hcr/r0=sinθdexp[(θd-θ0)ψm]-sinθ0

(3)

c为土体黏聚力；
θd为潜在滑移面在隧道顶部剪出口D点与滑体旋转中心O的连线与水平参考线的角度；
n0为土体非均质系数；
Kh为水平向地震力系数；
Kv为竖直向地震动系数；
W为土重；
σ1、σ3分别为大主应力和小主应力；
V(θ)为潜在滑移体运动速度；
黏聚力ci与最大主应力方向一致，与竖直方向呈i角；
δ为破坏面的切线与垂直于最大主应力线的夹角；
ch和cv分别为土体水平向和竖直向黏聚力图1 计算模型Fig.1 Model for calculation

Lr=L/r0

=sin(θ0+β)/sinβ-

sin(θh+β)exp[(θh-θ0)ψm]/sinβ

(4)

Ltr=Lt/r0

=cosθdexp[(θd-θ0)ψm]-

cosθ0+Lr+Hcrcotβ

(5)

式中：依据强度折减原理，ψm= tanφm= tanφ/Fs，其中Fs为土体强度折减系数；
若稳定性计算中此折减系数增加值边坡发生破坏时，此时的折减系数可称为边坡的安全系数，同样可得到cm=c/Fs。

1.2 土体强度非均质及各向异性

考虑隧道洞口段仰坡土体黏聚力强度指标表现非均质和各向异性。土体强度各向异性表示土体的黏聚力随方向发生改变的特性，图1中黏聚力ci与最大主应力方向一致，与竖直方向呈i角，δ为破坏面的切线与垂直于最大主应力线的夹角，角度δ的值为π/4+φ/2[17,32]。进一步参考Lo[33]的研究成果，各向异性土体黏聚力计算表达式为

ci=ch+(cv-ch)cos2i

(6)

式(6)中：k=ch/cv为土体强度的各向异性系数。

土体强度非均质性表示土体黏聚力随着深度的增加而线性递增的特性，如图1所示，假定土体非均匀性沿深度呈线性变化。当n0= 1.0时，土体强度表现为各向均质特性。土体同时受各向异性和非均匀性影响的黏聚力可表示为

[(1-n0)/Hr][sinθe(θ-θ0)tanφ-sinθ0]}

(7)

1.3 稳定性上限解

极限分析上限法认为对任一满足运动许可的破坏机构，当令其外力功所做功率等于内能耗散率时，可得到极限破坏荷载的一个不安全上限，即所确定的极限状态下的破坏荷载不小于真实极限荷载[18]，可表示为

(8)

式(8)中：σij为通过相关联流动法则与运动许可速度场中的应变率εij相关联的应力场；
vi为满足几何相容条件的速度场；
Ti、Fi分别为机构边界S上的面力和体积V中的体力。

基于极限分析上限理论，假设隧道洞口段仰坡失稳后以刚体形式运动，不考虑坡体内的局部破坏，当外力对潜在滑移体所作的功率大于坡体内部耗散功率时，仰坡将发生整体破坏。因此，在地震条件下，作用于潜在滑移体上的外力主要由滑体重力和地震力共同组成，坡体内部耗散功率主要是沿潜在滑移面的土体内能耗散功率。

1.3.1 潜在滑移体重力做功功率

滑动土体AEBCDA区域的外力功率为Wr，采用将A′EFA′区域土体重力做功功率减去A′AD区域土体重力做功功率及BFC区域土体重力做功功率。因此，AEBCDA区域的外力功率可表示为

(9)

式(9)中：γ为土体重度；
f1～f4为与θ0、θd、θh、β、φ相关的函数，可分别表示为

(10)

1.3.2 地震力做功功率

隧道仰坡同时受水平向及竖直向地震力作用，依据极限分析上限理论，地震力对坡体所做功率为

(11)

式(11)中：Kh为水平向地震力系数；
Kv为竖直向地震动系数；
Wev、Weh分别为竖直向地震力做功功率、水平向地震力做功功率；
f5～f8为与θ0、θd、θh、β、φ相关的函数，可分别表示为

(12)

1.3.3 潜在滑移面内能耗散率

潜在滑移面内能耗散率计算由图1中对数螺旋间断面CD的微分面积rdθ/cosφm与土体黏聚力ci及该间断面上切向速度vcosφm相乘，总的内能耗散率为沿间断面CD积分后计算可得

(13)

式(13)中：φm= arctan(tanφ/F)。

[(1-n0)/Hr][sinθe(θ-θ0)ψm-sinθ0]}

(14)

fd=fd1+fd2

(15)

式(15)中：fd和fd2可分别表示为

(16)

(17)

式(20)中：fd11、fd12及fd13可分别表示为

(18)

(19)

(20)

依据极限分析上限法原理，地震力作用下，考虑土体强度非均质及各向异性的隧道洞口段仰坡外力所做功率和内能耗散率应满足以下平衡方程：

(21)

由式(21)可推导得到隧道洞口仰坡临界高度Ns的表达式为

(22)

式(22)中：由临界高度所确定的仰坡高度Hcr是仰坡在临界状态下(安全系数Fs= 1.0)所能达到的最大高度。

进一步变换式(22)，得到可用于求解仰坡安全系数Fs的隐函数表达式为

=[(1+λKh)(f1-f2-f3-f4)+

Kh(f5-f6-f7-f8)]/(KfdH1rψ)-1

(23)

对于滑面不确定的隧道洞口段仰坡，坡体内存在众多潜在滑移面，其中得到最小仰坡临界高度所对应的滑移面为最危险滑移面。仰坡临界高度Ns及稳定性系数Nf可表示为参数θ0、θh、及θd的函数。基于MATLAB数值优化方法，便可确定相应于某一仰坡临界高度的潜在最危险滑移面位置及隧道拱顶破坏位置，为使隧道洞口仰坡破坏机制满足几何条件，仰坡临界高度最小上限解需满足约束条件：

(24)

何毅等[34]利用极限分析上限法求解得到各向同性及均质条件下边坡临界坡高随坡角变化规律，而土体各向异性系数n0= 1.0，非均质性系数k= 1.0可以认为是土体强度各向异性及非均质条件的特例。将本文方法退化为各向同性及均质条件与其稳定性解答对比，如图2所示。可以看出，边坡不同坡角β及内摩擦角φ影响下得到的边坡临界坡高Ns=γH/c解与文献[34]结果非常接近。但当β= 40°、φ= 20°时，本文方法得到的临界坡高略小于其解答，这是由于两者优化方法存在差异。

为进一步验证所提方法在地震力作用下的适用性，进一步将所得到的解答与Ling等[35]基于极限平衡法得到的不同地震力系数λ、Kh影响下的边坡稳定性系数解c/γH进行对比，如表1所示，发现

图2 静力条件下边坡临界坡高计算结果对比Fig.2 Comparisons of critical height influenced by static conditions

表1 地震力作用下边坡稳定性系数计算结果对比Table 1 Comparisons of stability numbers influenced by seismic force

应用本文方法得到的解答与文献[35]解答较吻合，且采用极限分析上限解普遍小于极限平衡解，说明极限分析上限解较极限平衡解更接近真解。因此，上述两则算例的对比结果验证了所提方法用于研究隧道洞口段仰坡稳定性的可行性。

3.1 临界坡高

稳定图可用于判定已有坡体高度是否满足边坡稳定性要求。针对隧道洞口段仰坡，通过式(22)可得到边坡临界坡高γHcr/c，若岩土体重度γ及黏聚力ci已知，应用绘制的稳定图即可得到保证隧道洞口段坡体稳定的临界坡高。图3、图4给出土体各向异性、非均质性及地震力作用下的隧道洞口段仰坡临界坡高变化规律。可以看出，非均质系数n0越大，γHcr/c越大，表明隧道洞口段仰坡的临界坡高Hcr越大。而各向异性系数k越大，γHcr/c越小，表明维持隧道洞口段仰坡稳定所需的临界坡高Hcr也越小。地震力对隧道洞口段仰坡临界坡高影响明显，表现为随着水平地震力系数Kh的增大，γHcr/c减小，说明地震力越强时，仰坡为维持稳定所需的临界坡高越小。此外，仰坡临界坡高与隧道洞口段仰坡坡角β密切相关，如当坡角较小时，不同水平地震力系数Kh影响下的γHcr/c值差异显著，而随着坡角增大，这种差异逐渐减小。图4中，λ为正值表明地震力作用与重力一致，即震动方向竖直向下，而λ为负值表明地震力作用与重力相反，即震动方向竖直向上，可以看出地震力竖直向上作用时的隧道洞口段仰坡临界坡高明显大于竖直向下的地震力作用时的临界坡高。

图3 坡角及地震力对隧道洞口段仰坡临界坡高影响规律Fig.3 Effect of slope angle and seismic force on critical height of slopes at tunnel entrance

图4 土体非均质及各向异性对临界坡高影响规律Fig.4 Effect of soil inhomogeneity and anisotropy on critical height of slopes at tunnel entrance

3.2 隧道拱顶破坏位置

依据图1建立的隧道洞口段仰坡破坏机制，潜在滑移面穿过隧道顶部，往往造成隧道拱顶发生破坏，因此，隧道拱顶破坏位置预测是进行隧道加固的重要环节。图5～图8给出土体各向异性、非均质性及地震力作用与潜在滑移面的初始位置L/Hcr及隧道拱顶破坏位置Lt/Hcr的关系曲线。可以看出，当土体各向异性系数n0及非均质性系数k增大时，潜在滑移面的初始位置L/Hcr及隧道拱顶破坏位置Lt/Hcr都呈增大的趋势，但随着k的增大，不同n0影响下的L/Hcr及Lt/Hcr差距逐渐缩小，直至k= 1.0时，数值基本趋于一致。同时，还可以明显发现隧道洞口段仰坡潜在滑移面初始位置分布曲线与隧道拱顶破坏位置曲线规律保持一致，说明潜在滑移面位置与隧道拱顶失稳位置密切相关，这是由于坡顶滑移位置及隧道拱顶破坏位置分别对应潜在滑移面的起始及剪出位置。

图5 土体各向异性及非均质性对潜在滑移面 初始位置L/Hcr影响规律Fig.5 Effect of soil inhomogeneity and anisotropy on initial position of potential sliding surface L/Hcr

图6 土体各向异性及非均质性对隧道拱顶破坏 位置Lt/Hcr影响规律Fig.6 Effect of soil inhomogeneity and anisotropy on failure location of tunnel vault Lt/Hcr

图7 坡角及地震力对潜在滑移面初始 位置L/Hcr影响规律Fig.7 Effect of slope angle and seismic force on initial position of potential sliding surface L/Hcr

图8 坡角及地震力对隧道拱顶破坏位置Lt/Hcr影响规律Fig.8 Effect of slope angle and seismic force on Lt/Hcr

从图7、图8中可以观察到，水平向地震力系数Kh增大且竖直向地震力作用方向向上时，L/Hcr及Lt/Hcr皆增大，说明地震力作用将导致隧道洞口仰坡发生深层滑动，即初始滑移面远离坡肩位置，隧道拱顶破坏范围更广。坡角β对L/Hcr及Lt/Hcr影响较大，表现为随着坡角β的增大，坡顶初始滑移位置呈先增大后减小的趋势，在β= 60° ～ 70°时，L/Hcr达到最大值，即坡顶初始滑移位置离坡肩最远。当坡角小于40°时，Lt/Hcr呈递增趋势，超过此角度，随着坡角的增大，Lt/Hcr呈递增趋势，即仰坡潜在滑移面穿过隧道拱顶的位置距离隧道洞口较近。

3.3 仰坡安全系数

安全系数作为评价边坡稳定性的重要指标，已在工程实践中得到普遍应用。然而，极限分析上限法无法直接求得边坡安全系数，基于此，结合强度折减理论及含安全系数Fs的隐函数方程式[式(15)]，通过绘制一种横坐标tanφ/Fs、纵坐标c/γHcrtanφ的稳定图间接计算隧道洞口段仰坡安全系数。计算实现过程如下：算例中设定隧道洞口段仰坡坡角β= 60°，各向异性系数n0= 0、0.2、0.4、0.6、0.8、1.0，非均质系数k= 0.5、0.6、0.7、0.8、0.9、1.0，λ= 0.5，Kh= 0.1。

图9 给出各向异性系数n0、非均质系数k影响下隧道洞口仰坡稳定图。可以看出，稳定图呈显著的非线性分布特征，非均质系数k影响下的稳定图曲线差异较小，各向异性系数n0的变化对稳定图曲线分布影响较明显。给定仰坡坡高Hcr=10.0 m，土体参数c= 40 kPa，φ= 20°，γ= 18 kN/m2，当k= 1.0,n0= 1.0;k= 0.6,n0= 1.0;k= 1.0,n0= 0.8时，可得到纵坐标c/γHcrtanφ= 0.611，进一步找到图9中该纵坐标值对应的横坐标tanφ/Fs值，分别为tanφ/Fs= 0.357、0.308、0.379，将φ= 20°代入可反算得到安全系数分别为Fs= 1.02、1.18、0.96。

针对土体强度各向异性及非均质性的不同分布情况，研究n0及k值变化对隧道洞口仰坡安全系数的影响规律，如图10所示。可以看出，在地震力作用下，随着各向异性系数n0的增大，仰坡安全系数Fs呈递增趋势，而非均质系数k越大，仰坡安全系数Fs越小。

图9 土体非均质及各向异性影响下仰坡稳定图Fig.9 Stability charts of slopes at tunnel entrance considering soil inhomogeneity and anisotropy

图10 隧道洞口仰坡安全系数与土体非均质及 各向异性关系Fig.10 Relationship between safety factors of slopes at tunnel entrance and soil inhomogeneity and anisotropy

应用极限分析上限法结合土体强度折减理论，构建了考虑土体强度非均质及各向异性的隧道洞口仰坡地震稳定性分析模型，推导出稳定性极限分析计算公式，得到仰坡临界坡高、隧道拱顶破坏位置的最优上限解及相应的稳定图，并建立一种以稳定图计算仰坡安全系数的方法。得出如下主要结论。

(1) 非均质系数n0越大，隧道洞口段仰坡的临界坡高Hcr越大，而各向异性系数k越大，维持隧道洞口段仰坡稳定所需的临界坡高Hcr也越小；
坡角较小时，不同水平地震力系数Kh影响下的临界坡高差异显著，而随着坡角增大，差异逐渐减小。

(2) 土体强度各向异性系数n0及非均质性系数k增大时，潜在滑移面的初始位置L/Hcr及隧道拱顶破坏位置Lt/Hcr皆增大，表明初始滑移面远离坡肩位置，隧道拱顶破坏范围更广。

(3) 通过引入算例绘制稳定图提出间接计算隧道洞口段仰坡安全系数的方法，发现在地震力作用下，随着各向异性系数n0的增大，仰坡安全系数Fs呈递增规律，而非均质系数k越大，仰坡安全系数Fs越小。