水利水电技术(中英文)㊀第52卷㊀2021年第3期
董鸣镝.基于对数螺旋破坏机制的三维边坡稳定性分析[J].水利水电技术(中英文),2021,52(3):169-176.
DONG Mingdi.Three dimensional slope stability analysis bad on logarithmic spiral failure mechanism[J].Water Resources and Hydropow-er Engineering,2021,52(3):169-176.
基于对数螺旋破坏机制的三维边坡稳定性分析
董鸣镝
(中国矿业大学(北京)，北京㊀100083)
收稿日期:2020-11-09
基金项目：国家重点研发计划(2018YFC0407102)
作者简介：董鸣镝(1970 )，男，博士研究生，主要研究方向为工程项目管理㊂E-mail: 摘㊀要：为研究三维边坡在对数螺旋破坏机制下的稳定性，基于极限分析上限法提出了非相关联流动法则下的速度场及安全系数的求解方法，可有效的解决相关联法则中摩擦角过大可能导致速度场不收敛的问题㊂基于Michalowski 构建了考虑体积应变的均质三维边坡的对数螺旋线破坏机制，并利用MATLAB 编制计算程序计算不同内摩擦角φ㊁不同边坡倾角β工况下边坡的稳定系数，分析内摩擦
角㊁边坡倾角对稳定系数和破坏模式的影响规律，引入强度折减技术计算边坡安全系数，通过与Michalowski 方法对比发现计算结果接近㊂研究结果表明：在三维边坡对数螺旋破坏机制下相对于一
定宽度的边坡适用，且经过计算稳定系数符合规律㊂研究成果为边坡的设计计算提供参考
㊂关键词：边坡工程；极限分析；三维边坡；非相关流动法则
doi :10.13928/jki.wrahe.2021.03.019开放科学(资源服务)标志码(OSID )：中图分类号:TU43
文献标志码:A
文章编号:1000-0860(2021)03-0169-08
Three dimensional slope stability analysis bad on logarithmic spiral failure mechanism云雾山风景区
DONG Mingdi
(China University of Mining and Technology-Beijing,Beijing㊀100083,China)
Abstract :In order to study the stability of three-dimensional slope under logarithmic spiral failure mechanism,bad on the
upper limit method of limit analysis,the solution method of velocity field and safety factor under the non-associated flow rule is京华时报
put forward,which can effectively solve the problem that the velocity field may not converge due to excessive friction angle in
the associated rule.Bad on Michalowski,the logarithmic spiral failure mechanism of homogeneous three-dimensional slope
considering volume strain is constructed,and the calculation program is compiled by MATLAB to calculate the slope stability
coefficient under different internal friction angles φand slope inclination angles β.The influence of inte
rnal friction angles and slope inclination angles on stability coefficient and failure mode is analyzed,and the strength reduction technology is introduced
to calculate the slope safety factor.Compared with Michalowski method,the calculation results are clo to each other.The
results show that the logarithmic spiral failure mechanism of three-dimensional slope is applicable to slope with a certain width,and the stability coefficient is in accordance with the law after calculation.The rearch results provide reference for slope
design and calculation.
Keywords :slope engineering;limit analysis;three-dimensional slope;irrelevant flow rule
董鸣镝//基于对数螺旋破坏机制的三维边坡稳定性分析
图1㊀垂直条分法滑动土体及条柱受力示意
Fig.1㊀Stress diagram of sliding soil and column by vertical slice method
0㊀引㊀言
㊀㊀边坡稳定性研究较常采用的方法主要有三种：极
限平衡法㊁极限分析法和有限元法㊂极限平衡法，即条分法，是最早被提出的边坡稳定性分析方法㊂
1936年，极限平衡法被正式用于边坡稳定性分析，这一方法被命名为瑞典条分法㊂瑞典条分法假设边坡滑裂面为圆弧，在边坡体表面划分直线，使之成为一系列竖直条块，忽略条块间的相互作用力，对各条块进行静力分析，建立滑坡体的静力㊁力矩平衡方程，解方程求得边坡安全系数㊂而后极限平衡法研究分支不断发展,1954年,JANBU [1]
将边坡滑裂面优化为非圆弧
线，但该方法并不满足力矩平衡方程㊂1955年,
BISHOP
[2]
研究时引入条块间的法向力，以土体抗剪
强度与滑裂面剪应力的比值作为安全系数，将平衡方程不断迭代求得结果㊂随着计算机技术的飞速发展，又产生了MONGENSTERN-PRICE 法[3]
(1965)和
SPENCER 法
[4]
(1967)㊂此后，随着有限元理的建立
和发展,ZIENKIEWICZ 将边坡的稳定性分析带入到有限元时代，后面发展成为目前应用最广的数值分析法
[5]
㊂由于过程简明，维计算方法结果往往偏保守,
误差甚至达15%[5]，而三维边坡稳定性分析能更加
贴近真实边坡破坏形态㊂尽管如此，相对于二维计算方法，三维计算方法尚未成熟，国内外诸多学者进行了这方面的研究㊂
随着理论体系的成熟，能较好的满足工程要求，二维边坡稳定性计算方法目前被广泛采用
[7-9]
㊂但是
二边坡三维稳定性分析方法主要有三维极限平衡法和三维极限分析法㊂与二维计算方法原理类似，三维极限平衡法将边坡划分为系列垂直条块，如图1所示，利用静力平衡方程求解边坡的安全系数
[6-8]
㊂
在三维极限平衡分析中，垂直条块间的力的大
小㊁方向㊁作用点均未知，且单个垂直条块前后左右均有作用力，同二维边坡相比，未知变量个数远大于平衡方程个数㊂根据Lma&Fredlund，由于求解三维边坡安全系数存在许多独立变量，需要引入大量假设，以有n 行m 列垂直条块的边坡为例，需要引入8mn 个假设[9]㊂也有其他学者采用不同的未知量假
设㊂由于垂直条块数量众多，方程数量庞大，算法的收敛性是个需要解决的问题[10-13]，另外搜索最优三
维滑动面也是一个巨大的挑战㊂
三维边坡稳定性分析，除了极限平衡法，另外一种较常用的方法是极限分析法㊂1975年,GIGER 将
极限分析法引入到三维竖直切坡稳定性分析中去，开创了三维边坡极限分析的先例，但其研究方法仅适用于三维竖直切坡，并不能分析一般边坡㊂基于Mohr-Coulomb 准则和相关联流动法则,MICHALOWSKI [14]
将三维边坡划分为一系列滑块，相邻滑块的交界面垂直于滑体的对称轴，并采用二维边坡的方法构建速度场，以此分析三维边坡稳定性㊂FARZANEH
[15]
改进
了上述方法，将其用于求解三维非均质边坡的稳定性问题，并讨论了用迭代算法求无约束和有约束问题的最小上限解㊂SUOBAR [16]用MICHALOWSKI 的方法计算三维挡土墙的被动土压力，由于此方法将三维速度场简化成二维，其应用性有限㊂CHEN [17]引入中性面概念，将中性面划分为一系列倾斜条块，并采用一定方法对中性面进行空间扩展，将边坡离散成为一系列三维倾斜条块，应用Mohr-Coulomb 准则和相关联流动法则构建速度场，最终求出三维边坡的安全系数，并在此基础上提出了一种应用单纯形自动搜索临界滑动面的方法㊂CHEN 在构建三维速度场时引入两个假设，破坏了解的严谨性，但其进行的大量算例表明两个假设引起的误差不大㊂孙平[18]对上述方法中由于相关联法则中摩擦角过大可能导致速度场不收敛
的问题进行了研究，并提出了非相关联流动法则下的速度场及安全系数的求解方法㊂MICHA-
LOWSKI [19]构建了考虑体积应变的均质三维边坡的对数螺旋线破坏机制，并提出相应研究思路㊂1㊀基于对数螺旋破坏机制
㊀㊀1975年,CHEN [20]在其著作
Limit Analysis and Soil Plasticity
中提出一种二维边坡旋转破坏机制㊂该机构将二维边坡视作刚体,
董鸣镝//基于对数螺旋破坏机制的三维边坡稳定性分析
假设其滑裂面为对数螺旋线，利用极限分析上限法推
导出稳定系数的解析解，并编制计算程序求得其稳定
系数㊂此后，对于二维均质边坡，诸多学者采用对数
螺旋线破坏模式进行稳定性分析㊂但实际工程中，边
坡都呈现三维形态，采用二维破坏模式低估了实际边
坡的复杂性，且采用二维破坏模式所得安全系数小于
实际情况[21-24]，影响实际加固和设计方案的选取㊂为解决上述问题，众多学者对三维边坡破坏模式
进行了研究，其中Michalowski引入对数螺旋线圆锥
破坏模式，对二维对数螺旋线模式进行空间扩展，并
考虑边坡体的体积应变，推导三维边坡稳定系数的解
析公式，编制相应计算程序㊂本文将对此破坏模式进
行详细介绍㊂
1.1㊀三维边坡极限分析理论
㊀㊀极限分析上限定理认为，在一个假设的满足速度边界条件和应变与速度相容条件的速度场中，由外荷载功率等于所消耗的内力功率所确定的荷载一定不小于实际破坏荷载㊂满足上述条件的速度场叫做运动许可速度场㊂因此，如果能找到一个运动许可速度场，则土体内的塑性流动必将发生或早已发生㊂根据虚功率原理推导出上限定理为：在所有的机动容许的塑性
变形速度场相对应的荷载中，极限荷载为最小㊂即ʏA T i u i d A+ʏΩF i u i d VȡʏΩσ∗ijε∗ij d V+ʏl c㊃Δv㊃d L
(1)式中,σ∗ij为由ε∗ij按塑性变形法则求出的应力;T i㊁F i分别为作用物体上的面力和体力;Δv为速度间断线两侧切向速度变化量;u i㊁ε∗ij分别为速度场中的速度和应变率;c㊁φ为抗剪强度指标中的粘聚力和摩擦角;d A,d V,d L为积分变量，分别对应积分面元，积分体积元和积分线元㊂
上式所确定的力系T i及F i大于或等于真实的极
限荷载㊂也就是说由上述方程确定的力系T i㊁F i是真实极限荷载的一个上限解，我们最终需要通过程序优化出来一个最小上限解㊂按照极限分析理论，能量耗散发生在土体塑性区内㊂然而，为了简化问题，对边坡进行稳定性分析时，经常假设边坡体由系列刚性块体组成，此时能量耗散只发生在速度间断线上，而速度间断线是各相邻刚性块体之间的窄过渡层㊂本文仍然采用刚性块体假设，不考虑边坡的体积应变㊂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)确定，至此，边坡的三维旋转机构构建完成
r=r0e(θ-θ0)tanφ(2)
rᶄ=rᶄ0e-(θ-θ0)tanφ(3)式中,r0和rᶄ0分别为图2所示稳定性分析机构中线段OA㊁OAᶄ的长度;θ0为线段OA与水平面的夹角;θ滑面AC上任意一点于机构旋转中心连线于水平线所成夹角;φ为边坡土体的摩擦角
㊂
图2㊀三维边坡旋转破坏机制
Fig.2㊀Mechanism of rotational failure of3D slope
滑面AC上任意一点于机构旋转中心连线于水平线所成夹角记为θ，每一个θ都有一个其对应的旋转断面的圆(如图2中阴影部分的圆所示)，则其半径为用R表示(圆心点与旋转中心O的连线的长度)，可由下式计算
董鸣镝//基于对数螺旋破坏机制的三维边坡稳定性分析
R=r-rᶄ2(4)㊀㊀所有旋转断面的圆心所构成的曲线的极坐标方程可表示为
r m=r+rᶄ2(5)㊀㊀以每个断面圆圆心为原点，做如上图2所示局部坐标系，即以圆心点与旋转中心的连线为y轴,y轴的垂直方向为x轴，则边坡体内任意一点的速度为
v i=(r mi+y i)ω(6)式中,v i为滑体中任意一点在绕中心O旋转过程中的速度;r mi为任意一点所对应断面圆圆心旋转中心的连线的长度;y i为任意一点在其所在局部坐标系中的纵坐标值;ω为滑体绕O点旋转的的角速度㊂该模型可用于有宽度限制的边坡稳定系数的计算㊂对于无宽度限制的边坡，该破坏机制不太
适用[25-26]㊂如图3所示，旋转体体积微元可由式(7)计算得到，其中dV表示体积微元的大小,r m微元体所对应断面圆圆心旋转中心的连线的长度,y表示微元体在其所在局部坐标系中的纵坐标值㊂
d V=(r m+y)dθd x d y
服装营销策略(7)
图3㊀旋转体体积微元
Fig.3㊀Volume element of revolution
㊀㊀则根据力做功功率的计算公式P=F㊃V，可求得滑坡体微元重力功率可表示为
d Wγ=γʏV v i cosθd V(8)式中,γ为边坡土体的容重;v i为微元体运动的速度其数值，由公式(6)计算可得;θ为体积微元与机构旋转中心连线O于水平线所成夹角㊂
通过滑坡体微元重力做功进行积分，求得边坡滑坡体重力功率可表示为
Wγ=2ωγ[ʏθBθ0ʏx1∗0ʏy∗a(r m+y)2cosθd y d x dθ+ʏθhθBʏx2∗0ʏy∗d(r m+y)2cosθd y d x dθ](9)
式中，参数θ0,θB,θh x∗1,x∗2,y∗,a,d为积分计算公式中的积分边界，其数值可由下式计算求得
d=
sin(β+θh)
sin(β+θ)r0e(θh-θ0)tanφ-r m(10)
a=
sinθ0
sinθr0-r m(11)
y∗=R2-x2(12)
x∗1=R2-a2(13)
x∗2=R2-d2(14)
θB=arctan
sinθ0
cosθ0-A(15) A=
sin(θh-θ0)
sinθh-
e(θh-θ0)tanφsinθh-sinθ0
sinθh sinβsin(θh+β)
(16)㊀㊀在求解公式(9)时，先对y进行积分，由图3可知,y的积分范围在区间(θ0,θB)内为(a,y∗)，在区间(θB,θh)内为(d,y∗)，其中x∗1㊁x∗2分别表示坡顶㊁坡面某一旋转断面的最大宽度㊂
由于不考虑体积应变，内能耗散仅发生在速度间断面上，不考虑边坡的体积应变，内能耗散发生在空间曲面AC上,Chen[20]在其著作 Limit Analysis and Soil Plasticity 中指出内部耗损率可以由该面的微分面积d S与粘聚力c以及与跨该面的切向间断速度(v cosφ)的连乘积沿整个间断面积分，即得到总的内部能量耗损率，其表达式为
D=ʏS c(v cosφ)d S=2ωc[ʏθBθ0ʏα∗10R(r m+ R cosα)2dαdθ+ʏθhθBʏα∗20R(r m+R cosα)2dαdθ](17)㊀㊀其中参数θ0,θB,θh为同上所述，积分变量dα为积微元体在局部坐标系中与原点连线同x轴之间的夹角,α1∗和α2∗为积分变量dα对应的积分边界，其数值可由公式(18)和公式(19)计算公式求得，其中参数α∗1㊁α∗2可由公式(10)和公式(11)计算求得
α∗1=arccos a R()(18)
权倾官海
α∗2=arccos d R()(19)
董鸣镝//基于对数螺旋破坏机制的三维边坡稳定性分析
2㊀稳定性与算例分析
2.1㊀稳定系数计算㊀㊀根据极限分析上限定理，外力做功功率与内能耗
散功率相等即
W γ=D
(20)
㊀㊀整理式(20)后可知，稳定系数γH /c 可表示成
θ0,θh ,rᶄ0/r 0的函数，求边坡的最小稳定系数γH /c 本质上是非线性规划问题，使用MATLAB 编制计算程序，能方便的求解稳定系数㊂若给出土体的c ,γ,
还可通过优化程序求得边坡的临界高度㊂
如图4所示，旋转体内插入一宽度为b 的块体，当插入块的宽度b 不断增大，三维边坡破坏机制宽度将不断增大㊂当b ңɕ时，有限宽度边坡向无限宽度边坡转换，最终变为二维平面应变模式㊂插入块的重力功率和内能耗散功率可由二维情况的重力功率和内能耗散功率乘以b 获得
什么人不用电
㊂
图4㊀三维边坡插入块体示意
Fig.4㊀Block diagram of 3D slope inrtion西湖龙井属于什么茶
求解边坡的稳定系数公式比较复杂，本文采用
MATLAB 软件编制程序进行计算㊂求解边坡的最小稳
定系数属于非线性约束二次规划问题,MATLAB 的
fmincon 函数是求解含约束优化问题的比较高效的函数之一㊂本文采用fmincon 函数进行计算㊂该函数在给定的约束条件内对稳定系数γH /c 进行优化，最终求得γH /c 的最小值㊂除了用fmincon 函数求解最小稳定系数，还可以采用枚举法，先根据约束条件划分比较粗糙的求解区间，并搜索该区间内的最小稳定系数，确定最小安全系数所在的目标区间后，细化该区间进行第二次循环，便可找到较优的稳定系数㊂枚举法虽然不能非常精确的求得最小稳定系数，经过比较，其精度完全可以满足计算要求，但是由于枚举法需要计算每个步长对应的解，即使增大步长，工作量依然巨大，导致计算时间较长，故本文采用fmincon 函数求解稳定系数㊂
2.1.1㊀稳定系数对比
为验证结果的正确性，需要与二维边坡稳定系数[55]进行对比，结果如表1所列㊂
表1㊀普通均质边坡二维与三维极限分析稳定系数Table 1㊀Stability coefficient of two-dimensional and
three-dimensional limit analysis of ordinary
homogeneous slope
内摩擦角
φ/(ʎ)
坡脚β/(ʎ)
稳定系数γH /c
二维极限分析
三维极限分析(刚体,B /H =10)
15
30
30
安徒生童话阅读
6.51 6.8045 5.86 6.1560 5.25 5.4875 4.57 4.7890 3.83
4.03
3021.7122.604512.0512.49608.648.9475 6.57 6.8390 5.02 5.25--4535.6336.556016.1116.51759.9610.2890
6.69 6.97㊀㊀由表1可知，当B /H =10时，三维边坡稳定系
数略大于二维情况稳定系数，这是由于当b 不断增大三维边坡逐渐向二维边坡模式转变，三维稳定系数不断趋近于二维稳定系数，同时也可表明，二维平面应变破坏机制比三维破坏机制计算结果危险[27-28]㊂
另将本文计算结果与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 随着边坡倾角β的增大而减小，这表明当边坡宽度一定时，随着边坡倾角的增大，边坡失稳的临界
