岩土边坡稳定性的刚体有限元上限分析法

更新时间：2023-04-24 14:25:01 阅读量：实用文档文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

岩土边坡稳定性分析推荐度：
相关推荐

刚体有限元法与极限分析相结合为求解边坡稳定性问题提供了新的方法。采用刚体有限单元离散边坡计算区域,同时构造运动许可速度场,在满足屈服条件、流动法则、虚功方程以及相应的边界条件的基础上,引入非线性数学规划方法求解最小安全系数。几个算例说明了该方法的正确性和可

第23卷第6期

岩石力学与工程学报 23(6)：898～905

2004年3月 Chinese Journal of Rock Mechanics and Engineering March，2004

岩土边坡稳定性的刚体有限元上限分析法

殷建华陈健李焯芬

(香港理工大学土木及结构工程系香港) (香港大学土木工程系香港)

摘要刚体有限元法与极限分析相结合为求解边坡稳定性问题提供了新的方法。采用刚体有限单元离散边坡计算区域，同时构造运动许可速度场，在满足屈服条件、流动法则、虚功方程以及相应的边界条件的基础上，引入非线性数学规划方法求解最小安全系数。几个算例说明了该方法的正确性和可行性。关键词工程地质，边坡稳定，极限分析，上限法，刚体有限元，非线性规划

分类号 TU 413.6+2 文献标识码 A 文章编号 1000-6915(2004)06-0898-08

UPPER LIMIT ANALYSIS OF STABILITY OF ROCK AND SOIL SLOPES

USING RIGID FINITE ELEMENTS

Yin Jianhua，Chen Jian1，Lee C F2

(1Department of Civil and Structural Engineering，The Hong Kong Polytechnic University， Hong Kong， China)

(2Department of Civil Engineering， The University of Hong Kong， Hong Kong， China)

Abstract In this paper，the development and application of a new upper limit method for 2D and 3D slope stability problems are presented. Rigid finite elements are used to construct a kinematically admissible velocity field. The proposed method formulates the slope stability problems as an optimisation problem based on the upper bound theorem. The objective function for determination of the minimum value of the safety factor has some unknowns, which are subject to a set of linear and non-linear equality constraints derived from an energy-work balance equation，the Mohr-Coulomb failure (yield) criterion，an associated flow rule and a number of boundary conditions. The objective function with constraints leads to a non-linear programming problem， which can be solved by a sequential quadratic algorithm. Four typical 2D and 3D slope stability problems selected from the literature are analysed with the presented method. The results of the presented limit analysis are compared with those produced by other approaches reported in the literature.

Key words engineering geology，slope stability，limit analysis，upper bound，rigid finite element，non-linear programming

1 前言

边坡稳定性分析是岩土工程中的重要研究课题。极限平衡法是工程实践中最为广泛采用的边坡稳定性分析方法，其原因在于它原理简单，计算方便，有比较成熟的应用经验。但是，这类方法是建

2002年7月23日收到初稿，2002年9月26日收到修改稿。 * 香港特区大学资助委员会资助项目(PolyU 5064/00E)。

立在静力平衡基础上的，对于滑动机构，需引入内力假设使之成为静定结构，因此，极限平衡法始终只是一种近似方法。

文[1]提出的基于塑性极限理论的极限分析方法能够为稳定性分析问题提供贴近理论的现实解答。文[2]全面阐述了用上限定理来求解地基承载力、土压力和边坡稳定性的原理和方法。这类方法

作者殷建华简介：男，1956生，博士，现任教授，主要从事软土工程、边坡分析与监测以及室人试验等方面的研究工作。E-mail：cejhyin@polyu.edu.hk。

第23卷第6期殷建华等. 岩土边坡稳定的刚体有限元上限分析法 899

具有坚实的理论基础，利用经典塑性力学的上、下限定理从下限和上限2个方向逼近真实解。在上限和下限分析中，其各自的关键所在是运动许可速度场和静力许可应力场的构造技术及其优化分析。由于较难找到合适的静力许可的应力分布从而获取下限解，因此，极限分析方法中最常用的是上限方法。文[3，4]应用塑性力学上限定理，从变形协调出发，建立运动许可速度场，根据外力功和内能耗散相平衡的原理确定安全系数，这种方法也称之为能量法。文[5，6]在能量法的基础上，探讨了采用非关联流动法则对二维土质边坡稳定安全系数的影响。最近，文[7，8]等改进和扩展了能量法，使其能够进行三维边坡的稳定性分析。

文[9～13]在结合有限单元法和数学规划方法求解二维和三维稳定性问题(主要是承载力问题)方面进行了卓有成效的工作。其基本思路是，利用有限单元法将岩土结构物离散，建立运动许可速度场和静力许可应力场，在满足屈服条件、流动法则、边界条件、平衡条件或虚功方程的前提下，借助线性或非线性数学规划方法得到极限荷载的严格上限解和下限解。

文[14]利用刚体有限元和极限分析下限定理对边坡进行了稳定性评价。刚体有限元的建模思想源于文[15]提出的刚体弹簧元模型，文[16，17]后来对其进行了发展和改进，相应地也存在着许多不同的名称，如称之为刚体弹簧元、刚体有限元、刚性有限元、刚体元、刚性元、界面元等等，本文统称之为刚体有限元。

本文介绍了一种应用极限分析的上限定理，借助于刚体有限元和非线性数学规划来进行二维和三维边坡稳定性的评价分析方法

[18]

。首先，用刚体有

限单元来离散边坡体，单元的刚体位移作为基本未知变量，速度间断场允许发生在相邻单元之间的界面上。为保证速度场运动许可，未知变量必须满足屈服条件、流动法则、虚功方程以及相应的边界条件，这些约束条件就构成了非线性规划问题，求其解即得边坡的最小安全系数。在本文中对这个非线性规划问题采用了一种序列二次规划法[19]

进行求解。

2 计算原理与方法

2.1 刚体有限元

刚体有限元法可以说是介于极限平衡条分法和有限单元法之间，兼有二者的优点。刚体有限元的

离散过程和支配方程均类似于有限元，只是在刚体有限元的离散模型中，单元被视为刚体，并以单元的形心位移(速度)作为基本未知量。相邻单元之间由界面连接，刚体单元本身不发生变形，其变形是通过单元之间的相对运动来体现，结构的变形能完全贮存在单元之间的界面上。

应用上限定理进行边坡稳定性分析，首先将计算区域按刚体单元进行离散，本文在二维情况下采用三角形单元，在三维情况下采用四面体单元。速度间断场允许发生在相邻单元之间的界面上。假定D为问题的维数(D = 2表示二维，D = 3表示三维)，则每个刚体单元形心处的速度vg为一个D(D+1)/2维矢量，即

vg={vi}T(i=1，L，D(D+1)/2)

单元内任一点p(x，y，z)的速度v(x，y，z)可表示为

v(x，y，z)=Nvg (1)

式中：

N为形函数矩阵，在二维条件下：

N= 10 (y yg) 01(x x (2) g)

在三维条件下：

N= 100

0(z zg) (y yg) 010 (z zg)

0(x x

g) (3) 001(y yg) (x xg)0

为了方便分析，引入了一个局部坐标系n-d-s，如图1所示，局部坐标轴n沿交界面外法线方向，局部坐标轴d沿交界面的倾向(即交界面的最陡下降方向)，而局部坐标轴s则沿交界面的走向方向。局部坐标系遵从右手螺旋法则。

图1 局部坐标系示意：n(法向)，d(倾向)，s(走向) Fig.1 Local coordinate system defined by n (normal

direction)，d (dip direction) and s (strike direction)

900 岩石力学与工程学报 2004年

以三维问题为例，图2表示了2个相邻的四面体单元，在整体坐标系下，其速度分别为v(1)和v(2)。以单元1的局部坐标系为例，交界面上任一点P的

速度用Δv(2 1)local

表示，则有Δv(2 1)(2 1)

(2 1)local=[Δvn，Δvd， Δv

(2 1)]T

(2 1)s

。其中，Δv

，Δv

(2 1)d

和Δv(2 1)s

分别表示

沿局部坐标轴n，d和s的速度分量。点P在单元1和2中的速度分别可表示为v(1)和v

(2)locallocal

，则其相对

速度为Δv

(2 1)local

(2)1)local

(local

。通过局部坐标和整体坐标的转换关系，可导出

Δv

(2 1)(2))1)(1)

local

(1local

=L(1)v

(2)

L(v (4)

式中：L(1)为由单元１的局部坐标n-d-s和整体坐标夹角方向余弦组成的转换矩阵，即

Fig.2 Three-dimensional velocity discontinuity

cos(n，x)cos(n，y)cos(n，z)

L(1)= cos(d，x)cos(d，y)cos(d，z) xs (5)

cos(s，)cos(s，y)cos(，z) 将式(1)代入式(5)，则有

Δv

(2 1)

2)(2)

)

local

=L(1)

[N(v

N(1)v

(1g

] (6)

写成一般形式为

Δv(2 1)local=A1vG (7) 其中，

(1)

L(1)

]

N(2)

1=[L 0

N(1)

v v(1)g G= v(2)

2.2 相关约束条件和边界条件

在岩土工程计算分析中常用Mohr-Coulomb屈服函数，其表达式为

|τ|=c′+σ′n

tan ′ (8) 与有限单元法计算高斯点应力不同，刚体有限元法计算交界面上的面力，可以直接将面力分量代入式(8)进行屈服破坏判断。这种处理方法实际上假定了单元交界面即为破坏面。这在一定程度上决定了刚体有限元法的计算结果与网格的剖分方式相关。

速度间断场允许发生在相邻三角形或四面体单元的共享边或面上。为使速度场运动许可，速度必须满足塑性流动法则。按照Mohr-Coulomb屈服函数和相关联流动法则，法向速度分量与切向速度分量之间必须满足如下关系：

Δvn= Δvttan ′ (9)

显然，由于绝对值的原因，Δvt可为零、负值或正值。从数学规划的角度来看，这是一个自由变量，任何一个自由变量都可以表示成2个非负数之差。因此，速度间断场上的相对切向速度分量Δvt可以用2个非负变量v+和v 来表示：

Δvt=v+ v (10)

其中，

Δvt={ΔvdΔvs}T

v+={v++T

1LvD 1} v ={v ，v T1，LD 1}

其中，

v+i ≥0

≥0 (i = 1，…，D－1) (11) i

本文采用了与文[11，12]同样的处理方法，|Δvt|可以表示为

D 1

Δvt=∑(v+i+v i) (12)

i=1这样，式(10)～(12)可以写成如下矩阵形式：

Δv=A2vd

v (13)

d≥0

其中，

v+ + T

d={v1v1LvD 1vD 1}

在二维情况下：

A tan ′ tan ′ 2=

1 1

第23卷第6期殷建华等. 岩土边坡稳定的刚体有限元上限分析法 901

在三维情况下：

tan ′ tan ′ tan ′ tan ′

A2= 1 100

001 1

运动许可的速度场在速度边界上必须满足相应的边界条件。假设单元k位于边界上，其速度已知为，则对机动容许的速度场必须满足以下等式：

= (14)

2.3 等效荷载

由于在刚体有限元模型中，以单元形心处的刚体位移作为基本未知量，因此，外力荷载应该转换为单元形心的等效荷载。在这里，可以借助自然坐标来简化其计算。

在一维问题中任意一点P的自然坐标Li(i= 1，2)定义为

L1=l2/l

L=l/l (15)

式中：l1，l2，l如图3所示，其与直角坐标的变换关系为

x=L1x1+L2x2

y=L (16)

1y1+L2y2

一维自然坐标有如下的性质：

∫∫Lab

a!b!1L2dl=L(a+b+1)!L (17)

图3 自然坐标 Fig.3 Natural coordinates

外力荷载一般有面力荷载和体力荷载。面力等效荷载的计算可以借助一维自然坐标。如图4所示，有一均布面力q=[0 q]T

，其等效形心荷载的计

算如下：

Q=∫

lNTqdl (18)

将式(16)，(17)代入式(18)，可得

Q=[0 qlT

qlAB(xg xc)]

(19)

式中：lAB为AB边的长度，xc为AB边中点的x方向坐标。

图4 三角形边上受均布面力

Fig.4 Uniformly distributed load on triangle edge

式(19)表明，作用在AB边上的均布面力的等效形心荷载沿y方向的分量为该边面力的合力，沿转角方向的分量则相当于该边面力的合力对形心点的力矩。这个成果也符合刚体作用力的平移定理。

体力的等效形心荷载可以类似于上述方法，借助面积坐标来计算。 2.4 虚功方程

根据虚功原理，对于任意一组运动许可的速度场，外力的虚功应当等于物体内能的耗散，即

∫

σ′ε&** *

ijd +∫Γ

*σΓ′ε***

ij′&Γ*dΓ=WV

+QV*+PV* (20) 式中：左边2项分别为产生于破坏体 *内和沿滑动面Γ*的内部耗散能；右边3项则分别代表了滑动体自重、所承受的面力以及孔隙水压力的等效形心荷载在虚速度V*上所做的功。

按照刚体有限元的假设，刚体单元本身没有应变，单元内部不存在内能耗散，因而式(20)左边第

1项等于零，即∫

*σ′ijε&*ijd *=0。系统的内能仅仅产

生于刚体单元之间的交界面上，可以表示为

∫

σΓ′**

Γ*

′ε&Γ*dΓ

=∫ SΔv

+σd

′n

Δvn)dS (21) 根据Mohr-Coulomb屈服准则和相关联流动法则，式(21)右边项可以写成

∫ S(Δv

+σn

ΔvS=d

′n)d∫ Sc′ΔvtdS= d

∫ D 1 Sc′(v++v

) d ∑ii dS (22) i=1

由式(22)，可以得到

∫ D 1+ ** Sc′(vd ∑i+v

i) dS=Wv*g+Qvg+Pvg (23) i=1

写成一般的矩阵形式为

902 岩石力学与工程学报 2004年

A3vd=A4vg (24)

式中：AT3={ci′li} (i=1，L，nD)，li为速度间断线(面)的长度或面积，nD为速度间断线(面)的总数；A4=W+Q+P。

2.5 目标函数

工程实践中，边坡的稳定性一般采用安全系数作为定量评价指标。边坡沿着某一滑裂面滑动的安全系数可以这样定义，将岩土的抗剪强度指标c′和

′降低为

c′e

=c′/F

tan ′ e=tan ′/F (25)

则岩土体沿着此滑裂面处处达到极限平衡。这种将强度指标的储备作为安全系数定义的基础方法是经过多年的实践，并被工程界广泛承认的一种作法。

塑性力学上限定理指出，通过按外力功和内能耗散相等来计算获得的外荷载一定比真实的临界荷载大或与其相等。在求解边坡稳定问题时，由虚功方程计算确定的安全系数将比真实解大或与其相等。因此，从数学角度看，边坡稳定问题的上限分析可以视为一个确定安全系数的极小值问题。从实用角度来看，则需应用从最优化原理发展的数值计算方法来确定最小安全系数。 2.6 最优化问题

以上介绍了将边坡稳定问题的上限分析转化为最优化问题求解的基本步骤，值得注意的是抗剪强

度折减参数c′e

和 e′出现在约束条件中，也就是说，作为未知量的安全系数以一种非线性形式出现在约束方程中。

对所有单元，组集以上方程，由此可以得到求解边坡稳定问题的最小安全系数的数学模型。求最小值F，其约束条件为

A1vG=A2vd

3vd=A4vg

(26) g=

v≥0

其中，

tan ′e

tan ′e tan ′e A3= 1

100 00

A4={ci′eli}

(i=1，L，nD)

在二维条件下，每个三角形单元包含3个未知

速度分量，每条速度间断线上包含2个非负未知变量。在三维情况下，每个四面体单元则包含有6个未知速度分量，而每个速度间断面上包含4个非负未知变量。在施加流动法则、虚功方程以及边界条件等约束条件后，未知变量必须满足一系列的线性和非线性等式以及线性不等式，这样就构成了一个非线性数学规划问题。在本文中，采用序列2次规划法进行求解。

3 算例

根据以上讨论的方法，编制了能够进行二维和三维边坡稳定性分析的软件UBRFEM。本文采用了4个算例对本方法及其程序进行了验证。其中，对

于二维边坡的算例，各采用了3个不同密度的网格进行计算与比较。

3.1 无重量承受垂直荷载的二维边坡

图5为一无重量承受垂直荷载的二维边坡。材料特性为c=98 kPa， =30°，边坡的倾角χ=45°。根据滑移线理论可得极限荷载q=111.44 kPa。采用能量法，相应这个荷载对于16个条块的划分模式，搜索得到最小安全系数为F=1.006。

由表1可得，对于如图6(a)所示的粗糙网格，

F=1.034；对于如图6(b)所示的中等密度网格，F= 1.012；对于如图6(c)所示的细密网格，F=1.003，

十分接近理论值F=1.0。

图5 一无重量承受垂直荷载的二维边坡(算例1) Fig.5 Weightless slope with vertical surface load

(example 1)

以上结果表明，本方法的计算结果与网格剖分形式相关。随着单元的加密，计算精度也随之提高。 3.2 非均质二维边坡

第23卷第6期殷建华等. 岩土边坡稳定的刚体有限元上限分析法 903

表1 土层的力学参数

Table 1 Mechanical parameters of soil strata

土层号

粘聚力c / kPa

内摩擦角 / (°)

容重γ / kN·m

－3

为1.476和1.435，而采用图8(c)中的细密网格得到安全系数为F=1.418，与能量法的结果非常接近。

1 0.0 2 5.3 3 7.2

38.0 23.0 20.0

19.5 19.5 19.5

10 m

(a) 粗糙网格

10 m

(b) 中等密度网格

10 m

图6 刚体有限元网格划分(算例1) Fig.6 RFEM meshes (example 1)

如图7所示的算例2为澳大利亚ACAD协会考核稳定分析程序的一个考题，其边坡的力学参数见表1。本方法所得的结果与陈祖煜采用能量法得到的安全系数进行了比较，计算结果见表2。陈祖煜提供的安全系数为F=1.416；采用如图8(a)和8(b)所示的2个不同密度的网格所获得的安全系数分别

图7 一非均质二维边坡(算例2) Fig.7 Inhomogeneous slope (example 2)

表2 安全系数计算结果 Table 2 Results of safety factor

本文方法

其他方法

算例

粗糙网格中等密度网格细密网格

EMU

理论解1 1.034 1.012 1.003 1.0061.0002 1.476 1.435 1.418 1.416

╱ 4 1.436

╱

╱ 1.402

图8 刚体有限元网格划分

(算例2) Fig.8

RFEM meshes (example 2)

904 岩石力学与工程学报 2004年

该算例同样说明，刚体有限元的分析结果与网格划分有关。加密网格能够显著地改善计算结果。 3.3 对称楔形体三维边坡

本算例为一个如图9所示的对称楔形体，其几何、物理参数见表3。其中，楔形体坡高10.2 m，岩石容重26.46 kN/m3。考虑不同的材料抗剪强度，分别采用c为5，15和20 kPa，为15°和30°来进行计算。

图9 一对称楔形体几何形状示意图(算例3) Fig.9 Symmetrical wedge in geometry (example 3)

表3 对称楔形体的几何、物理参数(算例3) Table 3 Geometry and unit weight for symmetrical

wedge (example 3)

部位倾向/(°)

倾角/(°)

左结构面 120 65 右结构面 240 65 顶面 180 0 坡面

180 90

表4对2种分析方法计算的安全系数进行了比较。由表4可知，2种方法所得到的计算结果相接近，并且得出了一致的结论，即对于给定的粘聚力

c值，安全系数随着摩擦角的增大而增加；同样，对于给定的值，安全系数随着c的增大而增加。图10为刚体有限元离散图。

表4 对称楔形体的安全系数计算结果(算例3) Table 4 Results of safety factor for symmetrical

wedge (example 3)

c /kPa 方法

= 15°

= 30°

EMU 0.675 1.278 本文方法 0.668 1.272 10

EMU 0.835 1.430 本文方法 0.845 1.433

EMU 1.173 1.749 本文方法 1.173 1.755

图10 刚体有限元离散(算例3) Fig.10 RFEM discretisation (example 3)

3.4

球体滑面的三维边坡

图11为一球体滑面的三维边坡，物理力学参数

为c=49.8 kPa，容重γ=9.8 kN/m3

，其平面离散模式如图12所示。

由表2可知，此例文[20]给出的所谓闭合解为

F=1.402，本文方法得到的安全为F = 1.436，比闭合解稍大。

图11 一球体滑面的三维边坡(算例4) Fig.11 Spherical slip surface in purely cohesive soil

(example 4)

图12 平面离散模式(算例4)

Fig.12 Plan view of discretisation pattern (example 4)

4 结论

本文提出了一种利用刚体有限元进行二维和三

第23卷第6期殷建华等. 岩土边坡稳定的刚体有限元上限分析法 905

维边坡稳定性分析的上限方法，将求解边坡的最小安全系数问题转化为一个非线性数学规划问题，这种问题可以用序列2次规划法求解。

通过4个算例对本文提出的方法和相应的程序进行了验证，得到的数值解与其他方法得到的解进行了比较，说明了本方法的正确性和可靠性。相比于文[10～12]提出的应用有限单元法进行上限分析的方法，本方法具有更为简单有效的优点。致谢本文的研究过程中陈祖煜教授和Professor

Dave Chan给予了很好的建议，在此表示衷心的感

谢。

参考文献

Drucker D C，Prager W.

Soil mechanics and plastic analysis or limit design[J]. Q Appl. Math.，1952，10：157～165 2

Chen W F.

Limit Analysis and Soil Plasticity[M]. New York：Elsevier Scientific Publishing Co.，1975

3 Donald I，Chen Z Y. Slope stability analysis by the upper bound

approach：fundamentals and methods[J]. Can. Geotech. J.，1997，34(6)：853～862 4

Chen Z Y. The limit analysis for slopes：theory，methods and applications[A]. In：Proceedings of the International Symposium on Slope Stability Engineering[C].

Rollerdam：Balema，1999，15～29 5

Wang Y J，Yin J H，Lee C F. The influence of a non-associated flow

rule on the calculation of the factor of safety of soil slopes[J]. Int. J. Numer. Anal. Meth. Geomech.，2001，25：1 351～1 359 6

Wang Y J. Stability analysis of slopes and footings considering

different dilation angles of geomaterial[PhD Thesis][D]. Hong Kong：Department of Civil and Structural Engineering，The Hong Kong Polytechnic University，2001

7 Chen Z Y，Wang X G，Haberfield C，et al.

A three-dimensional slope

stability analysis method using the upper bound theorem (part I)：theory and methods[J].

Int. J. Rock Mech. Min. Sci.，2001，38：369～378

8 Chen Z Y，Wang J，Wang Y J，et al.

A three-dimensional slope

stability analysis method using the upper bound theorem (part II)：numerical approaches，applications and extensions[J]. Int. J. Rock Mech. Min. Sci.，2001，38：379～397

Sloan S W. Lower bound limit analysis using finite elements and linear programming[J]. Int. J. Numer Anal. Meth. Geomech.，1988，12：61～77

10 Sloan S W. Upper bound limit analysis using finite elements and

linear programming[J]. Int. J. Numer Anal. Meth. Geomech.，1989，13：263～282

11 Sloan S W，Kleeman P W. Upper bound limit analysis using

discontinuous velocity fields[J]. Comput Appl. Mech. Eng.，1995，127：293～314

12 Lyamin A V，Sloan S W. Upper bound limit analysis using linear

finite elements and non-linear programming[J]. Int. J. Numer. Anal. Meth. Geomech.，2002，26：181～216

13 Kim J，Salgado R，Yu H S.

Limit analysis of soil slopes subjected to

pore-water pressures[J].ASCE

J. Geotech. and Geoenvir. Eng.，1999，125(1)：49～58

14 Zhang X. Slope stability analysis based on the rigid finite element