非均匀边坡稳定的光滑有限元极限分析

Limit Analysis of Smoothed Finite Element of the Stability of Inhomogeneous Slope

摘要: 非均匀层状边坡稳定是岩土工程中的重要问题之一。基于理想塑性上限和下限的极限分析方法是分析此类问题的常用方法之一。近期提出的基于混合常应力–光滑应变单元的极限分析方法(MCSE-LA)在克服体积锁定和计算性能(精度、效率以及收敛性)方面已被证明比传统有限元极限分析方法具有一定的优势，并且可以同时获得极限状态下的应力场和速度场。基于此，首先通过与经典有限元极限分析对比分析，验证MCSE-LA在非均匀边坡稳定分析中的可行性；其次，针对含软弱层边坡进行参数分析，发现软弱层的黏聚力对安全系数的影响更显著，而软弱层的摩擦角则对破坏机构特征的影响更明显；同时说明自适应的MCSE-LA可以较好地分析实际边坡的稳定性，为边坡工程的稳定性评价提供了一种有效的方法。

Abstract: The stability of inhomogeneous layered slopes is one of the crucial issues in geotechnical engineering. Among the commonly used methods for analyzing such problems is the limit analysis approach based on the upper and lower bounds of plasticity. Recently, a novel limit analysis method (MCSE-LA) based on the mixed constant stress-smoothed strain element has been proposed. It has been demonstrated to have certain advantages over traditional finite element limit analysis methods in overcoming volume locking and computational performance (accuracy, efficiency, and convergence), and can simultaneously provide stress and velocity fields at the limit state. Firstly, by comparing with classical finite element limit analysis, the feasibility of MCSE-LA in analyzing the stability of inhomogeneous slopes is verified. Secondly, parameter analysis is conducted for slopes containing weak layers, revealing that the cohesion of weak layers has a more significant impact on the safety factor, while the friction angle of weak layers has a more noticeable effect on the characteristics of failure mechanisms. It is also noted that the adaptive nature of MCSE-LA can effectively analyze the stability of actual slopes, providing an efficient method for evaluating slope stability in slope engineering.

1. 引言

非均质边坡的岩土体类型及其性质可能在较小的范围内差异很大 [1] ，导致此类边坡可能存在多级潜在滑动面，使边坡治理和监测范围很难确定 [2] 。因此，对非均质边坡的稳定性分析具有重要的现实意义。

以塑性上限和下限定理 [3] 为基础的极限分析方法是边坡稳定分析的一类常用方法之一，Chen [4] 系统地阐述了极限分析的基本原理及其在岩土工程中的应用，并利用变分法证明了二维均质边坡的最危险滑动面是一条对数螺旋线。在此基础上，Kumar和Samui [5] 提出了二维水平层状非均匀边坡稳定性的分析方法，假设每一层的潜在滑动面仍旧为对数螺旋线，不同土层的对数螺线滑动面具有相同的旋转中心，并且首尾相连构成一条过坡脚的分段对数螺旋线。类似地，孙志彬等 [6] 为了克服对数螺旋线在分析内摩擦角非均匀边坡中的局限性，提出了一种离散滑动面的上限破坏机构，将对数螺旋滑动面离散为一系列线段单元，采用点到点的方式生成运动许可的速度间断面，分析了内摩擦角、黏聚力随深度线性变化边坡的临界高度的变化规律。田润泽等 [7] 也在组合对数螺旋线破坏机构的基础上，提出了两层土坡临界高度的极限分析上限解。在上述方法中，滑动面上滑坡体的破坏结构均为绕一个中心的转动，对于更加复杂的边坡搜索最危险滑动面存在较大的困难。如存在软弱层的边坡，滑动面极有可能沿着软弱层，此时绕一个中心的转动破坏机构就不再适合 [8] 。另一方面，Michalowski [9] 提出了垂直条分的上限破坏机构和陈祖煜 [10] 提出了在斜条分上限破坏机构的基础上，通过建立条块之间的运动许可速度场，不再需要假设绕一个中心的转动机构，滑动面可以假设为任意形状，因此可以用于非均匀边坡等复杂的情况，并且这一方法已被推广到三维情况 [11] [12] [13] 。在上述两种思路的基础上，Huang等 [8] 针对含软弱薄层的非均匀边坡提出了转动–平移的上限破坏机构，在转动机构与平移机构之间设置速度间断，实现两种破坏机构的组成模式。Zuo等 [14] 则提出了一种多旋转中心的组合对数螺旋转动上限机构，滑动体被斜条分为一系列刚性块体，每个块体都绕一个独立的旋转中心转动，组成多块体的对数螺旋线转动破坏机构。

在以上的方法中通常需要假设边坡失稳时的破坏机构，即需要假设潜在滑动面位置、形状等几何特征以及滑坡体的运动特征，并且搜索最危险滑动面需要求解一个非线性数学规划问题。对于实际边坡，特别是复杂的非均匀边坡，这将是一个比较困难的问题。数值极限分析方法的提出有效地弥补了这些不足。目前数值极限分析方法已经被纳入《建筑边坡工程技术规范》(GB 50330-2013) [15] 作为评价破坏机构复杂的边坡稳定的方法。其中，有限元极限分析方法是较为常用的方法，如，Sloan和Kleeman [16] 提出的一类带间断面的特殊三节点三角形单元用于极限分析，克服了体积锁定问题，但其计算精度依赖于间断面的布局。Makrodimopoulos和Martin [17] 证明了六节点线性应变单元不仅可以克服体积锁定，还可以获得严格上限解。有限元极限分析已经成功地被用于分析非均匀边坡等复杂边坡的稳定分析 [18] [19] ，但此类高阶单元显著地增加了离散难度和计算自由度，降低了计算效率。

Lim等 [20] 根据凸规划对偶理论提出了根据基于静力平衡方程弱形式的“极限荷载乘子极大法()”与基于运动许可速度场的“极限荷载乘子极小法( $\alpha \to \min$ )”等价，均可以获得结构破坏的严格上限解。两种方法相比，具有变量少、形式简洁且无需塑性内能耗散率函数的解析表达式等优点。但在该方法中需要同时插值应力和速度场的混合单元 [21] 。

为解决上述困难，周锡文等 [22] 在光滑有限元 [23] 的基础上提出了一类新型混合常应力–光滑应变单元(Mixed Constant Stress-Smoothed Strain Element, MCSE)，一方面继承了光滑有限元克服体积锁定、精度和计算效率高的优点，还可以同时离散应力和速度场(位移场)，继承“方法”的优点。此类混合单元已被用于极限分析 [22] 和弹塑性增量分析 [24] ，已被证明 [25] [26] 此类单元不仅少于六节点三角形单元的计算规模，而且计算精度高于这类高阶单元。在MCSE极限分析中，可以根据凸规划对偶理论，从原问题中获得极限状态下的应力场和极限荷载乘子，从对偶问题中获得运动许可速度场，再结合基于最大塑性剪应变率的自适应加密算法，在塑性区细化网格，不仅可以提高计算精度，还可以较好地获得边坡破坏机构。

本文首先通过与不同方法的对比，说明自适应MCSE极限分析(A-MCSE-LA)在非均匀边坡稳定性的可行性。在此基础上分析含软弱层边坡的破坏机构特征及其影响因素。

2. MCSE极限分析方法

2.1. 极限分析的广义变分原理

根据极限分析的广义变分原理 [27] ，理想塑性体Ω达到极限状态后，运动许可速度场和静力许可应力场满足以下约束变分问题的驻值条件：

$\begin{array}{l}\underset{\stackrel{\dot{}}{u}}{{\displaystyle \min }}\underset{\alpha ,\sigma }{{\displaystyle \max }}\alpha +{\displaystyle {\int }_{\Omega }{\sigma }^{\text{T}}}\nabla \stackrel{\dot{}}{u}d\Omega -\alpha \left({\displaystyle {\int }_{S_{\sigma }}{t}^{\text{T}}\stackrel{\dot{}}{u}dS}+{\displaystyle {\int }_{\Omega }{b}^{\text{T}}}\stackrel{\dot{}}{u}d\Omega \right)\\ -\left({\displaystyle {\int }_{S_{\sigma }}{t}_0^{\text{T}}\stackrel{\dot{}}{u}dS+{\displaystyle {\int }_{\Omega }{b}_0^{\text{T}}}\stackrel{\dot{}}{u}d\Omega }\right)\\ \text{s}\text{.t}.F\left(\sigma \right)\le 0\end{array}$ (1)

式中： $\nabla$ 是梯度算子， $\sigma$ 是应力， $F\left(\sigma \right)$ 为屈服函数， $\stackrel{\dot{}}{u}$ 是速度，N是应力边界的外法线方向向量矩阵，b₀和t₀为初始不变外力荷载，b和t为可变外力荷载， $\alpha$ 是极限荷载乘子，即当可变荷载增大 $\alpha$ 倍后结构进入极限状态。

在一般有限元离散Ω的基础上，单元内任一点出的速度 $\stackrel{\dot{}}{u}$ 、应力 $\sigma$ 和应变 $\epsilon$ 可以由单元节点速度 $\widehat{\stackrel{\dot{}}{u}}$ 和应力 $\hat{\sigma }$ 插值表示，即：

$\stackrel{\dot{}}{u}\approx N_u\widehat{\stackrel{\dot{}}{u}},\sigma \approx N_{\sigma }\hat{\sigma },\epsilon =\nabla \stackrel{\dot{}}{u}\approx B_u\widehat{\stackrel{\dot{}}{u}}$ (2)

式中： $B_u=\nabla N_u$ 为应变梯度矩阵， $N_u$ 为速度插值形函数， $N_{\sigma }$ 为应力插值形函数。在式(2)的基础上可以得到式(1)的离散形式，并对速度场求变分后可以得到：

$\begin{array}{l}\max \alpha \\ \text{s}\text{.t}\text{.}{B}^T\hat{\sigma }=F_0+\alpha \cdot F\\ F\left(\sigma \right)\le 0\end{array}$ (3)

式中：

$\begin{array}{l}B={\displaystyle {\int }_{\Omega }{N}_{\sigma }^{\text{T}}}{B}_{u}d\Omega \\ F_0={\displaystyle {\int }_{S_{\sigma }}{N}_u^{\text{T}}{t}_{0}dS}+{\displaystyle {\int }_{\Omega }{N}_u^{\text{T}}}{b}_{0}d\Omega \\ F={\displaystyle {\int }_{S_{\sigma }}{N}_u^{\text{T}}tdS}+{\displaystyle {\int }_{\Omega }{N}_u^{\text{T}}}bd\Omega \end{array}$ (4)

优化问题(3)即是文 [21] 中的“ $\alpha \to \max$ 方法”，虽然在形式上是基于静力许可应力场的下限问题，但它是由式(1)离散后经一次变分得出的，故式(3)的最优解既满足运动许可速度场方程，又满足静力许可速度场方程。根据凸规划对偶理论，采用原对偶内点算法时，式(3)的对偶问题和基于速度场的上限问题( $\alpha \to \min$ )是等价的，因此采用原对偶内点算法可以在式(3)的对偶解中获取速度场和塑性乘子的信息。

2.2. MSCE极限分析的二阶锥规划模型

在两级光滑边域有限元(也即是Selective Edge-Based SFEM, sESFEM) [28] 的基础上，文 [22] 提出了MCSE单元(如图1所示)，故此类单元具有sESFEM免疫体积锁定和高精度的优点，详细分析参见文 [28] 。首先采用三节点三角形单元将问题域离散，构成一级网格，如图1(a)所示。在一级网格的基础上，如图1(b)所示的与P₁节点相邻的6个三角形单元，蓝色方形点为三角形中心(如P₂和P₄)，将其与三角形的三个顶点相邻组成三个子三角形单元。然后，以一级网格中的每条边为基础，如图1(b)中所示的第k条边，将与其相邻的两个子三角形单元作为应变光滑域，如图1(b)和图1(c)中所示的P₁、P₂、P₃、P₄围成的阴影部分，即是第k条边的光滑域 ${\Omega }_k^s$ ，其边界为 ${\Gamma }_k^s$ 。由此，在问题域内形成了 $N_s$ 个互不重叠的光滑边域。

Figure 1. MSCE mixed element

图1. MSCE混合单元

在光滑域 ${\Omega }_k^s$ 中，光滑应变为 [28] ：

${\tilde{\epsilon }}_k=\frac{1}{A_k^s}{\displaystyle \underset{i=1}{\overset{N}{\sum }}{A}_i}{\epsilon }_i=B_k^s\widehat{\stackrel{\dot{}}{u}}$ (5)

式中： ${\tilde{\epsilon }}_k$ 为光滑应变， $B_k^s$ 为光滑应变梯度矩阵，N = 1或2，为组成第k个光滑域的子三角形单元的个数，A_i为子三角形单元的面积， $A_k^s$ 为光滑域的面积。与sESFEM不同的是，在MCSE中，假设光滑域内的应力为常数，即：

$N_{\sigma }\hat{\sigma }=\hat{\sigma }={\sigma }_k=constant$ (6)

此外，若假设式(3)服从Mohr-Coulomb屈服准则，则可以表述为一系列的二阶锥约束方程，再利用MCSE单元将式(3)离散，则MCSE-LA最终可以表示为以下的二阶锥规划问题(SOCP) [22] ：

$\begin{array}{l}\max \left[\begin{array}{ccc}1& 0& {\stackrel{\dot{}}{u}}_0^{\text{T}}\end{array}\right]\left[\begin{array}{c}\beta \\ Z\\ r\end{array}\right]\\ \text{s}\text{.t}\text{.}\left[\begin{array}{ccc}-F& K_B& -E_u\\ 0& K_D& 0\end{array}\right]\left[\begin{array}{c}\beta \\ Z\\ r\end{array}\right]=\left[\begin{array}{c}{F}_0\\ k\end{array}\right]\\ Z\in {\mathcal{K}}_q\end{array}$ (7)

式中：Z为等效偏应力、平均应力与应力偏张量组成的向量，其余变量的具体表达及相应的解释可参考文 [22] 。根据凸规划对偶理论，式(7)的对偶问题为：

$\begin{array}{l}\min \left[\begin{array}{cc}{F}_0^{\text{T}}& k^{\text{T}}\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]\\ \text{s}\text{.t}\text{.}\left[\begin{array}{cc}-F^{\text{T}}& 0\\ K_B^{\text{T}}& K_D^{\text{T}}\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]+s_n^x=\left[\begin{array}{c}1\\ 0\end{array}\right]\\ s_n^x\in {\mathcal{K}}_q\end{array}$ (8)

式中： $y_1$ ， $y_2$ 和 $s_n^x$ 与Z向量互为对偶变量，并存在如下关系 [22] ：

$y_1=-U,y_2={\left[A_1^s{\lambda }_1\cdots A_k^s{\lambda }_k\cdots A_{Ns}^s{\lambda }_{Ns}\right]}^{\text{T}}$ (9)

式中：U为整体节点速度向量， ${\lambda }_k$ 为第k个光滑域内的塑性乘子。采用原对偶内点算法求解MCSE-LA时，可以同时求出极限状态下应力场及相应的极限荷载乘子，同时在其对偶问题中可获得运动许可速度场和塑性乘子，进一步可得出应变变化率及塑性内能耗散。

为提高MCSE-LA的计算精度，文 [22] 中以最大塑性剪应变率作为后验误差估计指标，然后再采用开源软件Triangle [29] ，根据每个单元的面积约束进行网格自适应加密，具体算法流程可参见文 [22] 。自适应的MCSE-LA方法的计算效率、精度和收敛性方面的性能在文 [22] 中验证。

3. 非均匀边坡稳定性分析

在文 [22] 中已经将MCSE-LA方法用于分析了典型二维边坡稳定性分析，本文首先进一步将此新方法与不同的方法进行对比分析，并分析软弱层力学性质对非均匀边坡稳定性的影响，最后采用此方法分析在一个实际边坡稳定性。以下问题中，安全系数按强度折减定义如下：

$c^{\prime }=\frac{c}{F_s},{\phi }^{\prime }=\arctan \left(\frac{\tan \phi }{F_s}\right)$ (10)

在本文使用的MCSE-LA方法中，采用的是加载的方式使边坡进入极限状态，即当边坡的重力极限载荷乘子 $\alpha$ 时，边坡即将失稳。本文采用文 [22] 中的二分法搜索边坡的安全系数F_s。

3.1. 两土层边坡稳定性分析

采用MCSE-LA的二阶锥规划方法分析图2所示的一个两土层边坡的稳定性。上层和下层的抗剪强度参数和容重分别为 $c_1,{\phi }_1,{\gamma }_1$ 和 $c_2,{\phi }_2,{\gamma }_2$ ，两土层之间分界面坡度为 ${\beta }_1$ ，边坡坡度为 $\beta$ ，坡高为H，两土层左侧高度分别为H₁和H₂。边界条件如图2所示，左右边界为对称边界，底部为固定边界。

Figure 2. A two-layer slope

图2. 两土层边坡

此类两土层边坡已在不少研究 [5] [7] [14] [30] [31] 中采用了不同的方法进行研究。采用MCSE-LA方法对文 [7] 中算例1的整体稳定进行分析，边坡高度H = 8 m，坡度为1:1.3，层高比H₁/H = 0.5，第一层土的强度参数为c₁ = 5.3 kPa， ${\phi }_1={23}^{\circ }$ ， ${\gamma }_1=19.5\text{kN}/{\text{m}}^3$ ，第二层土的强度参数为：c₂ = 7.2 kPa， ${\phi }_2={20}^{\circ }$ ， ${\gamma }_2=19.5\text{kN}/{\text{m}}^3$ 。图3(a)为本文使用的MCSE-LA自适应计算的网格及速度矢量图，图3(b)为计算获得的剪切塑性耗散图，由此图可以看出此两层边坡的失稳机构与文 [7] 中的分析结果基本一致。

Figure 3. Calculation results of MCSE-LA model (Unit: m)

图3. MCSE-LA模型计算结果(单位：米)

表1中列出了文 [7] 中不同方法的分析结果、采用OptumG2中的上限单元、下限单元和六节点高阶单元以及本文方法计算的结果。通过对比发现，本文方法与其它方法计算的结果基本相同，介于OputmG2计算的上限解与下限解之间，更接近下限解。综合图3和表1的结果可见，MCSE-LA方法可以有效地分析二维层状非均匀边坡的稳定性，不仅可以更准确地计算此类边坡的安全系数，还可以合理地判断边坡的失稳机构、滑动面的形状及位置。

Table 1. Comparison of safety factors with different methods

表1. 不同方法安全系数对比

3.2. 含软弱层的非均匀边坡

采用文中方法进一步分析文 [17] 中给出的一个含软弱夹层的非均质边坡，其几何模型、土体抗剪强度指标及其重度见图4和表2，显然中间土层属于软弱层。在MCSE-LA模型中边坡的边界条件与3.1中的相同，即两侧为对称边界，底部为固定边界。

Figure 4. A inhomogeneous slope with weak layer

图4. 含软弱层非均匀边坡

Table 2. Shear strength parameters of soils

表2. 土体抗剪强度参数

该模型已被前人采用不同的方法进行了研究分析，如Spencer法 [32] 、Bishop法 [18] 、有限元上限法 [18] 、有限元下限法 [18] 、刚体有限元上限法 [33] [34] 以及非连续面拓扑优化方法(DTO) [35] 等。本文使用MCSE-LA分析这个问题时，按两种折减方案进行分析，即：1) 将三层土体按式(10)全部进行折减的整体折减策略；及2) 仅将软弱层的土体2按式(10)进行折减的局部折减策略。图5(a)和图5(b)是两种折减方案的自适应加密网格，图5(c)和图5(d)分别为两种折减方案计算出的整体失稳机构和软弱层失稳机构，这两种失稳机构是由极限状态下剪切塑性内能耗散反映出来。在图5(c)和图5(d)还列出了两种折减方案所计算出的安全系数，整体折减时安全系数为0.42，只折减软弱层时安全系数为0.39。综合图5中的失稳模式及安全系数可知，两种折减方案获得的失稳机构没有明显区别，滑动面的位置和形状基本一致，相应的安全系数差别不大，折减软弱层的方案的安全系数更小，由此设计的边坡支护等治理方案则更安全。其次，通过失稳机构可以看出含软弱层的非均匀边坡的滑动面一般限制在软弱层中，其形状一般不是圆滑型或对数螺旋型，而是不规则形状。

Figure 5. Adaptive meshes and plastic dissipation (Unit: m)

图5. 自适应网格及塑性耗散(单位：米)

表3中列出了不同文献中的分析结果，可以看出，分析方法、失稳机构以及最危险滑动面的搜索算法均会影响安全系数的结果，但结果都在0.40附近。相比之下，基于有限元的极限分析方法不仅无需事先假定失稳机构，而且求极限荷载的优化问题一般都可以表示为线性规划、二阶锥规划等凸规划问题，如本文采用的MCSE-LA方法，最终转化为求解二阶锥规划问题，特别是再结合自适应分析后，还可以根据塑性剪切耗散准确地获得最危险滑动面的位置和形状，并且从表3中的计算结果看，本文的安全系数更接近上限解和下限解的中位值。并且根据文 [22] 的分析，MCSE-LA的计算精度和效率均高于六节点高阶单元。由此可见，自适应的MCSE-LA是评价非均匀边坡的一种有效的、较高精度的实用方法。

Table 3. Comparison of safety factors calculated by different methods

表3. 不同方法计算安全系数对比

在图4和表2所示模型的基础上，采用自适应的MCSE-LA方法进一步分析了软弱层土体黏聚力和内摩擦角对安全系数和破坏机构的影响，如图6和图7所示。图6(a)是在软弱层土体内摩擦角为5˚和重度为18.8 kN/m³情况下(其它两个土层取表2中数值)，软弱层黏聚力从2.94 kPa (即是土体1黏聚力的1/10)到29.4 kPa (此时在土体3以上部分为均匀边坡)时，安全系数的变化情况。从图中可见，软弱层的黏聚强度越高，则安全系数越大。相应地，图7(d)、图7(e)和图7(f)展示了黏聚力逐步增加过程中边坡破坏机构的变化，从图中可以看出软弱层的黏聚力越大，边坡越接近均匀边坡，其滑动面的形状越接近对数螺旋破坏机构，这与文献 [4] 中利用变分法证明的结论一致。

Figure 6. Relationship between strength parameters and safety factors of weak layer

图6. 软弱层强度参数与安全系数的关系

Figure 7. Relationship between strength parameters and failure mechanism of weak layer (Unit: m)

图7. 软弱层强度参数与破坏机构的关系(单位：米)

图6(b)是黏聚力为9.8 kPa、重度为18.8 kN/m³情况下(其它两个土层取表2中数值)，软弱层土体内摩擦角从0˚逐步增加到12˚时安全系数的变化情况。可以看出其变化规律与改变黏聚力时的规律类似，内摩擦角越大，则安全系数越大。相比较而言，黏聚力与安全系数之间关系呈现出较强的非线性，而摩擦角与安全系数的关系曲线则接近线性关系，可能说明黏聚力对安全系数的影响作用更强。图7(a)、图7(b)和图7(c)展示了内摩擦角变化时，边坡的破坏机构的变化。从图中可以看出，当摩擦角为0˚时，滑动面接近圆弧形，这说明纯黏性土边坡失稳时假设圆弧滑动面是合理的，当摩擦角为12˚时，其滑动面呈现出对数螺旋形状，尽管此时土层1与土层2的黏聚力不同，这可能反映出内摩擦角的变化对滑动面性质的影响更明显。

通过上述分析可知，自适应MCSE-LA可以较清晰地分析含软弱层非均匀边坡的稳定性问题，并且无需事先假设破坏机构，仅需要根据实际情况提供土层力学参数、几何模型、本构模型以及边界条件即可计算出可靠的安全系数和破坏机构。

4. 结论

本文将自适应MSCE-LA方法用于非均匀边坡的稳定性分析中，通过层状边坡、含软弱层的非均匀边坡的数值分析，得到以下结论：首先，基于混合常应力–光滑应变的数值极限分析方法可以有效地分析比较复杂的非均匀边坡的稳定性，在结合网格自适应分析后，无需事先假定边坡的破坏机构，可以准确地求出边坡稳定的安全系数，并自动获得边坡的最危险滑动面的位置和形状；其次，对于含软弱层的非均匀边坡，软弱层黏聚力的大小对安全系数的影响更为明显，而摩擦角的大小则对破坏机构的特征影响更为明显；纯黏性土的滑动面更接近圆弧，各层土的摩擦角接近时，无论各层的黏聚力是否接近，其破坏机构均接近对数螺旋线。