摘要：
为研究三维边坡在对数螺旋破坏机制下的稳定性，基于极限分析上限法提出了非相关联流动法则下的速度场及安全系数的求解方法，可有效的解决相关联法则中摩擦角过大可能导致速度场不收敛的问题。基于Michalowski构建了考虑体积应变的均质三维边坡的对数螺旋线破坏机制，并利用MATLAB编制计算程序计算不同内摩擦角φ、不同边坡倾角β工况下边坡的稳定系数，分析内摩擦角、边坡倾角对稳定系数和破坏模式的影响规律，引入强度折减技术计算边坡安全系数，通过与Michalowski方法对比发现计算结果接近。研究结果表明：在三维边坡对数螺旋破坏机制下相对于一定宽度的边坡适用，且经过计算稳定系数符合规律。研究成果为边坡的设计计算提供参考。
关键词：
边坡工程; 极限分析; 三维边坡; 非相关流动法则;
作者简介：
董鸣镝(1970—),男，博士研究生，主要研究方向为工程项目管理。E-mail:1113602@163.com;
基金：
国家重点研发计划(2018YFC0407102);
引用：
董鸣镝． 基于对数螺旋破坏机制的三维边坡稳定性分析［J］. 水利水电技术( 中英文) ，2021，52( 3) : 169-176.
DONG Mingdi． Three dimensional slope stability analysis based on logarithmic spiral failure mechanism［J］. Water Ｒesources and Hydropower Engineering，2021，52( 3) : 169-176.
图1 垂直条分法滑动土体及条柱受力示意
0 引 言边坡稳定性研究较常采用的方法主要有三种：极限平衡法、极限分析法和有限元法。极限平衡法，即条分法，是最早被提出的边坡稳定性分析方法。1936年，极限平衡法被正式用于边坡稳定性分析，这一方法被命名为瑞典条分法。瑞典条分法假设边坡滑裂面为圆弧，在边坡体表面划分直线，使之成为一系列竖直条块，忽略条块间的相互作用力，对各条块进行静力分析，建立滑坡体的静力、力矩平衡方程，解方程求得边坡安全系数。而后极限平衡法研究分支不断发展，1954年，JANBU 将边坡滑裂面优化为非圆弧线，但该方法并不满足力矩平衡方程。1955年，BISHOP 研究时引入条块间的法向力，以土体抗剪强度与滑裂面剪应力的比值作为安全系数，将平衡方程不断迭代求得结果。随着计算机技术的飞速发展，又产生了MONGENSTERN-PRICE法 (1965)和SPENCER法 (1967)。此后，随着有限元理的建立和发展，ZIENKIEWICZ将边坡的稳定性分析带入到有限元时代，后面发展成为目前应用最广的数值分析法 。由于过程简明，维计算方法结果往往偏保守，误差甚至达15% ,而三维边坡稳定性分析能更加贴近真实边坡破坏形态。尽管如此，相对于二维计算方法，三维计算方法尚未成熟，国内外诸多学者进行了这方面的研究。
随着理论体系的成熟，能较好的满足工程要求，二维边坡稳定性计算方法目前被广泛采用 。但是二边坡三维稳定性分析方法主要有三维极限平衡法和三维极限分析法。与二维计算方法原理类似，三维极限平衡法将边坡划分为系列垂直条块，如图1所示，利用静力平衡方程求解边坡的安全系数 。
在三维极限平衡分析中，垂直条块间的力的大小、方向、作用点均未知，且单个垂直条块前后左右均有作用力，同二维边坡相比，未知变量个数远大于平衡方程个数。根据Lma&Fredlund, 由于求解三维边坡安全系数存在许多独立变量，需要引入大量假设，以有n行m列垂直条块的边坡为例，需要引入8 mn个假设 。也有其他学者采用不同的未知量假设。由于垂直条块数量众多，方程数量庞大，算法的收敛性是个需要解决的问题 ,另外搜索最优三维滑动面也是一个巨大的挑战。
三维边坡稳定性分析，除了极限平衡法，另外一种较常用的方法是极限分析法。1975年，GIGER将极限分析法引入到三维竖直切坡稳定性分析中去，开创了三维边坡极限分析的先例，但其研究方法仅适用于三维竖直切坡，并不能分析一般边坡。基于Mohr-Coulomb准则和相关联流动法则，MICHALOWSKI 将三维边坡划分为一系列滑块，相邻滑块的交界面垂直于滑体的对称轴，并采用二维边坡的方法构建速度场，以此分析三维边坡稳定性。FARZANEH 改进了上述方法，将其用于求解三维非均质边坡的稳定性问题，并讨论了用迭代算法求无约束和有约束问题的最小上限解。SUOBAR 用MICHALOWSKI的方法计算三维挡土墙的被动土压力，由于此方法将三维速度场简化成二维，其应用性有限。CHEN 引入中性面概念，将中性面划分为一系列倾斜条块，并采用一定方法对中性面进行空间扩展，将边坡离散成为一系列三维倾斜条块，应用Mohr-Coulomb准则和相关联流动法则构建速度场，最终求出三维边坡的安全系数，并在此基础上提出了一种应用单纯形自动搜索临界滑动面的方法。CHEN在构建三维速度场时引入两个假设，破坏了解的严谨性，但其进行的大量算例表明两个假设引起的误差不大。孙平 对上述方法中由于相关联法则中摩擦角过大可能导致速度场不收敛的问题进行了研究，并提出了非相关联流动法则下的速度场及安全系数的求解方法。MICHALOWSKI 构建了考虑体积应变的均质三维边坡的对数螺旋线破坏机制，并提出相应研究思路。
1 基于对数螺旋破坏机制1975年，CHEN 在其著作“Limit Analysis and Soil Plasticity”中提出一种二维边坡旋转破坏机制。该机构将二维边坡视作刚体，假设其滑裂面为对数螺旋线，利用极限分析上限法推导出稳定系数的解析解，并编制计算程序求得其稳定系数。此后，对于二维均质边坡，诸多学者采用对数螺旋线破坏模式进行稳定性分析。但实际工程中，边坡都呈现三维形态，采用二维破坏模式低估了实际边坡的复杂性，且采用二维破坏模式所得安全系数小于实际情况,影响实际加固和设计方案的选取。
为解决上述问题，众多学者对三维边坡破坏模式进行了研究，其中Michalowski引入对数螺旋线圆锥破坏模式，对二维对数螺旋线模式进行空间扩展，并考虑边坡体的体积应变，推导三维边坡稳定系数的解析公式，编制相应计算程序。本文将对此破坏模式进行详细介绍。
1.1 三维边坡极限分析理论
极限分析上限定理认为，在一个假设的满足速度边界条件和应变与速度相容条件的速度场中，由外荷载功率等于所消耗的内力功率所确定的荷载一定不小于实际破坏荷载。满足上述条件的速度场叫做运动许可速度场。因此，如果能找到一个运动许可速度场，则土体内的塑性流动必将发生或早已发生。根据虚功率原理推导出上限定理为：在所有的机动容许的塑性变形速度场相对应的荷载中，极限荷载为最小。即
式中，σi*j 为由εi*j 按塑性变形法则求出的应力；T 、F 分别为作用物体上的面力和体力；Δv为速度间断线两侧切向速度变化量；u 、εi*j 分别为速度场中的速度和应变率；c、φ为抗剪强度指标中的粘聚力和摩擦角；dA, dV, dL为积分变量，分别对应积分面元，积分体积元和积分线元。
上式所确定的力系T 及F 大于或等于真实的极限荷载。也就是说由上述方程确定的力系T 、F 是真实极限荷载的一个上限解，我们最终需要通过程序优化出来一个最小上限解。按照极限分析理论，能量耗散发生在土体塑性区内。然而，为了简化问题，对边坡进行稳定性分析时，经常假设边坡体由系列刚性块体组成，此时能量耗散只发生在速度间断线上，而速度间断线是各相邻刚性块体之间的窄过渡层。本文仍然采用刚性块体假设，不考虑边坡的体积应变。
1.2 破坏机制的构建
如图2所构建的稳定性分析机构示意图，假设一均质边坡高度为H,宽度为B,倾角为β,土体黏聚力为c,内摩擦角为φ,重度为γ。当边坡失稳时，以角速度ω围绕旋转中心O旋转，曲线AC是式(1)确定的对数螺旋线，为坡体的滑裂面，该曲线通过坡趾C。数螺旋线A′C′确定了三维边坡滑坡体的上界，曲线AC与曲线A′C′相交于一点。以O为原点向边坡体做射线OA,交曲线A′C′于A′,以AA′为直径沿垂直于纸面方向做圆，该圆以O为旋转中心进行旋转，直径渐变，为上述射线与曲线A′C′和曲线AC两交点确定的线段。圆围绕O点旋转与边坡体相交形成的曲面为三维边坡的滑裂面，整个旋转机构为螺旋锥体，按照相关联流动法则，其顶角夹角为2φ,整个滑坡体的上边界A′C′和下边界AC由式(1)和式(2)确定，至此，边坡的三维旋转机构构建完成
式中，ro和r′o分别为图2所示稳定性分析机构中线段OA、OA′的长度；θo为线段OA与水平面的夹角；θ滑面AC上任意一点于机构旋转中心连线于水平线所成夹角；φ为边坡土体的摩擦角。
图2 三维边坡旋转破坏机制
滑面AC上任意一点于机构旋转中心连线于水平线所成夹角记为θ,每一个θ都有一个其对应的旋转断面的圆(如图2中阴影部分的圆所示),则其半径为用R表示(圆心点与旋转中心O的连线的长度),可由下式计算
所有旋转断面的圆心所构成的曲线的极坐标方程可表示为
以每个断面圆圆心为原点，做如上图2所示局部坐标系，即以圆心点与旋转中心的连线为y轴，y轴的垂直方向为x轴，则边坡体内任意一点的速度为
式中，vi为滑体中任意一点在绕中心O旋转过程中的速度；rmi为任意一点所对应断面圆圆心旋转中心的连线的长度；yi为任意一点在其所在局部坐标系中的纵坐标值；ω为滑体绕O点旋转的的角速度。该模型可用于有宽度限制的边坡稳定系数的计算。对于无宽度限制的边坡，该破坏机制不太适用 。
如图3所示，旋转体体积微元可由式(7)计算得到，其中dV表示体积微元的大小，r 微元体所对应断面圆圆心旋转中心的连线的长度，y表示微元体在其所在局部坐标系中的纵坐标值。
图3 旋转体体积微元
则根据力做功功率的计算公式P=F·V,可求得滑坡体微元重力功率可表示为
式中，γ为边坡土体的容重；vi为微元体运动的速度其数值，由公式(6)计算可得；θ为体积微元与机构旋转中心连线O于水平线所成夹角。
通过滑坡体微元重力做功进行积分，求得边坡滑坡体重力功率可表示为
式中，参数θc, θb, θhx1* ,x2* , y , a, d为积分计算公式中的积分边界，其数值可由下式计算求得
在求解公式(9)时，先对y进行积分，由图3可知，y的积分范围在区间(θ ,θ )内为(a,y ),在区间(θ ,θ )内为(d,y ),其中x1* 、x2* 分别表示坡顶、坡面某一旋转断面的最大宽度。
由于不考虑体积应变，内能耗散仅发生在速度间断面上，不考虑边坡的体积应变，内能耗散发生在空间曲面AC上，Chen 在其著作“Limit Analysis and Soil Plasticity”中指出内部耗损率可以由该面的微分面积dS与粘聚力c以及与跨该面的切向间断速度(vcosφ)的连乘积沿整个间断面积分，即得到总的内部能量耗损率，其表达式为
其中参数θo, θB,θh为同上所述，积分变量dα为积微元体在局部坐标系中与原点连线同x轴之间的夹角，α 和α为积分变量dα对应的积分边界，其数值可由公式(18)和公式(19)计算公式求得，其中参数α 、α 可由公式(10)和公式(11)计算求得
2.1 稳定系数计算
根据极限分析上限定理，外力做功功率与内能耗散功率相等即
整理式(20)后可知，稳定系数γH/c可表示成θo, θh, r′0/r0的函数，求边坡的最小稳定系数γH/c本质上是非线性规划问题，使用MATLAB编制计算程序，能方便的求解稳定系数。若给出土体的c,γ,还可通过优化程序求得边坡的临界高度。
如图4所示，旋转体内插入一宽度为b的块体，当插入块的宽度b不断增大，三维边坡破坏机制宽度将不断增大。当b→∞时，有限宽度边坡向无限宽度边坡转换，最终变为二维平面应变模式。插入块的重力功率和内能耗散功率可由二维情况的重力功率和内能耗散功率乘以b获得。
图4 三维边坡插入块体示意
求解边坡的稳定系数公式比较复杂，本文采用MATLAB软件编制程序进行计算。求解边坡的最小稳定系数属于非线性约束二次规划问题，MATLAB的fmincon函数是求解含约束优化问题的比较高效的函数之一。本文采用fmincon函数进行计算。该函数在给定的约束条件内对稳定系数γH/c进行优化，最终求得γH/c的最小值。除了用fmincon函数求解最小稳定系数，还可以采用枚举法，先根据约束条件划分比较粗糙的求解区间，并搜索该区间内的最小稳定系数，确定最小安全系数所在的目标区间后，细化该区间进行第二次循环，便可找到较优的稳定系数。枚举法虽然不能非常精确的求得最小稳定系数，经过比较，其精度完全可以满足计算要求，但是由于枚举法需要计算每个步长对应的解，即使增大步长，工作量依然巨大，导致计算时间较长，故本文采用fmincon函数求解稳定系数。
2.1.1 稳定系数对比
为验证结果的正确性，需要与二维边坡稳定系数 进行对比，结果如表1所列。
表1 普通均质边坡二维与三维极限分析稳定系数
由表1可知，当B/H=10时，三维边坡稳定系数略大于二维情况稳定系数，这是由于当b不断增大三维边坡逐渐向二维边坡模式转变，三维稳定系数不断趋近于二维稳定系数，同时也可表明，二维平面应变破坏机制比三维破坏机制计算结果危险 。
另将本文计算结果与Michalowski计算结果对比，列于表2至表4中。可以发现，本文计算结果与Michalowski计算结果相差无几，且绝大部分小于后者，而极限分析上限法所求为上限解，其值越小越具有意义，从与二维计算结果及Michalowski计算结果对比可知，本文计算可行。
2.1.2 规律分析
图4展示了不同内摩擦角及不同边坡倾角下γH/c的变化规律。以图4(a)为例，当内摩擦角φ和边坡倾角β均为定值时，边坡的稳定系数γH/c将随着边坡宽高比B/H的增大不断减小，当B/H<5时，γH/c随着B/H的增大急剧减小，当B/H>5时，γH/c减小幅度比较有限，将趋于稳定，这是由于随着B/H的增大，三维模型逐渐趋近于二维平面应变模型。
当内摩擦角φ和边坡宽高比B/H为定值时，γH/c随着边坡倾角β的增大而减小，这表明当边坡宽度一定时，随着边坡倾角的增大，边坡失稳的临界高度不断减小，即边坡倾角越大越容易失稳。
表2 本文方法与Michalowski三维极限分析稳定系数对比分析(φ=0°)
表3 本文与Michalowski三维极限分析稳定系数对比分析(φ=15°)
表4 本文与Michalowski三维极限分析稳定系数对比分析(φ=30°)
当宽高比B/H和边坡倾角β为定值时，γH/c随着摩擦角φ的增大而不断增大。这表明当边坡倾角、边坡宽度不变时，边坡失稳的临界高度会随着内摩擦角的增大而增大，即土体内摩擦角越大越不容易失稳。
图5 不同工况下的稳定系数
2.2 安全系数计算
边坡稳定性的评定指标有很多种，如临界高度Hcr、稳定系数N =γH/c,安全系数F ,而工程中最常用的边坡稳定性评价指标为安全系数。现实边坡往往具有一定的安全储备，通过折减土体强度参数，把潜在滑动土体的抗剪强度参数折减Fs倍，边坡恰好过渡到临界极限平衡状态。对滑动土体的抗剪强度进行折减，实际上就是对滑面土体的强度参数进行折减使滑面上的内部耗能减小并使其达到与外力做功相等的临界状态。折减的倍数反映了边坡的安全储备，称之为安全系数Fs。自强度折减技术与极限分析法结合以来，众多学者采用安全系数对二维边坡进行了稳定性分析，并取得了非常丰富的成果，然而，采用安全系数Fs分析三维边坡稳定的研究相对较少。本节将推导三维边坡安全系数公式，利用MATLAB编制求解公式，利用安全系数这一评价指标对三维边坡的稳定性进行分析
式中，Fs为剪切强度折减系数；c和φ为原始抗剪强度参数；cf和φf为折减后的抗剪强度参数。
采用这一方法时，Fs以隐式出现在求解的方程式中，计算时常需要进行迭代。将式(21)折减后的强度参数cf和φf代替式中(20)的强度参数c和φ进行迭代求解，可得Fs是关于位置参数θo, θh, r′0, r0的隐函数
为验证该三维边坡安全系数计算结果的正确性，选取典型边坡计算，对比FLAC3D计算结果，列于表5。本文三维计算案例中选取固定b/H值，由表5对比结果可知，本文计算结果与有限差分软件FLAC3D计算结果相差不大，证明了本文计算方法的有效性。
表5 普通边坡FLAC3D与三维极限分析安全系数对比
3 结 论本章基于极限分析上限法，构建了三维边坡破坏机制，推导了稳定系数、安全系数解析公式。通过计算不同工况下边坡的稳定系数、安全系数，主要得到如下结论：
(1)三维边坡对数螺旋破坏机制中，边坡倾角、内摩擦角一定时，稳定系数随着边坡宽高比B/H的增大先急剧减小，后趋于稳定；内摩擦角、宽高比一定时，稳定系数随着边坡倾角的增大而减小；宽高比、边坡倾角一定时，稳定系数随着内摩擦角的增大而增大。
(2)当插入块的宽度b不断增大时，稳定系数不断减小，最终趋近于二维平面应变稳定系数，但仍然比二维结果大，在表明该破坏机制中的边坡由三维状态逐渐向二维平面应变状态转变的同时，也说明同等条件下二维计算结果更加危险。
(3)同等工况下，三维边坡安全系数与FLAC3D计算结果相差无几，表明引入安全系数有效。
水利水电技术
水利部《水利水电技术》杂志是中国水利水电行业的综合性技术期刊（月刊），为全国中文核心期刊，面向国内外公开发行。本刊以介绍我国水资源的开发、利用、治理、配置、节约和保护，以及水利水电工程的勘测、设计、施工、运行管理和科学研究等方面的技术经验为主，同时也报道国外的先进技术。期刊主要栏目有:水文水资源、水工建筑、工程施工、工程基础、水力学、机电技术、泥沙研究、水环境与水生态、运行管理、试验研究、工程地质、金属结构、水利经济、水利规划、防汛抗旱、建设管理、新能源、城市水利、农村水利、水土保持、水库移民、水利现代化、国际水利等。
