考虑坡体后缘张拉应力折减的超载边坡稳定性分析

王汉敬

(中国土木工程集团有限公司 北京 100038)

1 引言

地震是自然界最具破坏性的地质灾害之一，不断有学者研究地震对土工构筑物的影响[1－4]。其中，对于当前研究边坡稳定性，通常采用线性的Mohr-Coulomb屈服准则(简称M-C准则)来评估其安全性，鲜有关于岩土体非线性强度特征及其对边坡稳定性影响的相关研究。但土工试验结果表明:采用线性的强度包络线，即将土体的黏聚力和内摩擦角用以计算抗剪强度，进而外推得出土体的抗拉强度，但此强度并非是土体真实的抗拉强度[5]。另一方面，通过抗拉强度的研究发现[6]，采用两种方法可以有效地分析边坡拉应力区的抗拉强度，使分析结果更加趋于真实的破坏条件，其一是对岩土体的抗拉强度进行一定程度的折减，即将线性破坏曲线优化为非线性强度包络线；其二则是在边坡后缘的拉应力区引入不同形态的裂缝。对于后者考虑裂缝的存在，Michalowski等[7－8]将裂缝的扩展过程与整体的滑落状态结合，研究结果表明:对于高陡边坡或受静水压力影响显著的边坡，裂缝的存在大大增加了失稳的风险。Zhao等[9]分析了线性准则下不同大小的抗剪强度参数对带裂缝边坡的地震稳定性的影响，研究显示地震荷载对边坡的稳定性影响尤为显著。Baker等[10－11]对拉剪区域进行了一定的折减，并结合变分法对强度包络线中的张拉区域进行相关的稳定性分析。康孝森等[12]基于非饱和土力学和最大伸长线应变理论推导出入渗影响下黄土边坡坡顶张拉裂缝极限深度控制方程。

实际上，超载边坡是广泛存在的，且目前关于超载边坡尚未有考虑土体拉应力区强度折减后的超载边坡动态全过程分析的相关研究。为研究抗拉强度折减对超载边坡在地震作用下稳定性的影响，本文在一定研究成果的基础上，首先通过极限分析上限法对超载边坡临界的破坏形态以及功能方程进行推导；然后对计算求解算法进行深度优化，结合随机搜索法(Monte Carlo Method)与二次序列规划程序(SQP)进行求解；最后针对真实地震中边坡破坏案例，选取实测地震动时程，分析不同超载系数以及张拉强度折减系数对边坡稳定性的影响。

2 土体强度包络线

对边坡稳定性进行分析时，土体一般采用理想均质材料，而M-C准则并不能精确反映土体的抗拉强度。大量的试验结果表明，非线性的强度破坏包络线能更真实地反映岩土体的抗拉强度。

经典的M-C准则在主应力空间可以表示为一个不规则六边形截面的锥体，如图1a所示。后来学者对经典的M-C屈服面进行多次修改，从而限制了抗拉强度，其屈服面如图1b所示。折减后的部分拉力截断单轴抗拉强度ft可表示为图2a中的应力圆C1。由图2中的莫尔圆转换至图1中的主应力空间是一个较为复杂的过程，图1b中截面ABCDEF上AB边上的应力状态由图2b中的应力圆C2表示，而AB上各个位置由其中主应力的大小确定(σ1≥σ2≥σ3)。例如，图2中的应力圆C1在σ2＝0的情况下对应于图1b中ABCDEF截面与σ3轴的交点。极限状态时的屈服面由图2中光滑的屈服面映射至图1中具有棱角的屈服面。采用此拉应力区的截断方式适用于极限分析中的解析运动学分析方法，但因为主应力空间中屈服面棱角处的法向量不唯一，其在数值方法中的适用性还需深入研究。传统的M-C强度包络线如图2a所示，M-C函数为简单的直线与圆弧部分的组合。在M-C计算中，三轴抗拉强度的大小为f3t＝c/tanφ，而单轴抗拉强度大小为ft＝2ccosφ/(1＋sinφ)。研究中当考虑抗拉强度的折减时，牵引矢量从直线变化为圆弧，且张拉应力区域的剪胀角δ从φ变化到π/2，当δ＝φ时代表线性的M-C屈服准则。为了量化展示折减后的抗拉强度，本文引入了抗拉强度折减系数ξ:

其中，ξ＝0～1。当ξ＝1时为线性M-C准则下的抗拉强度；当ξ＝0时表示完全不考虑抗拉强度时的破坏状态，如图2b所示。采用此拉伸截断方式令极限分析的解析运动学分析过程大大简化。

图1 屈服准则模型

图2 不同包络线示意

3 破坏机制与计算分析

3.1 破坏机制分析

图3为受地震荷载作用的临界均质超载边坡模型。块体BDEC以中心O为中心沿着滑裂面DEC旋转，而整个破坏为过坡趾的破坏，DEC部分为临界滑面。研究发现，张拉应力更容易出现在坡顶附近，而滑面上部DE部分由受拉伸区折减部分强度包络线而引起，下方较低的部分(EC)是由线性部分的包络线引起。DE部分剪胀角δ从D点的最大值δmax向E点的φ不断折减。图4则展示了DE面上的一个微元，可得:

沿DE积分可以得到:

图3 受地震荷载作用的临界均质超载边坡模型

图4 边坡上滑面计算单元

图3～图4中，r0和θ0分别为对数螺旋线起始点A到旋转中心半径的距离AO和AO与水平参考线的夹角(见图2)；θh为对数螺旋线终止点到旋转中心的半径与水平参考线的夹角；θm为OD与水平参考线的夹角；θtc为OE与水平参考线的夹角；δ为剪胀角，其详细分布形式可表示为:

极限分析上限法基于块体内外能量耗散率的平衡，功能方程为:

式中，Wr为土体自重状态下的能量耗散率；Ws为拟静力地震荷载作用下的能量耗散率；Wqt为超载作用的功率；D为滑体BDEC的内能耗散率。外部功率Wr和Wr为块体ABC的功率减去块体ADE的功率，分别表示为:

公式中带有下标的参数分别表示图3中区域OCA、AOB、BOC、OEA、AOD、DOE 的能量耗散功率；fw与fs均为包含变量β、θ0、θtc、θm、θh的函数方程。分析超载部分功率时，引入超载系数qt，fq为超载函数。根据主应变率，由张拉区与折减产生的内部能量耗散为，则整个滑裂面的做功可表示为:

式中，fd为包含变量θtc、θh、φ的函数；fdc为包含变量θtc、θm、φ的函数。将上述公式整理，即可得到屈服加速度Kc的表达式:

Kc为一个含θ0、θtc、θm、θh以及c/γH、坡脚β和内摩擦角φ等参数的函数，则屈服加速度Kc的获取为一系列变量组合下的最小值。

3.2 计算优化

为搜索Kc在θ0、θtc、θm、θh、c/γH、β、φ等参数下的临界值，需对Kc进行求导分析，具体限定以及计算式为:0＜θ0＜θm＜θtc＜θh＜π

在限定的基础上，进行屈服加速度的求解。但采取传统的计算方式存在一些缺陷，如遍历法计算速度慢、随机搜索法计算精度低以及优化分析易于陷入局部最优解等。本文采用随机搜索法和非序列二次优化法(SQP)结合进行分析求解，随机搜索法精度虽然较低但是可以很好地给出初值的区间范围，而得到较为合理的初值后可采用序列二次优化进行计算临界值搜索。部分验证结果见表1～表2。其中表1为参数设计，表2为与经典解的对照。可以看出计算结果与经典解非常吻合。

表1 边坡参数设计

表2 本文计算结果与经典解对比结果

4 屈服加速度计算及地震破坏案例分析

4.1 屈服加速度计算

本文针对不同程度的超载边坡进行稳定性分析，考虑超载系数qt＝0 ～0.3、坡脚β＝30°～90°、内摩擦角φ＝30°、c/γH＝0.1～0.25，以及拉应力区折减系数ξ＝0～1的破坏工况。图3中上部点线部分表示过坡趾下方的破坏(Below-toe failure)；其余部分线条表示过坡趾的破坏(Toe failure)；图中原点为过渡点。由图3可知，过坡趾破坏为摩擦土体的特征，而当内摩擦角较小时，对于坡脚较为缓和的边坡，则更容易发生过坡趾下方的破坏。但考虑地震荷载影响时，对于内摩擦角更大的边坡易发生过坡趾的破坏。而由图5可知，考虑张拉区域折减分析，边坡的屈服加速度折减迅速，而针对内摩擦角较小且坡度较大的边坡，破坏更为强烈。超载作用则会加速屈服加速度的折减。图中ξ＝1表示线性M-C屈服分析结果。

图5 不同超载系数下的边坡屈服加速度

4.2 地震破坏案例分析

选取2018年5月12日的里氏8级汶川地震为研究对象，记录台站为512MZQ，加速度放大系数为0.868，采样间隔为0.005 s。为进行全面分析，本文选取破坏工况为c/γH＝0.15、β＝55°、φ＝30°。图6所示为超载边坡在不同程度张拉折减系数下的永久位移。结果显示，考虑张拉折减的影响，将对边坡永久位移增加2倍甚至2倍以上，如qt＝0.3，边坡在M-C准则与考虑ξ＝0时产生的位移分别为123 cm与195 cm，考虑强度折减后永久位移增加了59%。而图7显示，当非线性系数ξ取0.4～1.0时，非线性准则对于边坡永久位移影响并不明显；而当ξ逼近0时，边坡永久位移迅速增加，影响显著。而结果曲率显示，考虑超载作用对于边坡永久位移的折减程度显著大于非超载边坡。

图6 汶川地震作用下考虑张拉折减边坡永久时程位移分布

图7 汶川地震波作用下不同张拉折减系数边坡永久位移对比

5 结论

(1)采用Monte-Carlo与SQP综合优化法，显著提高运算效率，且本文计算结果与经典解进行对比验证相吻合。

(2)张拉应力区域折减系数ξ对于临界加速度Kc的影响显著，随着超载作用不断增加，影响在不断放大，而坡度较大、内摩擦角较小的脆性边坡更容易受到破坏。

(3)考虑张拉折减与超载共同作用时，边坡产生的永久位移将超过传统M-C准则下两倍，且当抗拉强度折减程度不断增加时，边坡位移积累增加迅速，而超载作用则会加速位移的上升。