基于上限定理的含裂缝边坡稳定性图表分析方法

高明辉， 张雅贤， 王奇智， 赵建青， 孙 超， 周记名， 张娜娜

(1.河北省地质资源环境监测与保护重点实验室， 石家庄 050021； 2.中央军委后勤保障部安置住房保障中心， 北京 100036； 3.河北科技大学建筑工程学院， 石家庄 050018； 4.中国土木工程集团有限公司， 北京 100038)

大量道路边坡现场勘测发现，边坡表面分布着大量裂缝，而裂缝的存在对其整体安全性的影响显著，因此有众多学者对裂缝边坡展开了深入的研究[1-5]。并对裂缝的扩裂范围进行了限定，但常常忽略了裂缝位置存在的影响[6]，导致得到的临界裂缝破坏状态不是最不利裂缝诱发的。因此，近年来，不断有学者对裂缝扩展的形态和出现的区域进行理论上的创新和方法的改进。Utili[7]使用Chen[8]提出的极限分析上限法的同时，考虑了裂缝可能存在的不同深度和位置；Michalowski[9]则将裂缝的扩展与边坡的破坏机制相结合进行研究分析，而后以此研究为基础，进行了深入的探讨与分析。

导致裂缝边坡失稳和滑坡的因素很多，典型的如降雨[10]、地震[11]、交通负荷等或这些因素的共同作用。在边坡危险性评估分析时，稳定性系数Ns、安全系数Fs及动力荷载下的边坡位移大小，均可作为评估其是否失稳的指标，而安全系数Fs一直是预测坡体失稳的常用方法。强度折减法是一种简单、有效、系统的确定安全系数取值的技术，学者们在岩土工程设计计算中广泛采用了这一技术[12-14]。当使用强度折减法评估边坡稳定性，需要假想坡体两个重要的参数：黏聚力c和内摩擦角φ，同时折减到相同的安全系数。

对于常见的均质边坡设计，Baker等提出了基于极限平衡法的静力和动力条件下边坡稳定性分析的完整设计图；Michalowski等[15-16]提供了一系列基于极限分析方法考虑孔隙水压力和地震力的设计图表；Steward等[17]提出了计算5种不同破坏类型下土质边坡安全系数的两种设计图表。前人的研究表明，各类方法，如极限平衡法(limit equilibrium method)、有限元法(finite element method)或有限差分法(finite difference method)、极限分析法(limit analysis method)等，已被广泛使用以准确得到边坡的稳定性参数。以往研究中，较为常见的是传统的极限平衡法[18]，特别是生成需要大量安全系数数据设计图表时，与有限元方法相比，它的时间消耗更少。然而，传统的极限平衡法并未考虑到土体的本构关系。而极限分析上限法考虑了能量的守恒，有效解决了土体应力应变关系的问题，并且给出了精确的临界上限值，使得当前极限分析上限法成为分析动静态边坡稳定时经常采用的方法。

现基于经典的极限分析上限法，结合强度折减法对裂缝边坡的安全系数进行综合性分析，针对常见地震荷载条件下、裂缝中的静水压力以及坡顶部荷载作用等进行研究分析，将临界裂缝的形态及其对应的深度、位置展示，并与有限元软件Optum-G2的结果进行对比，以期为边坡安全性初步快速评估提供科学依据。

1 基本破坏机制以及公式推导

R0为摩尔圆半径；δ为剪胀角； f′c为单轴抗压强度；σ为应力； τ为剪切力；速度矢量δw； (δu, δv)为速度分矢量图1 边坡破坏强度包络线Fig.1 Slope failure strength envelope

理想的塑性边坡往往表现为达到极限应力时开始发生屈服。经典Mohr-Coulomb(M-C)屈服准则往往用于评估边坡安全性时首选的破坏准则，即屈服准则中采用线性的强度包络线。采用M-C屈服方法时，抗剪强度由土体的黏聚力c和内摩擦角φ来预测，如图1所示。而对土体在实测的压缩状态往拉伸方向外推得到其抗拉强度，所以并不是实际实验而得到的抗拉强度。而针对裂缝边坡，裂缝的存在可以有效地解决边坡抗拉强度的问题，因此用经典的M-C屈服准则对分析是合理的。而使用强度折减法进行边坡稳定性分析时候，需要同时将土体的黏聚力c和内摩擦角φ折减到一定的程度Fs,可表示为

(1)

式(1)中：cf、φf分别为折减后的黏聚力和内摩擦角。

极限分析上限法基于能量的守恒，破坏机制可以包含连续变形的区域，也可以由刚性块构成：平移(或转动)，由中间脆弱的剪切带隔开，通常被视为运动场中的速度不连续性。当用线性屈服条件描述土体的强度并用相关联流法则约束变形时，速度矢量δw与强度屈服面呈垂直状态，可将其分解成两部分：水平分量δu与竖直分量δv。通过应力矢量(σ,τ)和速度矢量(δu,δv)的乘积来表示非连续面上每单位面积的能量耗散率dc，可表示为

dc=cvcosφ

(2)

式(2)中：v为单位旋转速度。

在M-C屈服条件下，将边坡外部做功与内能能量耗散(滑裂面做功)相平衡，如式(3)所示，破坏机制如图2所示。边坡右侧BE为竖直裂缝，即预设的边坡后缘的产生拉应力区的裂缝位置，裂缝具有一定的位置x(与破址角的直线距离)及一定的深度y,而裂缝在不同深度、位置的限制下产生不同的安全系数，从而形成最小的安全系数，x、y可由式(5)、式(6)得到。裂缝形态如图3所示，计算目标为可能出现最不利影响时裂缝所处的位置x和所达到的深度y。

β为坡脚；θ0、θc、θh为旋转角；khG为水平加速度； kvG为竖直加速度；G为重力，q为上部荷载； r0、rh分别为旋转点到A点和D点的距离；H为边坡高度图2 裂缝边坡破坏机制Fig.2 Failure mechanism of slope with crack

L2为裂缝处到坡上缘的距离图3 临界裂缝位置及深度Fig.3 Critical crack position and depth

(3)

在实际破坏条件下，地震加速度并没有固定传递方向。而计算地震作用引起的外部功率的推导方法与边坡自身重量引起的外部功率的方法一致。在计算水平外力作用下的外部功的比率时，唯一的区别是速度方向。因此，将极限分析法与强度折减法相结合，可以得到裂缝边坡在地震条件下的安全系数计算公式为

[(f1-f2-f3-t1+t2+t3)(1+λkhG)+

(4)

(5)

sinθ0}-1+cotβ

(6)

考虑裂缝中静水压力εw影响时，引入无量纲系数K表示裂缝中的含水量，K=0时表示裂缝中未出现静水压力，而K=1时表示裂缝中充满了静水压力。则加入裂缝中的静水压力后的安全系数Fs、fr的计算公式分别为

(7)

(8)

式(8)中：γw为水的容重。

考虑边坡上部荷载作用时，引入无量纲系数qt，qt=q/γH。整理后的边坡安全系数Fs和fq的表达式分别为

[(f1-f2-f3-t1+t2+t3+qtfq)+

(9)

(10)

(11)

2 程序优化与结果分析

先前求解安全系数的计算方法存在显著的不足，常常存在搜索精度不高、迭代计算复杂以及无法获取全局最优解的问题。文章为解决这类问题，采用随机搜索法与序列二次优化法(sequential quadratic programming,SQP)结合的计算方法，即采用MONTE-CARLO算法在一定精度范围内确定一个初值，再使用SQP方法确定全局最优解，计算效率显著提高，且精准度更高。具体计算逻辑如图4所示。

图4 优化算法逻辑图Fig.4 Optimization algorithm logic diagram

考虑破坏工况为：坡脚30°～80°、内摩擦角10°～30°、参数组合c/γH为0.05～0.3、地震系数(0～0.2)g(g为重力加速度)、裂缝静水压力系数：0～0.3及上部荷载系数：0.1～0.15。为验证计算的准确性，将结果与经典结果[19]进行对比，对比参数如表1所示，计算结果如表2所示，结果显示，针对裂缝边坡，可以明显看到其相对于完整边坡(无裂缝)呈现出更加不利的安全系数，因此对于裂缝边坡的破坏稳定性研究，尤为关键。

表1 算例参数设计

表2 不同计算方法结果对比

2.1 考虑地震作用的影响

针对地震的影响，采用拟静力法进行研究分析，即针对khG=(0～0.2)g进行研究分析。横坐标表示边坡坡脚β，纵坐标表示计算求得的安全系数Fs，考虑内摩擦角的变化，以边坡参数组合c/γH为变量进行研究分析，研究地震力作用时，为仅考虑地震荷载影响的结果，因此并未加入裂缝中的静水压力以及坡顶上部荷载的影响。结果显示，地震荷载对于边坡安全系数的影响非常显著，图5(a)～图5(c)展示了不同地震荷载下的裂缝边坡的安全系数，结果显示考虑地震荷载的作用大大折减了Fs，如针对c/γH=0.05、β=30°、φ=20°的边坡，不考虑地震荷载时，其安Fs=1.24，坡体整体处于安全状态，而考虑地震影响时候，如khG=0.2g时，裂缝边坡的Fs=0.85，相比较静态降低了31%，且此时坡体处于滑落状态，因此地震作用力对于边坡安全性的影响尤为显著。而竖向地震动对于边坡稳定性也存在一定程度的影响，图5(d)展示了不同程度的竖向分析系数λ对于边坡安全系数的影响，结果显示，考虑竖向地震荷载的存在对于边坡的影响显著，考虑竖向地震系数λ为正时，即向下的竖向加速度，边坡的安全系数随之降低。在针对考虑地震力时的带裂缝边坡安全分析时，工程师可直接采用图5中一系列图表参数进行选取评估。

2.2 考虑边坡上部荷载作用的影响

将上部荷载的作用加入裂缝边坡的讨论中，考虑上部荷载系数qt(qt=q/γH)为0.05～0.15的工况，计算结果整理到图6中，与图5坐标轴一致，区别了地震荷载与上部荷载的作用，结果表明：上部荷载作用下，边坡的安全系数产生一定的折减，尤其针对c/γH较大时候的工况，而边坡的Fs同样随着上部荷载qt的增加而降低。上部荷载的作用将会对边坡的稳定性带来较大程度的影响，如针对c/γH，β=80°、φ=10°的边坡，当不考虑上部荷载作用时，根据图6(a)可见，边坡Fs=1.18，而考虑上部荷载的作用，当qt=0.15时，边坡Fs=0.95，降低了19.5%，折减程度明显，因此边坡稳定性分析时，需将上部荷载的作用进行讨论。

图5 地震荷载作用下裂缝边坡的安全系数Fig.5 Safety factor of slope with crack under earthquake load

图6 上部荷载作用下裂缝边坡的安全系数Fig.6 Safety factor of slope with crack under surcharge load

2.3 裂缝中静水压力的影响

将裂缝中的静水压力的作用加入裂缝边坡的讨论中，考虑无量纲的静水压力系数K为0.2～0.6的条件，结果如图7所示，在坡度较大的裂缝中，同样产生较大的折减，同样针对c/γH=3、β=80°、φ=10°的边坡，在考虑空隙水压力K=0.3时，边坡Fs=1.02，折减程度为14%，同样显著。因此对裂缝边坡设计时，不可忽略裂缝存在的静水压力的作用。

图7 静水压力作用下裂缝边坡的安全系数Fig.7 Safety factor of slope with crack under hydrostatic pressure

3 裂缝位置搜索及图表应用

3.1 裂缝位置的搜索

基于表1的工况，进行了3种条件下的边坡破坏临界裂缝的位置与深度，计算结果如表3所示。可以明显看出，更低的安全系数将导致距离破址更远的裂缝位置x/H与距离坡顶面更大的裂缝深度y/H。为更加形象展示计算结果，将滑裂面破坏形态进行展示，即分别考虑khG=0.2g，上部荷载系数qt=0.26以及裂缝静水压力系数K=0.2与仅考虑自身重力破坏的边坡临界滑裂面进行对比分析，并将其裂缝的位置、以及裂缝深度进行对比分析。结果如图8所示。可以明显看出，仅考虑自重条件下，裂缝的位置更加靠近坡顶面，而考虑不同条件的破坏工况时，裂缝出现在较远处。为验证裂缝位置、深度的准确性，采用有限元软件Optum G2计算结果进行了对比，对比两类工况：即khG=0.2g、qt=0.26的两类，从图9可知，计算得到的临界滑裂面与软件结果几乎一致，由此再次验证本文公式推导、计算结果均具备合理性。

图8 典型条件下临界状态滑裂面形态Fig.8 The morphology of sliding crack surface in critical state under typical conditions

图9 本文结果与Optum G2的破坏面对比Fig.9 Comparison between the current analytical results and Optum G2 results

3.2 图表选取运用

图5～图7可作为工程设计时安全系数的初步评估。为反映图表的有效性，针对一个案例进行示范。如假设实测β=45°、φ=20°、c=20 kPa、γ=17 kN/m3、H=12 m的边坡，在不考虑任何工况条件时，根据图5(b)显示，Fs=1.27时边坡处于安全状态，而针对khG=0.2g时可读取Fs=0.95，处于不安全的状态，因此设计时需要进行更多考虑。

4 结论

边坡后缘往往分布着裂缝，当前对于裂缝边坡的研究尚且缺乏，为准确分析裂缝边坡安全系数，将LAM、强度折减法、拟静力法综合研究推导。计算时，用蒙特卡洛随机搜索法与序列二次优化法结合进行数值计算。得出如下结论。

(1)采用蒙特卡洛与SQP算法相结合可以有效解决先前算法的缺陷，且结果精度更加准确，计算速度更快。

(2)三类典型破坏条件：地震、静水压力、上部荷载的存在均会对边坡的安全系数产生一定程度的影响，而对于坡度较高、内摩擦角较小的边坡影响更为显著。

(3)对临界破坏状态下裂缝的位置、深度还原，结果表明，相同参数条件下，产生较小的安全系数，而较小的安全系数易产生更深的临界裂缝深度以及更远的临界裂缝距离。

(4)针对不同条件的计算结果与传统文献经典解进行了对比分析，其精度优于传统分析结果，且计算结果与有限元结果基本吻合，因此绘制了大量安全系数的图表，可作为快速评估裂缝边坡安全系数的参考。