《工程科学学报》录用稿，htps:/doi.org/10.13374/i,issn2095-9389.2021.04.29.002©北京科技大学2020 基于有限元极限平衡法的三维边坡稳定分析 苏振宁”，邵龙潭)回 1)大连理工大学，工业装备结构分析国家重点实验室，大连116024 ☒通信作者，E-mail:shaolt@dlut.edu.cn 摘要本文提出了一种基于有限元弹塑性应力场和极限平衡状态的三维边坡稳定分析方法一一三维有限元极限平衡法。 首先，考虑三维空间中滑动方向，提出滑动面上一点在滑动方向上的极限平衡条件，并证明滑动面上土体整体达到极限平 衡状态与滑动面上土体各处在滑动方向上处于极限平衡状态等价。再通过刚体极限平衡假定计算主滑方向和滑动面上各点 滑动方向。最后，定义局部安全系数为抗剪强度与滑动方向上剪应力投影的比值，基于三维边坡整体极限平衡条件将局部 安全系数通过积分中值定理转变为整体安全系数。该方法计算简单，消除了剪应力比形式定全系数滑动面形状限制， 具备合理性与有效性。算例验证结果表明，该方法滑动方向假设合理，安全系数与严格三维极限乘 衡法结果一致。 关键词三维边坡稳定：有限元：极限平衡条件：安全系数：滑动方向 分类号TU473 finite element limit 0o3sad,。A Three-dimensional slope stability analysis equilibrium method SU Zhenning,SHAO Longtan Technology,Dalian,116024,China Corresponding author,E-mail:shaolt@dlut.edu.cn Abstract The finite element limit equilibrium method is a slope stability analysis method based on finite element stress analysis combined with limit equilibrium condition.The safety factor defined in the form of shear stress ratio in the three-dimensional space does not consider the influence of the sliding direction on the calculation results.In this paper,a three-dimensional finite element limit equilibrium method is proposed by considering the sliding direction.This method is different from the limit equilibrium method and the strength reduction method. This method analyzes the slope stability through the "true"stress state without reducing the material strength parameters.First,considering the sliding direction in the three-dimensional space,the limit equilibrium condition of a point on the slip surface in the sliding direction is proposed.An equivalent relationship is proved that the slip surface as a whole is in the limit equilibrium state and each points of the slip surface are in the limit equilibrium state in the sliding direction Then,the main sliding direction and the sliding direction of each point on the slip surface are calculated by the assumption of rigid body limit equilibrium.Finally,the local safety factor is defined as the ratio of the shear strength to the shear stress projection in the sliding direction.Based on the equivalent relationship of limitequilibrium condition of the three-dimensional slope,the local safety factor is transformed into the global safety factor through the integral median theorem.The method is simple to calculate,eliminates the limitation of the slip surface shape of the safety factor defined by the shear stress ratio form,and is reasonable and effective.The verification result of the calculation example shows that the sliding direction assumption of the method is reasonable,and the safety factor is consistent with the result of the strict three-dimensional limit equilibrium method. Key words three-dimensional slope stability;finite element method;limit equilibrium condition;safety factor;sliding direction 收离日期： 基金项目：国家自然科学基金资助(No.52079018)工业装备结构分析国家重点实验室自主课题基金(No.S14208)
基于有限元极限平衡法的三维边坡稳定分析 苏振宁 1)，邵龙潭 1) 1 1) 大连理工大学，工业装备结构分析国家重点实验室，大连 116024 通信作者，E-mail: shaolt@dlut.edu.cn 摘 要 本文提出了一种基于有限元弹塑性应力场和极限平衡状态的三维边坡稳定分析方法——三维有限元极限平衡法。 首先，考虑三维空间中滑动方向，提出滑动面上一点在滑动方向上的极限平衡条件，并证明滑动面上土体整体达到极限平 衡状态与滑动面上土体各处在滑动方向上处于极限平衡状态等价。再通过刚体极限平衡假定计算主滑方向和滑动面上各点 滑动方向。最后，定义局部安全系数为抗剪强度与滑动方向上剪应力投影的比值，基于三维边坡整体极限平衡条件将局部 安全系数通过积分中值定理转变为整体安全系数。该方法计算简单，消除了剪应力比形式定义安全系数滑动面形状限制， 具备合理性与有效性。算例验证结果表明，该方法滑动方向假设合理，安全系数与严格三维极限平衡法结果一致。 关键词 三维边坡稳定；有限元；极限平衡条件；安全系数；滑动方向 分类号 TU473 Three-dimensional slope stability analysis based on finite element limit equilibrium method SU Zhenning1), SHAO Longtan1) 1) State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian, 116024, China Corresponding author, E-mail: shaolt@dlut.edu.cn Abstract The finite element limit equilibrium method is a slope stability analysis method based on finite element stress analysis combined with limit equilibrium condition. The safety factor defined in the form of shear stress ratio in the three-dimensional space does not consider the influence of the sliding direction on the calculation results. In this paper, a three-dimensional finite element limit equilibrium method is proposed by considering the sliding direction. This method is different from the limit equilibrium method and the strength reduction method. This method analyzes the slope stability through the “true” stress state without reducing the material strength parameters. First, considering the sliding direction in the three-dimensional space, the limit equilibrium condition of a point on the slip surface in the sliding direction is proposed. An equivalent relationship is proved that the slip surface as a whole is in the limit equilibrium state and each points of the slip surface are in the limit equilibrium state in the sliding direction. Then, the main sliding direction and the sliding direction of each point on the slip surface are calculated by the assumption of rigid body limit equilibrium. Finally, the local safety factor is defined as the ratio of the shear strength to the shear stress projection in the sliding direction. Based on the equivalent relationship of limit equilibrium condition of the three-dimensional slope, the local safety factor is transformed into the global safety factor through the integral median theorem. The method is simple to calculate, eliminates the limitation of the slip surface shape of the safety factor defined by the shear stress ratio form, and is reasonable and effective. The verification result of the calculation example shows that the sliding direction assumption of the method is reasonable, and the safety factor is consistent with the result of the strict three-dimensional limit equilibrium method. Key words three-dimensional slope stability; finite element method; limit equilibrium condition; safety factor; sliding direction 1收稿日期: 基金项目: 国家自然科学基金资助（No. 52079018）工业装备结构分析国家重点实验室自主课题基金（No. S14208） 《工程科学学报》录用稿，https://doi.org/10.13374/j.issn2095-9389.2021.04.29.002 ©北京科技大学 2020 录用稿件，非最终出版稿
边坡稳定分析是岩土工程中的重要课题。通常三维空间中的边坡选取典型断面并假设断面处于平面 应变条件，从而进行二维稳定分析山。二维分析导致结果偏于保守，且不能反映边坡破坏的三维特性。许 多学者开展了三维边坡稳定分析方法的研究，有满足三个方向上的力和力矩平衡的严格三维极限平衡法P ,也有用于三维边坡稳定分析的强度折减法-1。 葛修润四指出极限平衡法和强度折减法都建立在强度折减的概念上，存在许多不合理之处，故在 “真实”应力状态的基础上提出了矢量和法并进行了一系列研究21。“真实”应力状态是指采用数值方 法根据未折减的土体强度参数计算得到的应力场。也有学者在真实应力状态下提出了安全系数为剪应力 比形式的三维稳定分析方法4，该形式的安全系数（式1）理论基础不牢固，仅在滑动面是球形或椭球 形时，安全系数才具有力矩比的物理意义。 t,ds (1) 式中，π为剪应力，为抗剪强度，ds为滑动面微元。 在二维边坡稳定分析中，邵龙潭和李红军剧将一点极限平衡条件扩展到沿滑动面的整体极限平衡条 件，根据该条件定义了安全系数，明确了其物理意义，建立了有限元极限乎衡法，并将其用于各类土工 结构的稳定分析。该方法理论体系严密，应力场真实，能适用于任意形状的滑动面。 本文将有限元极限平衡法扩展到三维，提出了三维滑动面上一点在滑动方向上的极限平衡条件，证 明了沿滑动面的整体极限平衡等价条件。给出了主滑方向的确定方法似及基于主滑方向的滑动方向计算 方法。基于极限平衡条件定义了局部和整体安全系数。最后通过算例验证了本文方法的合理性、有效性 以及对任意形状滑动面的适应性。 1极限平衡条件 1.1一点滑动方向上的极限平衡条件 极限平衡状态是研究对象即将失去而未失去平衡的代态，极限平衡法和强度折减法都通过折减强度 参数的方法，使边坡达到极限平衡状态。极限平衡法中的极限平衡状态为滑动面上每个点在滑动切平面 上的剪应力都等于抗剪强度，强度折减法的极限来衡状态判断则根据失稳判据主要分成3种1：(1)数 值解非收敛(2)塑性区贯通(3)位移突变。采角有限元极限平衡法对边坡进行稳定分析，因为应力场 真实（非极限平衡状态），所以需要对极限平衡状态进行预估。 在二维平面应变问题中，滑动面切平面快定剪应力方向和抗剪强度方向。剪应力方向沿滑动面切向 指向滑动方向。无论是非极限平衡状态还是极限平衡状态，剪应力只有大小变化没有方向变化，如图 1(a) 在三维中，滑动面上一 正应力和剪应力可根据该点处应力张量和滑动面切平面法向量计算： T=G.n a=(T.n).n T=T-Gn (2) 式中，σ为该点处应为张量，为滑动面切平面的法方向单位矢量，T为应力矢量，‘m为法向正应力矢量， π为剪应力矢量。当滑动面确定后，由非极限平衡状态向极限平衡状态的变化过程中，法向正应力矢量方 向不变，但剪应力矢量方向在滑动面的切平面上变化。根据陈惠发的研究，土体发生塑性应变时剪应 变与剪应力可假定为同向。所以假设极限平衡状态时剪应力方向与位移方向一致，即为滑动方向。需注 意的是，各点的滑动方向是在滑动面切平面上的方向，每个点切平面的法方向不同，滑动方向也不同。 若确定了滑动方向，非极限平衡状态下的滑动面应力如图1(b)所示。其中是极限平衡状态的抗剪强度矢 量。用公式表示，在极限平衡状态下，考虑在滑动方向上有： =T.d=(c+r)d (3)
边坡稳定分析是岩土工程中的重要课题。通常三维空间中的边坡选取典型断面并假设断面处于平面 应变条件，从而进行二维稳定分析[1]。二维分析导致结果偏于保守，且不能反映边坡破坏的三维特性。许 多学者开展了三维边坡稳定分析方法的研究，有满足三个方向上的力和力矩平衡的严格三维极限平衡法[2- 6]，也有用于三维边坡稳定分析的强度折减法[7-10]。 葛修润[11]指出极限平衡法和强度折减法都建立在强度折减的概念上，存在许多不合理之处，故在 “真实”应力状态的基础上提出了矢量和法并进行了一系列研究[12,13]。“真实”应力状态是指采用数值方 法根据未折减的土体强度参数计算得到的应力场。也有学者在真实应力状态下提出了安全系数为剪应力 比形式的三维稳定分析方法[14-16]，该形式的安全系数（式 1）理论基础不牢固，仅在滑动面是球形或椭球 形时，安全系数才具有力矩比的物理意义[17]。 d d f s F s （1） 式中，t 为剪应力，tf为抗剪强度，ds 为滑动面微元。 在二维边坡稳定分析中，邵龙潭和李红军[18]将一点极限平衡条件扩展到沿滑动面的整体极限平衡条 件，根据该条件定义了安全系数，明确了其物理意义，建立了有限元极限平衡法，并将其用于各类土工 结构的稳定分析[19]。该方法理论体系严密，应力场真实，能适用于任意形状的滑动面。 本文将有限元极限平衡法扩展到三维，提出了三维滑动面上一点在滑动方向上的极限平衡条件，证 明了沿滑动面的整体极限平衡等价条件。给出了主滑方向的确定方法以及基于主滑方向的滑动方向计算 方法。基于极限平衡条件定义了局部和整体安全系数。最后通过算例验证了本文方法的合理性、有效性 以及对任意形状滑动面的适应性。 1 极限平衡条件 1.1 一点滑动方向上的极限平衡条件 极限平衡状态是研究对象即将失去而未失去平衡的状态，极限平衡法和强度折减法都通过折减强度 参数的方法，使边坡达到极限平衡状态。极限平衡法中的极限平衡状态为滑动面上每个点在滑动切平面 上的剪应力都等于抗剪强度，强度折减法的极限平衡状态判断则根据失稳判据主要分成 3 种[20]：（1）数 值解非收敛（2）塑性区贯通（3）位移突变。采用有限元极限平衡法对边坡进行稳定分析，因为应力场 真实（非极限平衡状态），所以需要对极限平衡状态进行预估。 在二维平面应变问题中，滑动面切平面决定剪应力方向和抗剪强度方向。剪应力方向沿滑动面切向 指向滑动方向。无论是非极限平衡状态还是极限平衡状态，剪应力只有大小变化没有方向变化，如图 1(a)。 在三维中，滑动面上一点的法向正应力和剪应力可根据该点处应力张量和滑动面切平面法向量计算： ( ) n n T = σ n σ T n n τ T σ （2） 式中，σ 为该点处应力张量，n 为滑动面切平面的法方向单位矢量，T 为应力矢量，σn 为法向正应力矢量， τ 为剪应力矢量。当滑动面确定后，由非极限平衡状态向极限平衡状态的变化过程中，法向正应力矢量方 向不变，但剪应力矢量方向在滑动面的切平面上变化。根据陈惠发的研究[21]，土体发生塑性应变时剪应 变与剪应力可假定为同向。所以假设极限平衡状态时剪应力方向与位移方向一致，即为滑动方向。需注 意的是，各点的滑动方向是在滑动面切平面上的方向，每个点切平面的法方向不同，滑动方向也不同。 若确定了滑动方向，非极限平衡状态下的滑动面应力如图 1(b)所示。其中 tf是极限平衡状态的抗剪强度矢 量。用公式表示，在极限平衡状态下，考虑在滑动方向上有： ( ) τ f n T d σ τ d （3） 录用稿件，非最终出版稿
式中，d是一点滑动方向矢量，因为垂直于滑动面，d在滑动面上，所以o。·=0,保持左右两侧一致， 并考虑到与d方向相反，式(3)变为 -ty'd=t…d (4) 式(4)为一点在滑动方向上处于极限平衡状态的条件。 (a) (b) d 图1滑动面应力分析(a)二维(b)三维 Fig.1 Stress on slip surface(a)two-dimensional(b)three-dimensional 1.2滑动面整体极限平衡缘件 任意滑动面上土体整体达到极限平衡状态与滑动面上土体各处在滑动方向上处于极限平衡状态等价。 滑动面整体达到极限平衡状态定义为： das=「-t'dds (5) 式(4)在滑动面上处处成立与式(5)的等价， 证明如下： 如果滑动面上各处都在其自身的滑动方向上处于极限平衡状态，即满足式(3)，则对滑动面上微元 有 t·dds=-tr·dds (6) 对式(6)两侧进行积分，所以式(5)成立 反过来，如果式(5)成立， 转到更分符号内为标量，并且矢量点柔有结合律，则有 (+).dds =0 (7) 因为对于稳定或者处于极限 态的土体，每一点都必有 td≤-trd (8) 双@>0，则要衡7成立， 必须积分域上处处满足： (t+t)d=0 (9) 所以滑动面上每一点都满足式(4)。 需要注意，即使当滑动面上各点的剪应力矢量的大小都等于抗剪强度，整个滑动面也不是处于极限 平衡状态。滑动方向意味着一个运动许可方向，只有当各点的剪应力矢量方向指向滑动方向且大小等于 抗剪强度时，整体滑动面才处于极限平衡状态。 2滑动方向与主滑方向 滑动面上不同点有不同的滑动方向。滑动面上各点滑动方向的计算方法是本文方法的关键点之一。 需要根据现有的有限元计算结果，得到合理的滑动方向。 滑动方向的计算需要用到主滑方向的概念。Kalatehjari等2根据极限平衡法的基本假定之一，刚体滑 动体假定，阐明了三维边坡滑动存在唯一的主滑方向。主滑方向被定义为极限平衡状态下，滑坡体滑动
式中，d 是一点滑动方向矢量，因为σn 垂直于滑动面，d 在滑动面上，所以 n σ d 0，保持左右两侧一致， 并考虑到 tf与 d 方向相反，式（3）变为 f τ d τ d （4） 式（4）为一点在滑动方向上处于极限平衡状态的条件。 图 1 滑动面应力分析(a)二维(b)三维 Fig.1 Stress on slip surface (a) two-dimensional (b) three-dimensional 1.2 滑动面整体极限平衡条件 任意滑动面上土体整体达到极限平衡状态与滑动面上土体各处在滑动方向上处于极限平衡状态等价。 滑动面整体达到极限平衡状态定义为： d d f s s τ d τ d （5） 式（4）在滑动面上处处成立与式（5）的等价，证明如下： 如果滑动面上各处都在其自身的滑动方向上处于极限平衡状态，即满足式（3），则对滑动面上微元 有 d d f τ d s s τ d （6） 对式（6）两侧进行积分，所以式（5）成立。 反过来，如果式（5）成立，考虑到积分符号内为标量，并且矢量点乘有结合律，则有 ( ) d f s τ τ d 0 （7） 因为对于稳定或者处于极限平衡状态的土体，每一点都必有 f τ d τ d （8） 且ds 0，则要使式（7）成立，必须积分域上处处满足： ( ) f τ τ d 0 （9） 所以滑动面上每一点都满足式（4）。 需要注意，即使当滑动面上各点的剪应力矢量的大小都等于抗剪强度，整个滑动面也不是处于极限 平衡状态。滑动方向意味着一个运动许可方向，只有当各点的剪应力矢量方向指向滑动方向且大小等于 抗剪强度时，整体滑动面才处于极限平衡状态。 2 滑动方向与主滑方向 滑动面上不同点有不同的滑动方向 d。滑动面上各点滑动方向的计算方法是本文方法的关键点之一。 需要根据现有的有限元计算结果，得到合理的滑动方向。 滑动方向的计算需要用到主滑方向的概念。Kalatehjari 等[22]根据极限平衡法的基本假定之一，刚体滑 动体假定，阐明了三维边坡滑动存在唯一的主滑方向。主滑方向被定义为极限平衡状态下，滑坡体滑动 (a) (b) 录用稿件，非最终出版稿
方向的水平投影方向。主滑方向不同于滑动方向，主滑方向是定义在XOY平面上，每一个滑动面只对应 一个主滑方向。现对唯一主滑方向的合理性进行说明： 对滑动体进行垂直条分，根据假定每个条分柱体都是刚体。在极限平衡状态下，假设每个柱体运动 方向的水平投影有三种情况：平行、聚拢和分离。刚体假设阻止了柱体聚拢，而柱体分离则意味着柱体 间没有力的相互作用，力和力矩平衡被打破，这与实际情况不符。所以互相平行的柱体运动方向水平投 影是唯一的选项。 在三维极限平衡法中，很多方法对条柱底面的极限剪应力方向也有类似假定，比如要求剪应力平行 于XOZ平面23,2或剪应力在XOY平面内的投影方向保持一致562四。极限平衡法是通过折减材料强度参数 获得满足力和力矩平衡的极限状态从而计算安全系数的方法，底面剪应力方向就是极限状态时的滑动方 向。在非极限状态下，本文根据真实应力场计算滑动面上的剪应力没有这种方向一致性，但通过主滑方 向计算得到的滑动方向，正是对极限状态时滑动方向的预估，也是对极限状态时剪应力方向的预估。 滑动面上各点的滑动方向可以根据唯一的主滑方向计算。根据滑动方向处于滑动面上；滑动方向在 XOY平面的投影是主滑方向：滑动方向是单位矢量，可以联立方程计算滑动方向单矢量： d.n=0 (10) 式中山是主滑方向矢量，k是Z方向的单位矢量。三个方程可以求解滑动方向矢量d的三个元素。 通过非极限平衡状态的滑动面应力场对极限平衡状态推测得到主滑方向，可以将滑动面剪应力矢量 的和矢量的水平投影方向作为主滑方向，用公式表示为： xds(110) xd(110) (11) 主滑方向的选取具有一定的主观性，比如设补人员认为边坡在未来会受到某一方向的荷载，从而产 生沿某一特定方向的破坏，那么要考虑这一特定方向的安全系数，就可以将这一特定方向设为主滑方向 进行计算。 3局部和整体安全系数定义 安全系数是对研究对象距离极限乎衡状态的描述，安全系数等于1则认为研究对象处于极限平衡状 态。根据对一点在滑动方向上处于极限平衡状态的定义（式3），一点处的局部安全系数被定义为： F=T1d rd (12) 式(12)表示 点安全系数是抗剪强度与剪应力在滑动方向上的投影的比，其物理意义是一点在滑 动方向上距离极限乎衡状态的衡量，也是一点在滑动方向上强度的储备。 F(S)是在滑动面上的局部安全系数函数，在各点符合式(12)的定义，可使滑动面上土体微元均 达到极限平衡状态的，那么土体在滑动面S上整体达到极限平衡就变为 r-器 (13) 应用积分中值定理，式(13)右侧变为： (14) 则有
方向的水平投影方向。主滑方向不同于滑动方向，主滑方向是定义在 XOY 平面上，每一个滑动面只对应 一个主滑方向。现对唯一主滑方向的合理性进行说明： 对滑动体进行垂直条分，根据假定每个条分柱体都是刚体。在极限平衡状态下，假设每个柱体运动 方向的水平投影有三种情况：平行、聚拢和分离。刚体假设阻止了柱体聚拢，而柱体分离则意味着柱体 间没有力的相互作用，力和力矩平衡被打破，这与实际情况不符。所以互相平行的柱体运动方向水平投 影是唯一的选项。 在三维极限平衡法中，很多方法对条柱底面的极限剪应力方向也有类似假定，比如要求剪应力平行 于 XOZ 平面[23,24]或剪应力在 XOY 平面内的投影方向保持一致[5,6,22]。极限平衡法是通过折减材料强度参数 获得满足力和力矩平衡的极限状态从而计算安全系数的方法，底面剪应力方向就是极限状态时的滑动方 向。在非极限状态下，本文根据真实应力场计算滑动面上的剪应力没有这种方向一致性，但通过主滑方 向计算得到的滑动方向，正是对极限状态时滑动方向的预估，也是对极限状态时剪应力方向的预估。 滑动面上各点的滑动方向可以根据唯一的主滑方向计算。根据滑动方向处于滑动面上；滑动方向在 XOY 平面的投影是主滑方向；滑动方向是单位矢量，可以联立方程计算滑动方向单位矢量： u d n d d k k d d d k k d 0 1 （10） 式中 du是主滑方向矢量，k 是 Z 方向的单位矢量。三个方程可以求解滑动方向矢量 d 的三个元素。 通过非极限平衡状态的滑动面应力场对极限平衡状态推测得到主滑方向，可以将滑动面剪应力矢量 的和矢量的水平投影方向作为主滑方向，用公式表示为： d ( , , ) d ( , , ) s u s s s τ d τ 11 0 11 0 （11） 主滑方向的选取具有一定的主观性，比如设计人员认为边坡在未来会受到某一方向的荷载，从而产 生沿某一特定方向的破坏，那么要考虑这一特定方向的安全系数，就可以将这一特定方向设为主滑方向 进行计算。 3 局部和整体安全系数定义 安全系数是对研究对象距离极限平衡状态的描述，安全系数等于 1 则认为研究对象处于极限平衡状 态。根据对一点在滑动方向上处于极限平衡状态的定义（式 3），一点处的局部安全系数被定义为： f F τ d τ d （12） 式（12）表示一点安全系数是抗剪强度与剪应力在滑动方向上的投影的比，其物理意义是一点在滑 动方向上距离极限平衡状态的衡量，也是一点在滑动方向上强度的储备。 F s( )是在滑动面 s 上的局部安全系数函数，在各点符合式（12）的定义，可使滑动面上土体微元均 达到极限平衡状态的，那么土体在滑动面 s 上整体达到极限平衡就变为 d d ( ) f s s F s τ d τ d （13） 应用积分中值定理，式（13）右侧变为： d d ( ) f f G s s F F s τ d τ d 1 （14） 则有 录用稿件，非最终出版稿
dds t·dds (15) F。是整体安全系数，表征滑动面上各点局部安全系数的中值。在推导整体安全系数的定义中，对滑 动面形状和对称性未做要求，所以方法可以用于任意形状滑动面的安全系数计算。 二维情况下，因为t和d方向相同，T和d方向相反，式(10)退化为 rdl (16) 这是在基于应力的二维边坡稳定分析方法中被广泛使用的全局安全系数定义式，也是有限元极限平衡法 的安全系数定义。 三维情况下，对比安全系数定义式(1)和式(15)，如果假定π和方简相同和方向相反， 式(15)可以退化为式(1)。这种假定隐含剪应力矢量方向在非极限平衡状态和极限平衡状态相同的设 定，对倾向滑动方向的滑动面，剪应力矢量方向变化不大，但对倾向方向编离滑动方向的滑动面，剪应 力矢量方向在非极限平衡状态和极限平衡状态差距很大，这会导致误差， 这也是采用式(1)定义安全系 数的限制。 4计算方法 本文方法需要计算边坡区域的应力场，并基于应力计算边坡的安全系数。有限元法可以根据本构关 系计算应力场。在本文中，有限元本构使用线性弹性理想塑件模型，计算软件为商业有限元软件Abaqus 6.l4。理想的塑性屈服准则与Mohr-Coulomb破坏准则致。抗剪强度采用Mohr-Coulomb计算： T=o tano+c (17) 式中，口m为法向正应力，P为摩擦角，c为粘聚 在计算过程中需要对滑动面上的变量进行积分。数值积分需要对滑动面进行离散，对滑动面进行网 格划分。在这里选择将滑动面划分为乏角形网格，原因是三角形网格比四边形网格更容易地计算滑动面 的法线方向并更准确地描述滑动面，在边界处也更容易分割。通过三个步骤构造滑动面，如图2所示。 第一步：在XOY平面中创建三角形网格： 第二步：通过插值将三角形网格节点的z坐标从XOY平面调整到构造表面的位置： 第三步：获得坡面与三角形网格的交点，三角形网格被分割，模型内部的部分被保留，作为滑动表 面。 另一种方法是，茹果滑动面在有限元模型中已经被划分（比如滑动面是岩体结构面或确定的软弱 层)，则可以真接将滑动面位置的四面体单元表面的三角形网格提取出来，作为滑动面网格使用：如果 有限元采用的是六前体单元，则可以将四边形表面从中剖分，得到两个三角形网格单元。 (a)
d d f G s F s τ d τ d （15） FG 是整体安全系数，表征滑动面上各点局部安全系数的中值。在推导整体安全系数的定义中，对滑 动面形状和对称性未做要求，所以方法可以用于任意形状滑动面的安全系数计算。 二维情况下，因为τ 和d 方向相同， f τ 和d 方向相反，式（10）退化为 d d f l F l （16） 这是在基于应力的二维边坡稳定分析方法中被广泛使用的全局安全系数定义式，也是有限元极限平衡法 的安全系数定义[18]。 三维情况下，对比安全系数定义式（1）和式（15），如果假定 τ 和 d 方向相同，tf和 d 方向相反， 式（15）可以退化为式（1）。这种假定隐含剪应力矢量方向在非极限平衡状态和极限平衡状态相同的设 定，对倾向滑动方向的滑动面，剪应力矢量方向变化不大，但对倾向方向偏离滑动方向的滑动面，剪应 力矢量方向在非极限平衡状态和极限平衡状态差距很大，这会导致误差，这也是采用式（1）定义安全系 数的限制。 4 计算方法 本文方法需要计算边坡区域的应力场，并基于应力计算边坡的安全系数。有限元法可以根据本构关 系计算应力场。在本文中，有限元本构使用线性弹性理想塑性模型，计算软件为商业有限元软件 Abaqus 6.14。理想的塑性屈服准则与 Mohr-Coulomb 破坏准则一致。抗剪强度采用 Mohr-Coulomb 计算： n tan f c （17） 式中， n 为法向正应力， 为摩擦角，c 为粘聚力。 在计算过程中需要对滑动面上的变量进行积分。数值积分需要对滑动面进行离散，对滑动面进行网 格划分。在这里选择将滑动面划分为三角形网格，原因是三角形网格比四边形网格更容易地计算滑动面 的法线方向并更准确地描述滑动面，在边界处也更容易分割。通过三个步骤构造滑动面，如图 2 所示。 第一步：在 XOY 平面中创建三角形网格； 第二步：通过插值将三角形网格节点的 z 坐标从 XOY 平面调整到构造表面的位置； 第三步：获得坡面与三角形网格的交点，三角形网格被分割，模型内部的部分被保留，作为滑动表 面。 另一种方法是，如果滑动面在有限元模型中已经被划分（比如滑动面是岩体结构面或确定的软弱 层），则可以直接将滑动面位置的四面体单元表面的三角形网格提取出来，作为滑动面网格使用；如果 有限元采用的是六面体单元，则可以将四边形表面从中剖分，得到两个三角形网格单元。 录用稿件，非最终出版稿 (a)
(b) (c) 图2滑动面构造步骤.(a)三角形网格：(b)网格节点调整：(c)滑动面与坡面相交 Fig.2 Construction steps of the slip surface:(a)Triangle mesh;(b)Mesh node adjustment()Slip surface intersects the slope surface 滑动面三角形网格中心点应力张量由有限元网格节点处应力采量插值计算，插值方法选用三维线性 插值。本文方法安全系数计算流程如图3所示。 Calculate stress by finite element Calculate the stress tensor on the sliding surface by interpolation Calculate the normal stress and shear stress vector by Eq.2 Cuate the main sliding direction by Eq.11 录用 Calculate the sliding direction of each points by Eq.10 Calculate local and global safety factors by Eq.12 and Eq.15 图3安全系数计算流程 Fig.3 Flowchart of safety factor calculation 5算例 算例中将采用式(1)剪应力比形式定义的安全系数记作F1,将采用式(15)定义的安全系数记作 F2。 5.1算例1 算例I来自Fredlund和Krahn2对二维边坡的研究，有学者将这一算例拓展到三维情况并进行了研究。 该算例包含两种工况，一种是均质边坡，滑动面为椭球形；一种是带软弱薄层的边坡，滑动面为椭球型 和位于软弱薄层的平面组合。椭圆的长轴为133.8m,有限元模型y轴方向的长度为110m,其他相关参
图 2 滑动面构造步骤. (a)三角形网格；(b)网格节点调整；(c)滑动面与坡面相交 Fig.2 Construction steps of the slip surface: (a) Triangle mesh; (b) Mesh node adjustment; (c) Slip surface intersects the slope surface 滑动面三角形网格中心点应力张量由有限元网格节点处应力张量插值计算，插值方法选用三维线性 插值。本文方法安全系数计算流程如图 3 所示。 图 3 安全系数计算流程 Fig.3 Flowchart of safety factor calculation 5 算例 算例中将采用式（1）剪应力比形式定义的安全系数记作 F1，将采用式（15）定义的安全系数记作 F2。 5.1 算例 1 算例 1 来自 Fredlund 和 Krahn[25]对二维边坡的研究，有学者将这一算例拓展到三维情况并进行了研究。 该算例包含两种工况，一种是均质边坡，滑动面为椭球形；一种是带软弱薄层的边坡，滑动面为椭球型 和位于软弱薄层的平面组合。椭圆的长轴为 133.8 m，有限元模型 y 轴方向的长度为 110 m，其他相关参 (b) (c) 录用稿件，非最终出版稿
数如图4所示。横向边界条件为约束y轴方向位移（平面应变）。图5展示滑动面网格和有限元网格尺寸 对安全系数的影响。随着滑动面网格尺寸从4m减小到0.75，安全系数逐渐减小并稳定，有限元网格 尺寸从4m减小到1m,安全系数略有减小。图6和图7分别展示工况1的滑动方向和剪应力矢量方向， 通过对比可见两者方向基本一致，只在侧面靠近坡面处略有差别。因此表1中F1与F2在工况1下基本一 致。 (a) 20 15 a83 E-10P%-049 月-l22n8-26 20 10 (b) 非最出版稿 (c) 图4算例1.(a材料参数和断面：b)工况1：(c)工况2 Fig.4 Example 1:(a)Material parameters and cross-section;(b)Case 1;(c)Case 2 (a) (b) 2.30 230 1190 -Case】 --Case2 3 1B5 225 1.85 0.75 43 2 2.20 泰用稿件 1.80 220 1.80 1.75 215 1.76 075 3 2 110 1.70 210 1.70 1000 10000 10000 100000 Nunber of Slip Surface Elements Number of FEM Elements 图5单元数对安全系数的影响.(a)滑动面：(b)有限元 Fig.5 The influence of the number of units on the safety factor:(a)Slip Surface;(b)FEM
数如图 4 所示。横向边界条件为约束 y 轴方向位移（平面应变）。图 5 展示滑动面网格和有限元网格尺寸 对安全系数的影响。随着滑动面网格尺寸从 4 m 减小到 0.75 m，安全系数逐渐减小并稳定，有限元网格 尺寸从 4 m 减小到 1 m，安全系数略有减小。图 6 和图 7 分别展示工况 1 的滑动方向和剪应力矢量方向， 通过对比可见两者方向基本一致，只在侧面靠近坡面处略有差别。因此表 1 中 F1与 F2在工况 1 下基本一 致。 图 4 算例 1. (a)材料参数和断面；(b)工况 1；(c)工况 2 Fig.4 Example 1: (a) Material parameters and cross-section; (b) Case 1; (c) Case 2 图 5 单元数对安全系数的影响. (a)滑动面；(b)有限元 Fig.5 The influence of the number of units on the safety factor: (a) Slip Surface; (b) FEM (a) (b) (c) (a) (b) 录用稿件，非最终出版稿
20 (a) (b) 120 100 80 40 稿 图6工况1的滑动方向与主滑方向.(a)等轴测图心)俯视图 Fig.Sliding direction and main sliding direction of case 1(a)Isometric view(b)Top view 20、 10、 (a) 0 20 1的应为天靴广人 100 60 40 20 0 图7 工况1的剪应力矢量方向与主滑方向.(a)等轴测图：(b)俯视图 表1算例1的安全系数比较 Table 1 A comparison between the safety factors computed for example 1 Method Case 1 Case2 Rigorous limit equilibrium 1.707 Rigorous limit equilibrium 2.127 1.766 F 2.213 1.6074 F 2.226 1.730 5.2算例2 楔形体破坏是三维岩石边坡破坏的常见情况。如图8所示，这里使用两个典型的楔形破坏算例来验 证本文方法的合理性。算例包括具有对称结构表面的楔形体和具有非对称结构表面的楔形体。算例的材
图 6 工况 1 的滑动方向与主滑方向. (a)等轴测图(b)俯视图 Fig.6 Sliding direction and main sliding direction of case 1 (a) Isometric view (b) Top view 图 7 工况 1 的剪应力矢量方向与主滑方向. (a)等轴测图；(b)俯视图 Fig.7 Shear stress vector direction and main sliding direction of case 1: (a) Isometric view; (b) Top view 表 1 算例 1 的安全系数比较 Table 1 A comparison between the safety factors computed for example 1 Method Case 1 Case 2 Rigorous limit equilibrium[3] - 1.707 Rigorous limit equilibrium[6] 2.127 1.766 F1 2.213 1.607[14] F2 2.226 1.730 5.2 算例 2 楔形体破坏是三维岩石边坡破坏的常见情况。如图 8 所示，这里使用两个典型的楔形破坏算例[26]来验 证本文方法的合理性。算例包括具有对称结构表面的楔形体和具有非对称结构表面的楔形体。算例的材 (a) (b) (a) (b) 录用稿件，非最终出版稿
料参数和几何信息列于表2和表3中。通过图9和图10的对比可以明显看出，在非极限平衡状态下，剪 应力矢量方向与滑动方向不一致，所以导致表4中F1结果小于其他方法安全系数。Hoek-Bry解是一个静 力平衡解，一般被认为是该算例的理论解，本文方法和严格极限平衡法都与理论解一致，但F小于理论 解，有较大误差。 表2算例2的材料参数 Table 2 Mechanical parameters used in example 2 Symmetric wedge Asymmetric wedge Rock Structural surfaces Rock Structural surfaces Weight X(kN.m) 26 26 20 Poisson ratio u 0.49 0.49 0.49 0.49 Friction angle ( 45 20 45 30 Cohesion c(kPa) 1000 20 1000 50 表3算例2的几何信息 Table 3 Geometric information used in example 2 Symmetric wedge Asymmetrie wedge Dip direction() Dip() Dip() Left structural surface 115 45 120 40 Right structural surface 245 240 60 Upper surface 180 o 180 0 Slope surface ，非最终 180 60 180 60 (a) 录用稿 (b) 图8算例2.楔形体破坏()对称楔形体(b)非对称楔形体 Fig.8 Example 2:wedge failure(a)Symmetric wedge(b)Asymmetric wedge
料参数和几何信息列于表 2 和表 3 中。通过图 9 和图 10 的对比可以明显看出，在非极限平衡状态下，剪 应力矢量方向与滑动方向不一致，所以导致表 4 中 F1结果小于其他方法安全系数。Hoek-Bray 解是一个静 力平衡解，一般被认为是该算例的理论解，本文方法和严格极限平衡法都与理论解一致，但 F1小于理论 解，有较大误差。 表 2 算例 2 的材料参数 Table 2 Mechanical parameters used in example 2 Symmetric wedge Asymmetric wedge Rock Structural surfaces Rock Structural surfaces Weight (kN•m-3) 26 20 26 20 Poisson ratio 0.49 0.49 0.49 0.49 Friction angle (°) 45 20 45 30 Cohesion c(kPa) 1000 20 1000 50 表 3 算例 2 的几何信息 Table 3 Geometric information used in example 2 Symmetric wedge Asymmetric wedge Dip direction(°) Dip(°) Dip direction(°) Dip(°) Left structural surface 115 45 120 40 Right structural surface 245 45 240 60 Upper surface 180 10 180 0 Slope surface 180 60 180 60 图 8 算例 2. 楔形体破坏(a)对称楔形体(b)非对称楔形体 Fig.8 Example 2: wedge failure (a) Symmetric wedge (b) Asymmetric wedge (a) (b) 录用稿件，非最终出版稿
≥Sliding direction Main sliding directior 100 50 200 100 150 (a) 50 100 0 出版稿 Sliding direction Main sliding direction 200 (b) 150 100 50 0 图9非对称楔形体的滑动方向与主滑方向.(a)等轴测图b)俯视图 Fig.9 Sliding direction and main sliding direction of an asymmetrical wedge (a)Isometric view(b)Top view vector direction 录用 150 200 100 150 (a) 100
图 9 非对称楔形体的滑动方向与主滑方向. (a)等轴测图(b)俯视图 Fig.9 Sliding direction and main sliding direction of an asymmetrical wedge (a) Isometric view (b) Top view (a) (b) (a) 录用稿件，非最终出版稿