Limit analysis method for slope stability based on discretization of rigid blocks
-
摘要: 针对边坡稳定分析这一岩土力学的经典问题，放弃通过极限平衡等式方程和引入假定来直接求解边坡安全系数的传统思路，在建立边坡极限状态机构一般性构建方法的基础上，依据塑性力学上下限定理，以界面力或速度为主要变量，以应力许可静力场或运动许可位移场的基本要求为约束条件，在不引入任何假定的前提下，将求解边坡安全系数的问题转化为求解安全系数上限、下限解的优化数学问题，通过上下限值逐渐逼近安全系数，形成了完整统一的边坡稳定极限分析方法体系。该方法体系由于没有引入任何假定，不仅使边坡稳定分析方法具备了更为严密的理论基础，而且可方便地推广到三维领域。通过6个代表性算例验证了该方法计算结果的准确性、可靠性以及工程应用的合理性、可行性。Abstract: The classic slope stability analysis problem is solved by abandoning the traditional and direct solutions for safety factors of slopes, such as the limit equilibrium equation and the introduction of assumptions. Instead, the general construction method for constructing a slope limit state mechanism is adopted, and the plastic (upper and lower) limit theorems are taken into account. Moreover, taking the interfacial forces and velocities as the main variables and the basic requirements for statically admissible stress field or kinematically admissible displacement field as the constraints, and on the precondition that all assumptions are omitted, the solution for the safety factor problem of slopes is converted into an upper- and lower-bound optimization problem. Furthermore, a complete and united limit analysis method for slope stability is established by gradually shifting the upper- and lower-limit values toward the real solution for the safety factor. To conclude, the proposed method can provide a robust theoretical basis for slope stability analysis due to the omission of assumptions and address the bottleneck resulting from the method of expanding the 2D slope stability analysis to its 3D form. The accuracy and realiability of the calculated results as well as the rationality and feasibility of its engineering applications are validated through 6 representative examples.
-
Keywords:
- slope stability /
- limit analysis /
- upper and lower bound theorem /
- optimization model /
- solution
-
0. 引言
边坡稳定分析是岩土力学的经典问题之一，也是工程界和学术界高度关注的热点。传统的边坡稳定分析方法，采用强度储备的安全系数定义，对滑坡体进行条分离散，通过建立极限平衡方程来求解边坡稳定安全系数。由于未知量（条间力的大小和方向、作用点位置）数目大于可用方程数量，在力学上属于静不定问题，无法直接获得解答[1-4]。
在二维边坡稳定分析领域，传统的方法通常将滑体划分为n个条块，需引入n-1个假定条件，才能使得问题静定可解[2]。200多年来，国内外学者通过引入不同的假定，提出的边坡稳定分析方法有数十种之多，如瑞典圆弧法、毕肖普法、简布法、摩根斯坦–普瑞斯法和萨玛尔法等[5-9]。这些假定的引入不仅破坏了边坡稳定分析方法的理论严密性，而且这些假定在实际应用中是有适用条件的，如果选用的稳定分析方法不当，会导致错误的计算结果，给工程带来严重问题。
将二维方法拓展到三维领域时，这些问题就更为突出。如果参照二维方法将滑体分为n行m列个条柱，则至少需要引入8nm-1个假定[1]，才能使问题静定可解。由于未知量和需引入的假定急剧增加，多数情况下已经失去了实际应用意义。因此，到目前为止三维稳定分析方法仍很不成熟[10-20]。
本文结合水利水电建设中重大边坡工程的实践，逐渐建立了基于塑性力学极限分析理论[21-30]的边坡稳定分析方法体系。在边坡稳定分析中，塑性力学极限分析上下限定理可表述如下：
（1）下限定理：对于一个滑动机构，任意构建一个应力场，只要这个应力场满足静力平衡方程且不违反屈服准则，就是一个静力许可的应力场，那么根据这个应力场通过静力平衡所计算得到的极限荷载，一定不大于边坡所能承受的真实极限荷载，是极限荷载的下限，对应边坡安全系数的下限。
（2）上限定理：对于一个滑动机构，任意构建一个速度场，只要这个速度场满足变形协调条件且不违反流动准则，就是一个运动许可的位移场，那么根据这个位移场通过内能耗散与外力做功相等所求得的极限荷载，一定不小于边坡所能承受的真实极限荷载，是极限荷载的上限，对应边坡安全系数的上限。
依据上下限定理，本文放弃通过极限平衡等式方程和引入假定来直接求解边坡安全系数的传统思路，在建立边坡极限状态机构一般性构建方法的基础上，以边坡极限状态机构界面力或速度为主要变量，以应力许可静力场或运动许可位移场的基本要求为约束条件，在不引入任何假定的前提下，将求解边坡安全系数的问题转化为求解边坡稳定安全系数上限、下限解的优化数学问题，通过上下限值逼近安全系数，形成了完整统一的边坡稳定极限分析方法体系。
1. 极限破坏机构一般性构建方法
边坡极限状态破坏机构的构建是开展边坡稳定性分析的前提和基础。传统的条分法大多采用了简化方法来考虑条块内分界面的影响，这种影响有时是不能忽略的。本文提出的极限状态机构构建方法与传统的条分法不同，采用了能全面考虑边坡地层分布和地质构造特征的块分系统，比传统条分法具有更广泛的适用范围；同时，建立了反映边坡所有界面交切关系的拓扑矩阵以及基于拓扑关系矩阵的块体搜索算法，实现了复杂边坡块分模型的快速生成与自动调整，不仅大大降低了稳定分析前处理工作量，而且所形成的拓扑关系矩阵也为稳定分析优化数学模型，特别是复杂三维模型的自动建立提供了必备条件。传统的条分法可以看作是该方法的一种特例，当不考虑每一条块内部分区界面，采用垂直界面进行块体划分时，块分系统就退化到传统的条分法。图 1为包含各种分区边界的二三维边坡模型及编号系统，图 2为基于矩阵表达式的二维、三维块体搜索算法。
2. 边坡极限分析上下限解优化模型
2.1 下限解优化数学模型
根据极限分析下限定理，在极限破坏机构一般性构建方法基础上，提出了具有统一数学表达形式的二维和三维极限分析下限解优化数学模型。模型中以块体界面作用力的大小、方向和作用点为主要变量，以求解安全系数最大值为目标函数，以满足下限定理许可静力场要求的静力平衡和屈服准则为约束条件，在不引入任何假定的前提下将边坡稳定分析问题转化为求解边坡安全系数下限值的优化数学问题。
（1）稳定安全系数
边坡稳定安全系数（F）采用与传统方法一致的强度储备法定义：
ce=c/F ，fe=tanϕ/F ，} (1) 式中，c，φ为边坡材料的黏聚力和内摩擦角，ce，fe为边坡材料进行折减后的黏聚力和内摩擦系数。
求解安全系数（F）的下限值即为优化数学模型的目标函数。
（2）变量体系
对于离散后的边坡块分系统，待求解的变量除了安全系数（F）外还包括：空间块体每个界面上的法向力（Nj）和切向力（可用两个分量表示，其量值分别为Tlj和Tmj）以及作用点的空间位置rj (rxj，ryj，rzj)。
为了描述界面力的方向特性，需要在每个空间多边形POj上建立一个局部坐标系，该坐标系以界面的法线向量npj为N轴正向，以空间多边形任意一条边界内法线向量tmj=mij为Tm轴正向，i可取为多边形POj的任一边，通常取i=1；Tl轴正向则根据N轴、Tm轴通过右手螺旋法则确定tlj=mij×npj。对于包含该界面的块体，还需进一步根据界面在该块体中的内法线向量nik与界面的法线向量npj的关系来最终确定N，Tm和Tl轴的正向：
Nik=Nj⋅npj⋅flag ，Tlik=Tlj⋅tlj⋅flag ，Tmik=Tmj⋅tmj⋅flag ，} (2) 式中，flag={1,nik=npj−1,nik=−npj。
对于图 3中的两个相邻块体Ba，Bb，两者共用多边形界面POj，POj在Ba，Bb中的局部编号分别为PO2a，PO4b。对于块体Bb而言，界面PO4b的内法线向量n4b与上文定义的界面POj的法线向量npj同向，此时PO4b界面上的局部坐标与POj临时局部坐标相同；而对于块体Ba，界面PO2a的内法线向量n2a则与界面POj的法线向量npj相反，在这种情况下PO2a界面上的局部坐标也应随之反转。根据上述描述，相邻块体Ba，Bb，共用界面POj上的界面力在两个块体中的定义可分别用式（3）和式（4）表示：
N2a=Nj⋅npj ，Tl2a=Tlj⋅tlj ，Tm2a=Tmj⋅tmj ，} (3) N4b=Nj⋅(−npj) ，Tl4b=Tlj⋅(−tlj) ，Tm4b=Tmj⋅(−tmj) 。} (4) ![]()
图 3 相邻块体Ba，Bb共用界面上力的表示方法Figure 3. Common interface forces for adjacent blocks Ba and Bb通过上述处理，每条界面上的界面力只需要一组变量Nj，Tlj和Tmj，就可以准确描述共用界面上作用力与反作用力的关系。
由上述分析可以看出，文中提出的下限解优化数学模型中变量总数为6q′+1，q′为除了临空面之外所有多边形界面的总和。需要指出的是，在本文方法中，除上述变量外，还有如块体重力Gk、外部荷载Ql等已知量会出现在相关的约束方程中。
进一步，如果令每个界面POj局部坐标系中的Tmi=0，并将三维分析变量体系中的各个矢量退化为二维矢量，则三维问题的变量体系就自动退化为二维问题。
（3）约束方程
对于任意一块体（Bk），其力和力矩的平衡方程可以表示为
(h2∑i=1(Nik+Tlik+Tmik)+h3∑i=1Qik+Gk)⋅lp=0， (5) [h2∑i=1(rik−PGk)×(Nik+Tlik+Tmik)+h3∑i=1(PQik−PGk)×Qik]⋅lp=0， (6) 式中，h2为块体Bk包含的界面数量，Gk为块体Bk的重力，PGk为块体重心，Qik为作用于块体Bk上的外荷载，PQik为外荷载作用点，h3为作用在块体Bk上外荷载数量，lp{lx=(1,0,0)ly=(0,1,0)lz=(0,0,1)为投影轴的方向向量。当所有矢量为二维矢量，且lp{lx=(1,0)ly=(0,1)时，上述方程自动退化为了二维静力平衡约束方程。
每一个块体在边界面POj上的作用力应不违反摩尔库仑屈服准则，如下式所示：
√Tlj⋅Tlj+Tmj⋅Tmj⩽ (7) 对于任意界面 P{O_j} 上的作用力，为保证其作用点 {r_j} 位于多边形 P{O_j} 内部，应同时满足以下两个条件：
条件1，作用点在界面上，
({r_j} - a{p_j}) \cdot n{p_j} = 0 \text{，} (8) 条件2，作用点位于界面内，
({r_j} - a_j^i) \cdot m_j^i \geqslant 0{\text{ }}(i = 1,2, \cdots ,{h_1} - 1) \text{，} (9) 式中， a{p_j} 为多边形 P{O_j} 形心， a_j^i 为多边形 P{O_j} 第i个边界的中点。
同时，为反映岩土体不受拉的特性，还需增加式（10）所示的约束条件：
{N_j} \geqslant 0 。 (10) （4）优化数学模型
综合上述分析，可将边坡稳定分析下限解归结为以安全系数最大值为目标函数，以满足极限分析下限定理静力许可应力场要求为约束条件的优化数学模型，优化模型如图 4所示，这是一个超大自由度、高度非线性的复杂优化问题。
2.2 上限解优化数学模型
采用与下限解类似的方法，提出了具有统一数学表达形式的二维和三维极限分析上限解优化数学模型。上限解所采用的变量体系和约束方程与下限解完全不同，具体来讲，上限解以块体界面速度的大小和方向为主要变量，以满足上限定理运动许可位移场要求的速度协调、不违反流动准则和功能平衡方程为主要约束条件，最终建立的优化数学模型与下限解具有类似的形式。
（1）变量体系
上限解对应的变量体系主要包括安全系数（F）、块体速度以及界面相对速度。
空间块体每个界面的相对速度包括法向速度和位于该界面的切向速度（可用两个分量表示），其量值分别为v{N_j}，vT{l_j}和vT{m_j}，如图 5所示。局部坐标及相邻两个界面间的速度描述方式与下限解方法类似，每条界面上的界面速度只需要一组变量v{N_j}，vT{l_j}和vT{m_j}，就可以准确描述共用界面上的相对运动关系。此外为计算内能耗散还需引入虚拟的界面作用力，界面作用力定义与下限解相同。上限解优化数学模型的变量总数为3p+1，p为块体总数。
（2）约束方程
对于任一界面 P{O_j} ，其速度协调方程为
\left.\begin{array}{l}v{N}_{j}=({V}_{j}^{1}-{V}_{j}^{2})\cdot n{p}_{j}\text{ }\text{，}\\ vT{l}_{j}=({V}_{j}^{1}-{V}_{j}^{2})\cdot t{l}_{j}\text{ }\text{，}\\ vT{m}_{j}=({V}_{j}^{1}-{V}_{j}^{2})\cdot t{m}_{j}\text{ }\text{，}\end{array}\right\} (11) 式中，V_j^1，V_j^2为共用界面 P{O_j} 相邻两个块体的速度。
根据流动法则，界面上法向速度与切向速度的夹角不超过内摩擦角，相应的流动法则约束方程为
v{N_j} \leqslant \frac{{\tan {\varphi _j}}}{F}\sqrt {vTl_j^2 + vTm_j^2} 。 (12) 当所有矢量均为二维矢量且VTm_j^{} = 0时，式（11），（12）自动退化为二维问题的速度协调和流动法则约束方程。
根据外力做功与内能耗散相等的要求可建立整个滑坡体的功能平衡方程式：
\begin{aligned} & \sum\limits_{i=1}^{h_{\mathrm{ak}}}\left(V N_k^i \cdot N_k^i+V T l_k^i \cdot T l_k^i+V T m_k^i \cdot T m_k^i\right) \\=& G_k \cdot\left(V N_k+V T l_k+V T m_k\right)+\\ & \sum\limits_{i=1}^{h_q} Q_k^i \cdot\left(V N_q^i+V T l_q^i+V T m_q^i\right), \end{aligned}\text{，} (13) 式中， \sum\limits_{i = 1}^{{h_{{\text{ak}}}}} {(VN_k^i \cdot N_k^i + VTl_k^i \cdot Tl_k^i + VTm_k^i \cdot Tm_k^i)} 为界面的内能耗散， {h_q} 为界面总数； {G_k} \cdot (V{N_k} + VT{l_k} + VT{m_k}) 为重力做功，{G_k}为块体重力； \sum\limits_{i = 1}^{{h_q}} {Q_k^i \cdot (VN_q^i + VTl_q^i + VTm_q^i)} 为外荷载所作的功， {h_q} 为外荷载数量， Q_k^i 为外荷载大小。显然式（13）对于二维和三维问题都是通用的。
另外，为保证块体之间互不嵌入，需增加以下约束条件：
v{N_j} > 0 。 (14) （3）优化模型
集成上述分析成果，可将边坡稳定分析上限解归结为以安全系数最小值为目标函数，以满足极限分析上限定理运动许可位移场要求为约束条件的优化数学模型，如图 6所示，与下限优化模型一样也是一个超大自由度、高度非线性的复杂优化问题。
3. 优化模型的高效求解方法
对于上述复杂的优化模型，已很难采用已有的商用软件直接求解。为解决这一难题，在对非线性程度较高的约束方程进行降维处理的同时提出了一种基于逻辑判断的二分迭代优化求解新算法。
3.1 屈服准则约束方程降维数学表征方式
在优化模型约束方程中，下限解屈服准则约束方程是非线性程度最高的。对屈服准则约束方程通过以下两方面进行了线性化处理。
（1）屈服准则约束方程几何变换降维处理方法
如前所述，下限解中每一个块体在边界面 P{O_j} 上的作用力应满足式（7）所示的屈服准则约束条件，为高度非线性的四次方程。式（7）的几何意义为：切向力 {T}_{j}(T{m}_{j}，T{l}_{j}) 在空间多边形 P{O_j} 上位于以界面力作用点 {r_j} 为圆心、以 {N_j}\frac{{\tan {\varphi _j}}}{F} + \frac{{{c_j}A}}{F} 为半径的圆内，如图 7（a）所示。
为了降低式（7）的非线性程度，将上述圆用其内接等边多边形进行简化逼近，如图 7（a）所示，内接等边多边形第i条边对应的内法线向量可通过式（15）计算得到，如图 7（b）所示。
m{l^i} = ( - \cos {\theta ^i}, - \sin {\theta ^i}){\text{ }}(i = 1,2, \cdots {h_4}) \text{，} (15) 式中， {h_4} 为内接等边多边形的边界数量， {\theta ^i} 为 m{l^i} 与 Tm 轴负方向的夹角，如下所示：
{\theta ^i} = (2i + 1){\text{π }}/{h_4} 。 (16) 根据空间凸多边形点集的定义，式（7）可转化为
({T_j} - a{l^i}) \cdot m{l^i} > 0{\text{ }}(i = 1,2, \cdots ,{h_4}) \text{，} (17) 式中， a{l^i} 为内接等边多边形第i条边界上任一定点。
根据内接多边形几何关系可知 a{l^i} \cdot m{l^i} = - \left( {{N_j}\frac{{\tan {\varphi _j}}}{F} + } \right. \left. {\frac{{{c_j}{A_j}}}{F}} \right) \cdot \cos ({\text{π }}/{h_4}) ，则式（17）可进一步转化为
T{m_j}\cos {\theta ^i} + T{l_j}\sin {\theta ^i} \leqslant \left( {{N_j}\frac{{\tan {\varphi _j}}}{F} + \frac{{{c_j}{A_j}}}{F}} \right) \cdot \cos ({\text{π }}/{h_4}){\text{ }}(i = 1,2, \cdots ,{h_4}) 。 (18) 通过上述转换，屈服准则约束方程由四次方程简化为一系列二次方程。
（2）屈服准则约束方程变量替代降维处理方法
对于上述通过几何变换后的屈服准则约束方程，通过引入水平惯性力（惯性力系数 \eta ）可进一步将二次方程线性化，其具体思路为在静力平衡方程中增加水平惯性力这个虚拟外力 \eta {G_k} ，即
\left( {\sum\nolimits_{i = 1}^{{h_2}} {(N_k^i + Tl_k^i + Tm_k^i) + \sum\nolimits_{i = 1}^{{h_3}} {Q_k^i} + {G_k} + \eta {G_k}} } \right) \cdot lp = 0 。 (19) 此时不再去直接求解边坡安全系数F，而是给定不同的强度折减系数 {F_{\text{c}}} 值（ {F_{\text{c}}} 为常量），优化求得 {F_{\text{c}}} 对应的水平惯性力系数 \eta ，进一步通过迭代求得 \eta = 0 时（相当于没有虚拟外力）的边坡稳定安全系数F。按照这一思路，式（19）变成了一组线性方程，如下所示：
\operatorname{{\mathit{Tm}}}_j \cos \theta^i+T l_j \sin \theta^i \leqslant$ $\left(N_j \frac{\tan \varphi_j}{F(\eta)}+\frac{c_j A_j}{F(\eta)}\right) \cdot \cos \left(\pi / h_4\right) \quad\left(i=1,2, \cdots, h_4\right) 。 (20) 对于上限解中非线性程度最高的流动准则约束方程，采用了与下限解相同的做法，通过在功能平衡方程中引入水平惯性力（式（21）），同样可以将直接求解安全系数的问题变为迭代求解 \eta = 0 的问题。
\sum\nolimits_{i = 1}^{{h_q}} {(V{N_j} \cdot {N_j} + VT{l_j} \cdot T{l_j} + VT{m_j} \cdot T{m_j})} = \sum\nolimits_{k = 1}^{{h_p}} {{V_k}} \cdot {G_k} + \sum\nolimits_{l = 1}^{{h_l}} {{Q_l}} \cdot V{Q_l} + \sum\nolimits_{k = 1}^{{h_p}} {\eta {G_k}{V_k}} \cdot lx \text{，} (21) 式中， lx 为坐标轴方向矢量（1，0，0）。
此时流动准则约束方程中的强度折减系数 {F_{\text{c}}} 为一常量，式（12）由四次方程降为了二次方程式：
V{N_j} \leqslant \frac{{\tan {\varphi _j}}}{{{F_{\text{c}}}}}\sqrt {vTl_j^2 + vTm_j^2} 。 (22) 3.2 基于逻辑判断的二分迭代优化求解算法
如前所述，采用本文提出的边坡稳定极限分析方法求解边坡安全系数，最终归结为一个超大自由度、高度非线性的优化求解问题，为了解决这一难题，进行了降维处理。在降维处理过程中将直接求解边坡安全系数F转变成迭代求解 \eta = 0 的问题，这一过程采用了二分迭代法，如图 8所示。
在二分迭代的每一步迭代过程中，都需要优化求解给定 {F_{\text{c}}} 对应的 \eta 值。虽然采取了降维处理，但 \eta 的优化求解也不是一件容易的事，对于复杂问题仍然存在收敛慢甚至不收敛的问题。
为此，根据 \eta 与 {F_{\text{c}}} 的单调关系以及安全系数F对应于 \eta = 0 的特性，提出了一种适用于本问题的独特算法，即基于逻辑判断的二分迭代优化求解算法。该算法在每一步二分迭代优化求解 \eta 时并不要求其完全收敛，只需通过正负号的逻辑判断就可调整二分迭代的搜索区间，逐步逼近 \eta = 0 时的安全系数F。这一算法不仅规避了优化求解收敛难的问题，而且大大提高了计算效率。图 9以一个典型的三维边坡为例，给出了二分迭代逻辑判断优化求解的具体过程。
为直观说明本文方法的优势，以图 10（a）所示较为简单的楔形体为例，与直接优化求解法进行了对比，如图 10（b）所示。对比结果表明，本文方法比直接优化求解法的计算效率提高了20倍以上，效果十分显著。图 10（c）则显示了本文方法的整个迭代过程，从计算过程来看，本文方法不仅收敛很快而且迭代过程十分稳定。本文所提出的高效优化求解方法，从根本上解决了该方法体系在实际应用中的瓶颈问题，使得三维方法的推广应用真正成为可能。
4. 边坡稳定极限分析方法体系
在上述研究基础上形成了边坡稳定极限分析的完整体系，如图 11所示。该方法体系与传统边坡稳定分析方法不同，采用优化方法分别求出边坡安全系数的上限和下限，通过上下限去逼近安全系数，直接求得安全系数或将安全系数控制在一个较小的范围内。对于一些简单二维问题，上限解和下限解可以完全闭合；而对于一些复杂的三维问题，上限解和下限解可以给出安全系数的取值范围，这比传统方法获得的一个无法预知误差的单一结果具有很大的优势。
5. 算例验证及工程应用
5.1 二维边坡算例验证
图 12（a）为具有经典理论解的边坡计算模型和计算结果，根据索科洛夫斯基滑移线理论计算得到的极限荷载q=1091.417 kPa，理论上讲，如果将上述分布荷载施加于边坡上，则边坡处于极限状态，其安全系数为1.00；图 12（b），（c）分别为均质和非均质边坡ACDAS考题的计算模型和计算结果。应用本文方法并采用与上述算例相同的破坏模式进行了对比分析。
从图 12中计算结果的对比可以看出，采用本文方法计算得到的上限解和下限解均与理论解或裁判推荐答案相同，从一定程度上说明了本文计算方法的合理可靠性；同时也表明，对于一些简单的二维问题，极限分析上限解和下限可以完全闭合。
5.2 三维边坡算例验证
图 13（a）为根据ACDAS二维均质边坡考题经平面拉伸后建立的三维边坡模型（实际上是一个平面应变问题）及其计算结果。本文方法计算得到的上限解和下限解均为1.00，与ACDAS二维考题裁判答案（1.00）一致，从理论上讲是完全正确的。
图 13（b）是一个在边坡稳定分析领域经常引用的具有椭球形底滑面的三维边坡算例，很多学者给出了不同分析方法的稳定分析结果。采用本文方法计算得到的安全系数上下限分别为2.18和2.22，将边坡安全系数限值在一个较小的范围内，其误差范围约为（2.20±0.02）。通过与其它分析方法计算结果的对比可以看出，本文方法计算得到的下限解大于Zhang[31]通过极限平衡法计算得到的结果2.12，（Chen等[2]研究认为极限平衡法是一种具有下限特性的解），小于Chen等[2-3]通过EMU3D计算得到的上限解为2.26，无疑本文方法计算结果更加合理，而且给出了安全系数的误差范围，比其他方法只获得单一安全系数具有更大的优势。
5.3 工程实例
三维工程边坡算例选取了较为复杂的ALG库区变形体，前期勘察结果表明变形体顶部已出现大量与坡面平行的裂缝，且张开较大（部分达到m级的张开度），说明边坡整体安全裕度不高，接近失稳破坏临界状态。变形体整体情况如图 14所示，变形体主要结构面及岩体强度参数如表 1所示。
表 1 ALG变形体主要材料工程力学特性参数Table 1. Parameters of ALG deformation bodyF1 J2 J3 岩石 c
/kPa\varphi
/(°)c
/kPa\varphi
/(°)c
/kPa\varphi
/(°)c
/kPa\varphi
/(°)\rho
/(kg·m-3)20 16.7 70 25.6 70 25.6 175 32 2500 对图 14（b）所示的A—A、B—B、C—C和D—D 4个剖面采用EMU法进行了二维边坡稳定性分析。根据二维稳定分析结果，剖面A—A、B—B计算得到的安全系数分别为1.34，1.27，说明边坡具有一定安全裕度，而剖面C—C、D—D计算得到的安全系数则分别为0.80，0.53，说明边坡已经失稳破坏。显然，根据这样的结果无法对边坡的整体稳定性做出判断，需要考虑实际边坡的三维效应，通过三维分析才能对边坡的稳定情况给出切合实际的评价。
图 14（c）为三维分析所采用的计算模型，图中同时给出了本文方法计算得到的稳定分析结果。从分析结果可以看出，该边坡安全系数上下限分别为1.09和1.05，其均值为1.07略大于1.0，表面边坡已处于临界状态，与实际现状完全吻合，这一结果进一步验证了本文方法的合理性，同时也说明了对于一些三维效应明显的边坡开展三维稳定分析是十分必要的。
6. 结论
（1）提出了边坡极限状态破坏机构的一般性构造方法；基于塑性力学极限分析理论，在不引入任何假定的前提下，建立了二三维边坡稳定安全系数上下限解的统一优化数学模型。
（2）提出了优化模型中高度非线性约束方程的降维数学表达方式，同时引入了复杂优化模型的逻辑判断迭代求解方法，解决了大自由度优化模型求解的收敛难题，提高了优化效率和精度，使得该方法应用于三维边坡稳定分析领域成为可能。
（3）集成上下限方法，建立了通过上下限逼近安全系数真实解、可准确确定计算误差的二三维边坡稳定极限分析方法体系。
（4）本文方法的合理性得到了经典算例的验证，并在实际工程中得到了成功应用。
-
![]()
图 3 相邻块体 {B_{\text{a}}} ， {B_{\text{b}}} 共用界面上力的表示方法
Figure 3. Common interface forces for adjacent blocks {B_{\text{a}}} and {B_{\text{b}}}
表 1 ALG变形体主要材料工程力学特性参数
Table 1 Parameters of ALG deformation body
F1 J2 J3 岩石 c
/kPa\varphi
/(°)c
/kPa\varphi
/(°)c
/kPa\varphi
/(°)c
/kPa\varphi
/(°)\rho
/(kg·m-3)20 16.7 70 25.6 70 25.6 175 32 2500 
下载: 导出CSV
-
[1] 周健, 崔积弘, 贾敏才, 等. 静力触探试验的离散元数值模拟研究[J]. 岩土工程学报, 2007, 29(11): 1604–1610. doi: 10.3321/j.issn:1000-4548.2007.11.002 ZHOU Jian, CUI Ji-hong, JIA Min-cai, et al. Numerical simulation of cone penetration test by discrete element method[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(11): 1604–1610. (in Chinese) doi: 10.3321/j.issn:1000-4548.2007.11.002
[2] CHEN Z Y, WANG X G, HABERFIELD C, et al. A three-dimensional slope stability analysis method using the upper bound theorem part I: theory and methods[J]. International Journal of Rock Mechanics and Mining Sciences, 2001, 38(3): 369–378. doi: 10.1016/S1365-1609(01)00012-0
[3] CHEN Z Y, WANG J, WANG Y J, et al. A three-dimensional slope stability analysis method using the upper bound theorem Part II: numerical approaches, applications and extensions[J]. International Journal of Rock Mechanics and Mining Sciences, 2001, 38(3): 379–397. doi: 10.1016/S1365-1609(01)00013-2
[4] 陈祖煜, 弥宏亮, 汪小刚. 边坡稳定三维分析的极限平衡方法[J]. 岩土工程学报, 2001, 23(5): 525–529. doi: 10.3321/j.issn:1000-4548.2001.05.001 CHEN Zu-yu, MI Hong-liang, WANG Xiao-gang. A three-dimensional limit equilibrium method for slope stability analysis[J]. Chinese Journal of Geotechnical Engineering, 2001, 23(5): 525–529. (in Chinese) doi: 10.3321/j.issn:1000-4548.2001.05.001
[5] JANBU N. Application of composite slip surface for stability analysis[C]// Proceedings of Europan Conference on Stability of Earth Slopes, Stockholm, Sweden, 1954.
[6] BISHOP A W. The use of the slip circle in the stability analysis of slopes[J]. Géotechnique, 1955, 5(1): 7–17. doi: 10.1680/geot.1955.5.1.7
[7] MORGENSTERN N R, PRICE V E. The analysis of the stability of general slip surfaces[J]. Géotechnique, 1965, 15(1): 79–93. doi: 10.1680/geot.1965.15.1.79
[8] SPENCER E. A method of analysis of the stability of embankments assuming parallel inter-slice forces[J]. Géotechnique, 1967, 17(1): 11–26. doi: 10.1680/geot.1967.17.1.11
[9] HOEK E, BRAY J W. Rock slope engineering[M]. London: The Institution of Mining and Metallurgy, 1977.
[10] SARMA S K. Stability analysis of embankments and slopes[J]. Journal of the Geotechnical Engineering Division, 1979, 105(12): 1511–1524. doi: 10.1061/AJGEB6.0000903
[11] WANG X G, LIN X C, LI X, et al. 3D slope stability analysis method based on Pan's maximum principle[J]. Landslides, 2020, 17(5): 1163–1176. doi: 10.1007/s10346-019-01303-7
[12] WANG X G, LIN X C, LI X, et al. 3D slope stability analysis method based on Pan's maximum principle[J]. Landslides, 2020, 17(5): 1163–1176. doi: 10.1007/s10346-019-01303-7
[13] 冯树仁, 丰定祥, 葛修润, 等. 边坡稳定性的三维极限平衡分析方法及应用[J]. 岩土工程学报, 1999, 21(6): 657–661. doi: 10.3321/j.issn:1000-4548.1999.06.005 FENG Shu-ren, FENG Ding-xiang, GE Xiu-run, et al. 3D limit equilibrium method for slope stability and its application[J]. Chinese Journal of Geotechnical Engineering, 1999, 21(6): 657–661. (in Chinese) doi: 10.3321/j.issn:1000-4548.1999.06.005
[14] 郑榕明, 朱禄娟, 谷兆祺. 非对称旋转破坏的三维Bishop边坡稳定算法[J]. 岩土工程学报, 2002, 24(6): 706–709. doi: 10.3321/j.issn:1000-4548.2002.06.007 CHENG Yung-ming, ZHU Lu-juan, GU Zhao-qi. Rotational failure of 3D non-symmetric slope predicted by Bishop's method[J]. Chinese Journal of Geotechnical Engineering, 2002, 24(6): 706–709. (in Chinese) doi: 10.3321/j.issn:1000-4548.2002.06.007
[15] 张均锋, 王思莹, 祈涛. 边坡稳定分析的三维Spencer法[J]. 岩石力学与工程学报, 2005, 24(19): 3434–3439. doi: 10.3321/j.issn:1000-6915.2005.19.005 ZHANG Jun-feng, WANG Si-ying, QI Tao. Three-dimensional spencer method for slope stability analysis[J]. Chinese Journal of Rock Mechanics and Engineering, 2005, 24(19): 3434–3439. (in Chinese) doi: 10.3321/j.issn:1000-6915.2005.19.005
[16] 李同录, 王艳霞, 邓宏科. 一种改进的三维边坡稳定性分析方法[J]. 岩土工程学报, 2003, 25(5): 611–614. doi: 10.3321/j.issn:1000-4548.2003.05.020 LI Tong-lu, WANG Yan-xia, DENG Hong-ke. An improved method for three-dimensional slope stability analysis[J]. Chinese Journal of Geotechnical Engineering, 2003, 25(5): 611–614. (in Chinese) doi: 10.3321/j.issn:1000-4548.2003.05.020
[17] 朱大勇, 钱七虎. 三维边坡严格与准严格极限平衡解答及工程应用[J]. 岩石力学与工程学报, 2007, 26(8): 1513–1528. doi: 10.3321/j.issn:1000-6915.2007.08.001 ZHU Da-yong, QIAN Qi-hu. Rigorous and quasi-rigorous limit equilibrium solutions of 3D slope stability and application to engineering[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(8): 1513–1528. (in Chinese) doi: 10.3321/j.issn:1000-6915.2007.08.001
[18] 郑宏. 严格三维极限平衡法[J]. 岩石力学与工程学报, 2007, 26(8): 1529–1537. doi: 10.3321/j.issn:1000-6915.2007.08.002 ZHENG Hong. A rigorous three-dimensional limit equilibrium method[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(8): 1529–1537. (in Chinese) doi: 10.3321/j.issn:1000-6915.2007.08.002
[19] 陈昌富, 朱剑锋. 基于Morgenstern-Price法边坡三维稳定性分析[J]. 岩石力学与工程学报, 2010, 29(7): 1473–1480. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201007025.htm CHEN Chang-fu, ZHU Jian-feng. A three-dimensional slope stability analysis procedure based on Morgenstern-Price method[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(7): 1473–1480. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201007025.htm
[20] 郭明伟, 葛修润, 李春光, 等. 边坡和坝基抗滑稳定分析的三维矢量和法及其工程应用[J]. 岩石力学与工程学报, 2010, 29(1): 8–20. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201001001.htm GUO Ming-wei, GE Xiu-run, LI Chun-guang, et al. Three-dimensional vector sum method employed in slope and dam foundation stability analyses and its applications to practical engineering[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(1): 8–20. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201001001.htm
[21] 凌道盛, 戚顺超, 陈锋, 等. 一种基于Morgenstern-Price法假定的三维边坡稳定性分析法[J]. 岩石力学与工程学报, 2013, 32(1): 107–116. doi: 10.3969/j.issn.1000-6915.2013.01.016 LING Dao-sheng, QI Shun-chao, CHEN Feng, et al. A limit equilibrium method based on Morgenstern-Price method for 3D slope stability analysis[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(1): 107–116. (in Chinese) doi: 10.3969/j.issn.1000-6915.2013.01.016
[22] MCCOOK D K. Limit analysis and soil plasticity[J]. Soil Science Society of America Journal, 1976, 40: 4. doi: 10.2136/sssaj1976.03615995004000040002x
[23] LE C V, NGUYEN-XUAN H, NGUYEN-DANG H. Upper and lower bound limit analysis of plates using FEM and second-order cone programming[J]. Computers & Structures, 2010, 88(1/2): 65–73. http://www.onacademic.com/detail/journal_1000034067936510_8d83.html
[24] SLOAN S W, KLEEMAN P W. Upper bound limit analysis using discontinuous velocity fields[J]. Computer Methods in Applied Mechanics and Engineering, 1995, 127(1/2/3/4): 293–314. http://www.researchgate.net/profile/Scott_Sloan/publication/222487057_Upper_bound_limit_analysis_using_discontinuous_velocity_fields/links/53d72b280cf228d363ead982
[25] LYAMIN A V, SLOAN S W. Upper bound limit analysis using linear finite elements and nonlinear programming[C]// Civil-Comp Proceedings, Finite Elements: Techniques and Developments. Stirlingshire, UK: Civil-Comp Press, 2002: 61–77.
[26] YIN J, CHEN J, LEE C. Upper limit analysis of stability of rock and soil slopes using rigid finite elements[J]. Chinese Journal of Rock Mechanics & Engineering, 2004, 23(6): 898–905.
[27] LIU F T, ZHAO J D. Upper bound limit analysis using radial point interpolation meshless method and nonlinear programming[J]. International Journal of Mechanical Sciences, 2013, 70: 26–38. doi: 10.1016/j.ijmecsci.2013.01.017
[28] YU S B, ZHANG X, SLOAN S W. A 3D upper bound limit analysis using radial point interpolation meshless method and second-order cone programming[J]. International Journal for Numerical Methods in Engineering, 2016, 108(13): 1686–1704. doi: 10.1002/nme.5273
[29] YUAN S, DU J N. Effective stress-based upper bound limit analysis of unsaturated soils using the weak form quadrature element method[J]. Computers and Geotechnics, 2018, 98: 172–180. doi: 10.1016/j.compgeo.2018.02.008
[30] HUANG M S, FAN X P, WANG H R. Three-dimensional upper bound stability analysis of slopes with weak interlayer based on rotational-translational mechanisms[J]. Engineering Geology, 2017, 223: 82–91. doi: 10.1016/j.enggeo.2017.04.017
[31] ZHANG X. Three-dimensional stability analysis of concave slopes in plan view[J]. Journal of Geotechnical Engineering, 1988, 114(6): 658–671. doi: 10.1061/(ASCE)0733-9410(1988)114:6(658)
-
期刊类型引用(10)
1. 周创兵，姜清辉，姚池，位伟，胡冉. 水利水电工程高陡边坡变形与稳定性研究进展及挑战. 水力发电学报. 2025(01): 1-17 . 
百度学术
2. 郑颖人，张金良，尹德文，邵颖，苏凯，吴昊，张智沛. 基于下限模型的临界滑动面理论与数值解法. 岩土工程学报. 2025(02): 438-442 . 本站查看
3. 张标，瞿凡，蒋毅. 饱水软化效应下砂质泥岩边坡三维稳定性上限分析. 土木工程学报. 2025(04): 108-123 . 百度学术
4. 叶欣玲，司富生，胡华，鲁蕊. 里底水电站库区巴迪滑坡蓄水前后稳定性评价. 西北水电. 2024(02): 19-25 . 百度学术
5. 楼晓明，孙逸玮，张蓟. 软土地基沟渠开挖诱发远处围堰失稳的实例分析. 水利水电技术(中英文). 2024(S1): 151-159 . 百度学术
6. 卢树盛，向能武，孙云志. 高边坡勘察新技术在水利工程中的应用现状、挑战与展望. 中国水利. 2024(14): 44-50 . 百度学术
7. 鲁松林. 基于简化Janbu法的夹炭质泥岩软弱层的溶蚀峰丛斜坡边坡稳定性分析. 价值工程. 2024(24): 125-127 . 百度学术
8. 邱锐，乐巧丽. 基于PLAXIS强度折减法的应变软化边坡稳定性分析. 科技与创新. 2024(18): 104-105+109 . 百度学术
9. 童新建，赵世航，刘传新. 基于流固耦合理论的航道边坡稳定性分析. 黑龙江科技大学学报. 2024(05): 722-728 . 百度学术
10. 葛虎胜，刘峰，苏怀斌，王宇，廖九波，方庆红. 高陡边坡开挖卸荷位移演变规律及边坡稳定性研究. 矿业研究与开发. 2023(06): 83-88 . 百度学术
其他类型引用(6)
计量
- 文章访问数: 234
- HTML全文浏览量: 45
- PDF下载量: 110
- 被引次数: 16
