深圳某黏土边坡稳定的数值极限分析

Numerical Limit Analysis of Stability for a Soil Slope in Shenzhen

摘要: 数值极限分析方法是边坡稳定性评价的主要方法之一，特别是针对复杂非均质边坡的问题具有一定的优势。文中使用OptumG2软件分析了深圳某典型非均质边坡的安全系数和失稳模式，该软件中采用了基于混合应力六节点三角形单元的数值极限分析方法，将边坡稳定性问题转化为二阶锥规划问题进行求解。在基于塑性剪切耗散的自适应方法的基础上，不仅可以获得比较准确的边坡稳定安全系数，而且还可以确定边坡的失稳模式。在此基础上，通过数值参数分析，首先验证了该软件分析实际复杂边坡的可行性，其次通过对比两种不同的放坡治理方案的优劣，发现在安全系数相同的情况下二级放坡方案挖方量相对较少，但二级放坡改变了边坡的应力状态，导致放坡后的边坡的失稳模式可能出现局部失稳模式，也可能是两条对数螺旋线复合的失稳模式，可能会增加边坡治理的难度。由此可见，采用文中方法不仅可以为解决实际复杂边坡工程稳定性评价提供有效的途径，还为滑坡治理提供了方法。

Abstract: The numerical limit analysis is one of the main methods for slope stability evaluation, and it has certain advantages especially for problems with complex inhomogeneous soil slopes. The safety factor and failure mechanism of a typical inhomogeneous slope in Shenzhen is analyzed by using OptumG2 software, and a numerical limit analysis based on a mixed stress six-node triangle element is adopted in the software, and the slope stability problem is recast as a second order cone programming. By using the adaptive analysis based on plastic shear dissipation, not only can an accurate safety factor of slope be obtained, but the failure mechanism of the slope can also be determined. On this basis, through numerical parameter analysis, the feasibility of the software in analyzing actual complex slopes is first verified. Secondly, by comparing the advantages and disadvantages of two different slope reduction schemes, it is found that the excavation volume of the secondary slope reduction scheme is relatively small under the same safety factor, but the secondary slope reduction changes the stress state of the slope, resulting in the instability mode of the slope after the slope reduction may appear in a local instability mode, or it may be a composite instability mode of two logarithmic spirals, which may increase the difficulty of slope control. It can be seen that the method in this paper can not only provide an effective way to solve the stability evaluation of actual complex slope engineering, but also provide a method for landslide control.

1. 引言

塑性极限分析上限和下限方法是边坡稳定分析的常用方法之一，Chen [1]系统地阐述了极限分析的基本原理及其在边坡稳定分析中的应用，并利用变分法证明了二维均质边坡的最危险滑动面是一条对数螺旋线，随后此类方法被用于各种边坡的稳定性分析，判断边坡的失稳模式及安全系数[2]-[4]。作者[5]简要地总结了此类方法在分析非均质边坡稳定性时的不足，在这类方法中通常需要假设边坡失稳时的破坏机构，即需要假设潜在滑动面位置、形状等几何特征以及滑坡体的运动特征，并且搜索最危险滑动面需要求解一个非线性数学规划问题。对于实际边坡，特别是复杂的非均匀边坡，这将是一个比较困难的问题。

基于有限元的数值极限分析方法[6] [7]可以有效地克服上述困难，不少研究也证明了此类方法的成果可以用于分析非均匀边坡等复杂边坡的稳定分析[8] [9]。Krabbenhoft等[10] [11]根据凸规划对偶理论提出了基于静力平衡方程弱形式的“极限荷载乘子极大法($\alpha \to \max$ )”与基于运动许可速度场的“极限荷载乘子极小法($\alpha \to \min$ )”等价，均可以获得结构破坏的严格上限解。两种方法相比，具有变量少、形式简洁且无需塑性内能耗散率函数的解析表达式等优点。在该方法中需要同时插值应力和速度场的混合单元[12]，如六节点混合应力单元(mixed six-node triangular element, mixed T6-FEM) [11]和混合常应力–光滑应变单元(mixed constant stress-smoothed strain element, MCSE) [12]。其次该方法还可以根据凸规划的对偶理论，从原问题中获得极限状态下的应力场和极限荷载乘子，从对偶问题中获得运动许可速度场，基于最大塑性剪应变率的自适应加密算法，在塑性区细化网格，不仅可以提高计算精度还可以较好地获得边坡破坏机构。文[5]中将混合应力有限元极限分析方法应用于分析非均质边坡的稳定性，说明此类方法不仅可以准确地求出非均质边坡稳定的安全系数，还可以在结合自适应分析的基础上自动地获取边坡失稳时的滑动面位置和形状，并发现软弱层的黏聚力的大小对安全系数的影响更为明显，而摩擦角的大小则对破坏机构的特征影响更为明显，其次纯黏性土的滑动面更接近圆弧，各层土的摩擦角接近时，无论各层的黏聚力是否接近，其破坏机构接近对数螺旋线。

由此可见，基于自适应混合应力有限元极限分析方法可用于分析实际工程中复杂边坡的稳定性。并考虑到数值极限分析方法已经被纳入《建筑边坡工程技术规范》(GB50330-2013) [13]作为评价破坏机制复杂的边坡稳定的方法。近年来，基于数学规划方法的有限元极限分析和弹塑性分析已经被集成在一款新型软件OptumG2 [14]中。该软件因其操作简单、易于建模、支持自适应计算且收敛性强等优点，已经在多个领域得到应用，如静力和地震作用下挡土墙的土压力系数分析[15]；非饱和土边坡的强度折减稳定分析[16]；含双桩基础边坡的地震稳定性和破坏模式分析[17]；考虑各向异性和空间变异性的边坡随机有限元分析[18]；基坑开挖稳定性分析[19]以及降雨入渗条件下含裂缝边坡的非饱和土边坡稳定性分析[20]等。因此，OptumG2软件是解决复杂边坡稳定性分析的有效途径之一。在边坡工程实践中，作者常常遇到复杂边坡稳定性评价和处理问题。例如，在深圳市龙岗区的一个非均质边坡中，采用传统极限平衡的条分法时常常遇到最危险滑动面难以确定的问题。根据上述综述分析，本文将采用OptumG2软件中的数值极限分析方法进行建模分析，以验证其在解决此类问题时的可行性。

此外，根据《建筑边坡工程技术规范》(GB50330-2013)，当工程场地有放坡条时，且无不良地质作用时宜优先采用坡率法[13]。本文通过参数分析，对比一级放坡和二级放坡两种坡率法方案，证明了当边坡高度较高时可采用多级放坡的坡率法。进而为复杂边坡稳定性分析和滑坡治理提供参考。

2. 工程概况

该边坡位于深圳市龙岗区风门路以东，风门坳科技园后侧，坡面裸露，已经发生过崩塌现象。根据钻探结果显示，场地内勘探深度范围内揭露的地层主要由4大层：上覆第四系人工填土层(Q4ml)、第四系坡积层(Qdl)、第四系残积层(Qel)、下伏基岩白垩系花岗岩风化带(K1zh)。

Figure 1. A typical engineering geological profile of slope in Fengmenao Science Park, Longgang District, Shenzhen

图1. 深圳市龙岗区风门坳科技园后侧边坡典型工程地质剖面

Table 1. Material properties of soils

表1. 土体抗剪强度参数

地层	黏聚力c/kPa	内摩擦角φ/(˚)	重度γ/kN·m^−3
含砾黏土	20.0	16	18.5
砾质黏性土	24.0	20	19.8
全风化花岗岩	25.0	28	21.5

选取图1所示的一个典型地质剖面作为分析对象，第一层为含砾黏土，黄色、褐黄色、褐红色，稍湿。可塑状，以黏性土粒为主，含有少量石英颗粒，含量约为5%~10%，遇水易软化。第二层为砾质黏性土，黄色、褐黄色、褐红色，可塑~硬塑状，矿物颗粒主要为石英，石英颗粒粒径在2 mm~5 mm之间，含量约为30%，岩芯呈土柱状，遇水易软化崩解，该层所有钻孔均有揭露，揭露厚度在8.70 m~22.60 m之间，平均厚度为15.91 m。图1中点1至点6是典型边坡剖面的边界，点10→9→8→7→5→11→12→13→14构成的折线作为第一层含砾黏土与第二层砾质黏性土的分界面，点15→16→17→18→19→20→21构成的折线作为第二层与第三层全风化花岗岩的分界面。根据勘察报告，各层土体的抗剪强度和重度的力学参数见表1。

3. 基于OptumG2的极限分析模型

3.1. OptumG2中的Mixed-T6-FEM极限分析模型

OptumG2 [14]是一款先进的岩土工程设计和分析有限元分析软件，其主要功能包括边坡稳定性分析、渗流和地下水流动、基坑开挖、堤坝、基础和挡土墙的建模与分析。OptumG2的优势在于其高效精确的计算能力、易于建模的用户界面、自适应网格功能以及能够提供上下限解的独特分析方法。此外，OptumG2还支持多种材料模型和分析类型，包括弹塑性分析和极限分析，使其成为解决复杂工程问题的理想工具。

Figure 2. Mixed-T6-FEM

图2. 六节点混合应力单元

在OptumG2软件中，采用了图2所示的mixed-T6-FEM单元进行离散，即基于六个节点进行运动许可速度场的二次插值，而基于三个高斯积分点进行应力线性插值，即单元内任一点出的速度$\dot{u}$ 、应力$\sigma$ 和应变$\epsilon$ 可以由单元节点速度$\widehat{\dot{u}}$ 和应力$\widehat{\sigma }$ 插值表示，即：

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

式中：$B_u=\nabla N_u$ 为应变梯度矩阵，$N_u$ 为速度插值形函数，$N_{\sigma }$ 为应力插值形函数。在式(1)的基础上可以得求极限荷载系数的数学优化问题，即：

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

式中：

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

根据凸规划的对偶理论[10]，采用原对偶内点算法时，式(2)的对偶问题和基于速度场的上限问题($\alpha \to \min$ )是等价的，因此采用原对偶内点算法可以在式(2)的对偶解中获取速度场和塑性乘子的信息。上述过程可在OptumG2中建模求解，并以塑性剪切耗散作为自适应控制变量，较高精度地计算复杂工程的安全系数，同时较好地获得工程结构的失稳模式。

3.2. 参数分析与放坡方案对比

根据多年的工程实践，以强度指标的储备作为安全系数定义的方法被工程界广泛认可，即：

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

在计算安全系数时，一般并非直接求解式(2)，而是以重力最为可变荷载，计算给定抗剪强度指标c和$\phi$ 时的极限荷载乘子$\alpha$ ，若$\alpha >1$ 则增加折减系数的数值，若$\alpha <1$ 则减小折减系数的数值，直到$\alpha \text{=}1$ 时，即为所求的安全系数F_s。

Figure 3. Adaptive mesh and shear dissipation of slope with limit analysis

图3. 边坡极限分析自适应网格及剪切耗散

首先在OptumG2中选择基本MC材料模型，按表1中强度参数和图1中的几何模型建立数值极限分析的模型，并选择关联流动法则，选择K0分析方法计算初始地应力，不考虑拉张截断的影响。边坡数值模型的左右边界设定为对称边界条件，即仅约束法向位移，底部为固定边界条件。在工况管理器中选择6-高斯点单元，初始单元数量为10,000个，并以剪切耗散作为自适应控制变量进行6次的自适应迭代。计算完成后自适应网格和剪切耗散，如图3所示，计算的安全系数为0.533，说明该边坡处于不稳定的状态，这与实际工程情况是吻合的。图3(b)的剪切耗散，反映了边坡失稳时塑性区贯通，也即是滑动面的位置，可以看出该失稳机构属于对数螺线类型，这与作者[5]采用MCSE-LA分析非均质边坡的结果一致。

Figure 4. Two sloping options

图4. 两种放坡方案

通常对不稳定边坡治理最简单、最直接的方案就是放坡，即文[13]中所指坡率法。即将边坡坡角降低，以增强边坡的稳定性。图4中为两种不同的放坡方案，图4(a)为一级放坡，即沿图中红色虚线5→22进行开挖，削去5→22→4→5所围成的三角形区域，从而将坡角降低；图4(b)中则为二级放坡方案，即沿5→22→23→24红色虚线进行分级开挖，削去5→22→23→24→4→5所围成的多边形区域。分级放坡与一级放坡相比，一方面可以控制坡角，另一方面可以控制坡高。分级放坡后实际上是降低了每一级边坡的高度，由此可以较少的降低坡角，有可能减少开挖土量。下文将采用数值极限分析模型的参数分析，对比这两种放坡的效果。

Figure 5. Comparison of safety factors of different slope angles with excavation volume

图5. 不同放坡坡角的安全系数与挖方量对比

图5中是图4(a)的一级放坡方案中随坡角不断减小边坡稳定安全系和土挖方量的变化情况。从图中可以看出随着挖方土量不断的增加，坡角不断减小后，边坡的安全系数不断增加，当坡角$\beta >26.6˚$ 时，安全系数逐渐增大超过1，此时挖方量则超过了3300 m³。不同的放坡方案有可能减少土方开挖量，为此，这里考虑图4(b)的二级放坡方案，建立数值极限分析模型进行不同放坡方案的参数对比分析。

Figure 6. Failure mechanisms of two slope cutting schemes when the safety factor is 1

图6. 安全系数为1时两种切坡方案的失稳模式

Figure 7. Comparison of excavation between two slope cutting options

图7. 两种放坡方案的挖方对比

图6对比了图4中两种放坡方案边坡整体安全系数等于1时的失稳模式。根据塑性剪切耗散自适应结果可以看出两种切坡方案的失稳模式略有不同，方案1的失稳模式明显符合对数螺线的特征，而方案2则接近两个对数螺线的复合模式。相对于方案1而言，方案2中分级放坡后，上部土体对下部土体而言成为了超载作用，导致下部失稳时的滑动面的形状和位置发生了变化，与上部土体中滑动面不一致。图7中对比了边坡安全系数为1时两种方案的挖方土量，方案1的挖方量约为6153 m³，方案2的挖方量约为4227 m³，其中阴影部分为两种放坡方案的挖方土量之差约为1926 m³。可见，两种放坡方案在安全系数相等的情况下，方案2中采用二级放坡的方式挖方量相对较少，可以降低边坡治理的工程成本。

Figure 8. Different failure mechanisms of two-stage bench excavation scheme

图8. 两级台阶挖方方案的不同失稳模式

然而二级放坡中边坡的失稳模式相对复杂，如图8所示，三种模型放坡的尺寸不同，两级边坡的坡度也不同，虽然三种模型的安全系数接近，但失稳模式却不同。图8(a)和图8(c)属于局部失稳，而图8(b)为复合对数螺线的整体失稳模式。此种现象可能是由于不同的放坡中开挖导致应力状态改变规律不同。图8(a)中下部坡面附近土体的应力状态较其它区域更接近屈服面，导致此区域先发生滑动。而图8(b)的这样的区域则出现在上部土体。由此可见，在方案2中虽然挖方土量减少了，但边坡的失稳模式变得相对复杂，对挖方施工和边坡治理中其它工序带来了一定的难度。

4. 结论

本文将OptumG2中自适应mixed-T6-FEM极限分析用于深圳龙岗区风门坳科技园后侧非均质边坡的稳定性分析，通过数值计算的参数分析，发现该边坡处于不稳定状态，失稳模式呈对数螺线特征。在此基础上，对比了不同的放坡方案的稳定性和挖方土量，在一级放坡方案中随着边坡坡度的降低，边坡稳定性不断增强。同时发现在安全系数相同的情况下，二级放坡方案的挖方量相对较少，但两级放坡后边坡的失稳模式相对复杂，可能出现局部失稳模式，也可能是两条对数螺旋线复合的失稳模式，可能会增加边坡治理的难度。

致 谢

作者感谢OptumCE提供OptumG2程序(学术版)的免费访问权限以进行此项研究。