基于极限分析上限法双层土坡稳定性分析

王 迪， 王宏权， 王晓飞， 徐惠民， 张 飞

(1.岩土力学与堤坝工程教育部重点实验室， 江苏 南京 210098；2.河海大学 土木与交通学院， 江苏 南京 210098； 3.广东省南粤交通大丰华高速公路管理中心， 广东 梅州 514300； 4.安徽省交通规划设计研究总院股份有限公司， 安徽 合肥 230088； 5.南京市水利规划设计院股份有限公司， 江苏 南京 210022)

边坡稳定问题是土力学三大经典问题之一，是岩土力学与工程研究中的重点领域。关于边坡稳定性问题，国内外学者提出了各种分析方法(极限平衡法、极限分析法、数值模拟法)，并被广泛应用于实际工程。边坡稳定性分析最早是基于Coulomb的挡土墙土压力计算方法和Rankine的主动土压力与被动土压力计算方法，建立了最初的极限平衡法。随后，条分理论的提出对边坡稳定性极限平衡分析法的发展有重要贡献，基于土体条分的各种假设建立了诸多方法[1-7]。由于土体塑性理论的建立，极限定理的思想被用于边坡稳定性分析中，提出了极限分析法。伴随着计算机技术和数值模拟方法的快速发展，各种数值分析方法被提出并引入到边坡稳定性分析中。但极限平衡法[8-9]由于将土条简化为刚体，忽略了岩土体的应力-应变关系，没有考虑土体发生破坏时的变形协调条件，计算得到的安全系数仅是假定滑动面上的平均安全度，结果并不十分理想。数值模拟法虽然具有简易的操作性和很强的可视化效果，但在本构模型、屈服准则以及失稳判据的选取方面缺少统一性，不同的模拟软件本身也对计算精度有着不同程度的影响。极限分析法的优势在于，它可以不必一步步地对各个阶段的弹塑性状态进行分析而直接求得破坏荷载。从理论上，极限分析法相比极限平衡法更为严格，极限分析法分为上限解法和下限解法，上限解需要建立运动许可的速度场(即破坏面)，下限解需要建立静力许可的应力场。若上限解与下限解相等时，则就认为获得了“真解”。

目前，边坡稳定极限分析法的研究主要针对均质土坡，但实际工程中的边坡问题常是非均质的，关于非均质边坡的稳定性研究仍需要进一步的改进与完善。Chen等[10]采用对数螺旋线的转动破坏机制，建立能量平衡方程得到均质边坡稳定性的上限解，随后该方法又被Chen等[11]拓展到非均质各向异性的边坡分析中。Chen[12]对分层位置在坡脚处的双层边坡进行分析，给出了限制深度下的稳定系数。然而，Chen[12]研究中只考虑了分层位置在坡脚处的情况，自然界中成层边坡的分层位置不一定都在坡脚处，对于任意分层位置双层土坡，目前尚未开展相关研究。此外，边坡失稳破坏时其破坏面的分布，一直是工程界比较关心的问题，但是针对不同参数变化对边坡滑动面分布规律的研究相对较少，尤其是对于常见的双层边坡。

针对以上研究的不足和需要，本文基于Chen[12]提出的双层边坡稳定性的极限分析上限解法，推导了任意分层位置的双层土坡在不同破坏机制下临界高度的计算公式，依托MATLAB软件平台，通过编制不同破坏机制下的稳定系数的计算程序，分析了临界高度与分层深度的关系，并得到不同条件下双层土坡临界高度随分层深度改变的线性变化规律。同时探究了双层土坡临界滑动面位置的分布规律，建立双层土坡的稳定性分析方法，为工程设计提供一定的参考依据。

1 任意分层位置双层土坡稳定性极限分析方法

极限分析法[13-14]是基于理想弹塑性体的基本定理-上限定理和下限定理求解物体极限应力状态的分析方法。对于岩土工程问题，通过一些简单的假定和边界条件的简化[15-16]，极限分析法提供了一种寻求合理真解范围的有效途径。

1.1 基于坡底破坏机制双层土坡临界高度计算公式

图1为双层土坡坡底破坏机制示意图，假设双层土坡破坏机制为旋转机构，破坏面为对数螺旋面，通过建立对数螺旋面上的能量耗损率与旋转土体重力所做外功率的等量关系，得到双层土坡极限分析上限解，进而用于双层土坡的稳定性评价。

图1 坡底破坏机制示意图

假设滑动破坏区A-B-C-D-A绕旋转中心O做刚体旋转，而对数螺旋面A-E′-D以下土体保持静止，因此，A-E′-D面是一个薄层的速度间断面。设边坡的高度为H，双层土分界面E-E′以上土体高度为αH，基准线OA、OD和OE′的倾角分别为θ0、θh和θm，基准线长度分别为r0、rh和rm，边坡的坡角为β，∠BDC为β′，AB的长度为L。考虑土体重度变化不大，假定双层土坡上下两层土体的重度均为γ；为便于构建统一滑动面，上下两层土体的内摩擦角都为φ；上下两层土体的黏聚力不同，分别为c1和c2。

对数螺旋面的方程为：

r(θ)=r0e(θ-θ0)tanφ

(1)

因此，基准线OD的长度为：

rh=r0e(θh-θ0)tanφ

(2)

在所考虑的滑动破坏区域A-B-C-D-A内，土体重力所做的外功率W，详见Chen在《Limit Analysis and Soil Plasticity》[12]一书，在此仅给出其表达式。

(3)

内部能量耗损主要发生在速度间断面A-E′-D上。能量耗散率的微分，可以由该面的微分面积rdθ/cosφ与滑动面上的黏聚力ci以及该速度间断面上的切向间断速度Vcosφ的乘积计算求得，因此，双层土坡的总的内部能量耗散率可以沿着整个速度间断面积分A-E′-D求得，即：

(4)

令外功率W与内能耗散率D相等，即式(3)与式(4)相等，可以得到：

(5)

化简式(5)，即可得到双层土坡的无量纲的临界高度表达式：

f(θ0，θh，β′，H/r0,L/r0)

(6)

其中：f为目标函数；θ0，θh，β′，H/r0,L/r0为待搜索优化的变量。

1.2 基于坡面破坏机制双层土坡稳定性极限分析方法

当滑动面通过双层土坡坡脚以上时，土坡将发生坡面破坏，见图2。此时，通过引入高度变量h(h

图2 双层土坡坡面破坏机制

这里，引入变量n，并假设：

(7)

此时，由几何关系得：

(8)

A-B-D滑动破坏区域土重所做的外功率为：

(9)

内部能量耗损主要发生在速度间断面A-E′-D上。能量耗散率的微分，可以由该面的微分面积rdθ/cosφ与滑动面上的黏聚力ci以及该速度间断面上的切向间断速度Vcosφ的乘积计算求得，因此，双层土坡的总的内部能量耗散率可以沿着整个速度间断面积分A-E′-D求得，即：

(10)

令外功率W与内能耗散率D相等，即式(9)与式(10)相等，则有：

化简这个等式，则有：

(11)

其中：f为目标函数；θ0，θh，β′，h/r0，L/r0为待搜索优化的变量。

1.3 双层土坡上限解优化计算过程

通过上面基于双层土坡的坡底破坏和坡面破坏机制建立的能量平衡方程，得到双层土坡的无量纲的临界高度γH/c1的表达式：

(12)

其中：f是目标函数。公式(12)给出了双层土坡临界高度γH/c1的一个上限。为了得到γH/c1的上限解，需要对目标函数f进行优化计算，求其最小值，从而获得γH/c1的最小上限解。

目标函数f有五个待优化的参数，为了提高计算效率，本文将采用陈祖煜[17]院士提出的随机搜索方法。对于给定的边坡参数值(坡高H，坡角β，土体黏聚力比值c2/c1，以及上层土的厚度αH)，通过对优化参数的初始值，初始搜索步长以及精度进行设定，进行迭代计算求解，从而得到目标函数f的最小值。优化过程中，角度(θ0，θh，β′ )的精度为0.01°，比值(H/r0，L/r0，h/r0)的精度为0.001。

2 双层土坡临界高度与分层深度关系

依托MATLAB数学计算平台，针对不同破坏机制编制相应的边坡稳定程序进行优化计算，综合考虑各个影响因素的相互影响，通过整理计算结果，绘制了不同因素下双层土坡的临界高度与分层深度的关系图(见图3)。

图3 双层土坡临界高度图(φ = 20°)

从图3可以看出，对于均质土坡(α=1)，当已知边坡的几何参数和土性参数时，其临界高度为一定值，这与前人研究结果一致。对于双层土坡，当边坡的坡度β、坡内土体的内摩擦角φ以及上下两层土的黏聚力比值c2/c1一定时，双层土坡的临界高度γH/c1随分层位置α的变化呈现不同的规律性。

3 不同因素对双层土坡最危险滑动面的影响

土坡发生失稳破坏时存在三种不同的破坏机制，不同破坏机制下其滑动面的位置具有明显的差异性，而土坡的破坏机制取决于其几何参数和坡内土体的力学参数。因此，针对双层土坡在不同的几何参数和力学参数影响下，研究其临界滑动面位置的分布规律，旨在为工程实践提供一定的参考。图4和图5为内摩擦角φ=5°时双层土坡在不同的坡角β、上下土层不同分界面位置α及其黏聚力比值c2/c1下，其临界滑动面的位置分布图。

图4 α = 0.50时双层土坡滑动面位置分布图(φ= 5°)

从图4中可以看出，随着黏聚力比值c2/c1的增大，双层土坡临界滑动面的位置越来越接近于坡面，且对于极缓坡度(β不大于15°)的边坡，主要发生坡底破坏。随着坡角β和黏聚力比值c2/c1的增大(即下层土较上层土坚硬程度的增加)，边坡由坡底破坏逐渐发展为坡脚破坏，这在图4(b)中表现的尤为明显。随着坡角β的继续增加，边坡主要发生坡脚破坏。然而当坡角增大到一定程度(坡角β不小于60°)时，随着黏聚力比值c2/c1的增加，边坡将在两层土的分界面处发生坡面破坏，由图4(c)中可以清晰的看出，当c2/c1≥ 5时，这些边坡发生了坡面破坏。

图5为α=1.0时双层土坡滑动面位置分布图。坡底破坏主要出现在较缓边坡，较陡边坡主要发生坡脚破坏，且随着下层土坚硬程度的增加，双层土坡的破坏机制发生改变；对于上硬下软型土坡，随着坡角β的增加，对于较陡边坡也可发生坡底破坏。

图5 α=1.0时双层土坡滑动面位置分布图(φ= 5°)

对比图4可以看出，当分层深度α发生变化时，相同坡角β和黏聚力比值c2/c1下的双层土坡的破坏形式发生了改变，这在图4(b)和图5(b)中表现的较为明显，同时也说明了分层深度α的变化改变了双层土坡的临界失稳机制。对于该边坡，随着分层深度α的增加，土坡由坡脚破坏转变为坡底破坏。

同理，通过研究内摩擦角φ=10° 时双层土坡在不同的坡角β、上下土层不同分界面位置α及其黏聚力比值c2/c1下，其临界滑动面的位置分布，当内摩擦角φ增大到10°时，对于很缓的边坡(如β<15°)，其破坏形式基本上均为坡脚破坏，且随着黏聚力比值的增加，其滑动面的位置基本保持不变，但对于垂直边坡，当黏聚力比值很大(c2/c1≥5)时，双层土坡将发生坡面破坏。

随着分层深度的变化，当α为0.5时，当黏聚力比值很大(c2/c1≥5)时，除坡度很缓的边坡(β<15°)，双层土坡均可发生坡面破坏，其他规律与之前具有很好的一致性。

当分层深度在坡脚处(α=1.0)时，与之前规律不同的是，对于很陡(β≥75°)的边坡，无论是上软下硬型双层土坡还是上硬下软型双层土坡，均发生坡脚破坏，其滑动面位置基本上不随黏聚力比值的变化而发生改变。对于其他坡度的上硬下软型双层土坡，则有会发生坡底破坏。

4 结 论

本文推导了双层土坡不同破坏机制下临界高度的计算公式，编制了相关程序进行双层土坡稳定性计算；在不同几何参数和力学参数条件下，对双层土坡滑动面的分布规律进行了探究，得出以下结论：

(1) 对于上硬下软型土坡(c2/c1<1)，双层土坡的临界高度γH/c1随分层深度α的增加而线性增加，也就是说临界高度γH/c1随着上层较硬土层厚度的增加而线性增加。对于上软下硬型土坡(c2/c1>1)，双层土坡的临界高度γH/c1随分层深度α的增加而减小，即临界高度γH/c1随着上层较软土层厚度的增加而减小。

(2) 坡底破坏主要出现于内摩擦角较小(φ<10°)的较缓(β<45°)边坡，但对于上硬下软型双层土坡，在内摩擦角较大且分层深度在坡脚位置处时，也可发生坡底破坏；坡面破坏主要发生在分层深度在坡面中部及中部以下位置处的上软下硬型土坡；坡脚破坏在不同分层深度和不同力学参数下的双层土坡中均有出现。

(3) 双层土坡临界滑动面随着黏聚力比值c2/c1的增加越接近与坡面，随着分层深度的增加和坡角β的增加越来越容易出现坡面破坏，且临界滑动面随内摩擦角φ和坡角β的增加，不同黏聚力比之下的滑动面的位置愈加接近。