CN110083883A - 基于极限分析下限定理的三维边坡安全系数求解方法 - Google Patents

Abstract

本发明公开了一种基于极限分析下限定理的三维边坡安全系数求解方法，该方法以块体界面作用力的大小、方向和作用点为主要变量，以满足下限定理许可静力场要求的静力平衡(同时满足力和力矩平衡)和屈服准则为约束条件，将边坡稳定分析问题转化为求解安全系数最大值的非线性优化问题；同时对发明公开的方法计算结果的准确性、可靠性通过一个代表性算例进行了验证，均取得了令人满意的结果。本发明为实际边坡工程的三维稳定性分析提供了新的手段。

Description

基于极限分析下限定理的三维边坡安全系数求解方法

技术领域

本发明涉及一种岩土工程边坡稳定性评价方法，具体地说，本发明涉及一种根据极限分析下限定理求解三维边坡安全系数的方法。本发明属于岩土工程边坡稳定性极限分析领域。

背景技术

目前，在边坡稳定性分析领域，因缺少成熟有效的三维分析方法，通常都将边坡稳定性分析简化为平面问题来考虑，而忽视了实际工程的空间效应。这样的做法在多数情况下是可以接受的，但对于某些规模巨大、空间效应十分显著的边坡，如仍采用二维分析的方法显然是不合适的。

国内外不少学者曾尝试将边坡二维极限平衡法的思想推广应用到三维边坡稳定性分析中，通过引入各种假定，使复杂的三维问题变得静定可解，相应不同的假定形成了各具特点的分析方法。与二维极限平衡法相比，三维边坡稳定性分析方法需引入大量假定，对于离散成n行m列条块的三维滑体,共需引入8×m×n个假定。这些假定的引入不仅大大削弱了这类方法的理论基础，而且这些假定的合理性在实际工程中很难得到满足，多数情况下无法得到令人满意的分析结果，甚至会导致错误的结论，这正是已有三维极限平衡法在实际工程没有得到推广应用的主要原因。

为了克服三维极限平衡法的不足，近年来已有不少学者探索将塑性力学极限分析理论与有限元结合，提出了边坡稳定性极限分析有限元法。该方法以塑性力学上下限定理为基础，虽然可以弥补传统的三维极限平衡法理论基础不够严密的缺陷，但同时因该方法需要对研究对象进行有限元离散，并以每一个单元节点的应力或位移作为未知量，不仅工作量大、计算效率低，而且由于未知量多、自由度高，给计算求解带来了很大的困难。现阶段相关研究成果还主要停留在一些相对简单的算例验证上，将其推广应用于复杂的实际工程中，特别是三维领域的边坡稳定性分析领域还将面临巨大的挑战。

因此，进一步深入研究三维边坡稳定性分析方法，减少不合理假定对稳定性分析结果的影响、提高分析方法适用范围和精度是非常必要的。

发明内容

鉴于上述原因，本发明的目的是提供一种基于极限分析下限定理的三维边坡安全系数求解方法。该方法以块体界面作用力的大小、方向和作用点为主要变量，以满足极限分析下限定理许可静力场要求的静力平衡(同时满足力和力矩平衡)和屈服准则为约束条件建立约束方程，以安全系数最大值为目标函数，形成三维边坡稳定分析下限解法的优化算法模型，将三维边坡稳定性分析问题转化为求解安全系数最大值的非线性优化问题；通过该优化算法模型求得三维边坡安全系数。

为实现上述目的，本发明采用以下技术方案：一种基于极限分析下限定理的三维边坡安全系数求解方法，它包括如下步骤：

S1：根据三维边坡的实际情况，确定该三维边坡的多边形界面，将三维边坡分割成若干个块体，通过多边形界面集合和块体集合定义该三维边坡；

S2：以块体界面作用力的大小、方向和作用点为变量体系，建立满足极限分析下限定理许可静力场要求的约束方程，该约束方程包括静力平衡方程、不违反屈服准则的约束不等式；

S3：以安全系数最大值为目标函数，以满足极限分析下限定理许可静力场要求的约束方程为约束条件，形成三维边坡稳定分析下限解法的优化算法模型；

S4：通过该优化算法模型求解三维边坡安全系数。

本发明提出的以界面PO_j为基础的变量体系，极大简化了约束方程未知量数量，约为以节点为基础的变量体系变形数量的1/10，为处理复杂三维边坡问题提供了基础。本发明根据空间多变形定义将在界面PO_j上屈服条件由“莫尔圆”转化为内接多边形的，使得原有四次方约束方程转化为二次方约束方程，极大降低了约束方程的非线性。本发明公开的基于极限分析下限定理的三维边坡安全系数求解方法，为实际边坡工程的三维稳定性分析提供了新的手段。

附图说明

图1是本发明基于极限分析下限定理的三维边坡安全系数求解方法流程图；

图2是本发明实例虚拟边坡模型；

图3是构建图2所示虚拟边坡模型的块体示意图；

图4是构建虚拟边坡模型的空间多边形界面及其边界的内法向量示意图；

图5是构建虚拟边坡模拟的块体共用边界上作用力与反作用力示意图；

图6是极限分析下限定理屈服准则线性化示意图；

图7是本发明根据极限分析下限定理优化求解三维边坡安全系数的方法流程图；

图8是利用本发明评价图2所示的虚拟边坡稳定性的实施例图。

具体实施方式

下面结合附图及实施例进一步阐明本发明，应理解这些实例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

目前在边坡稳定性分析领域，因缺少成熟有效的三维分析方法，通常都将边坡稳定性分析简化为平面问题来考虑，而忽视了实际工程的空间效应。也有学者尝试将边坡二维极限平衡法的思想推广应用到三维边坡稳定分析中，通过引入各种假定，使复杂的三维问题变得静定可解，相应不同的假定形成了各具特点的分析方法。这些假定的引入不仅大大削弱了这类方法的理论基础，而且这些假定的合理性在实际工程中很难得到满足，多数情况下无法得到令人满意的分析结果，甚至会导致错误的结论。

因此，本发明提供了一种基于极限分析下限定理的三维边坡安全系数求解方法，如图1所示，该方法包括如下步骤：

S1：根据三维边坡的实际情况，确定该三维边坡的多边形界面，将三维边坡分割成若干个块体，通过多边形界面集合和块体集合定义该三维边坡；

S2：以块体界面作用力的大小、方向和作用点为变量体系，建立满足极限分析下限定理许可静力场要求的约束方程，该约束方程包括静力平衡方程、不违反屈服准则的约束不等式；

S3：以安全系数最大值为目标函数，以满足极限分析下限定理许可静力场要求的约束方程为约束条件，形成三维边坡稳定分析下限解法的优化算法模型；

S4：通过该优化算法模型求解三维边坡安全系数。

本发明评价三维边坡稳定性的原理是：根据极限分析下限定理，以块体界面作用力的大小、方向和作用点为主要变量，以满足极限分析下限定理许可静力场要求的静力平衡(同时满足力和力矩平衡)和屈服准则为约束条件，以安全系数最大值为目标函数，形成三维边坡稳定分析下限解法的优化算法模型，将边坡稳定分析问题转化为求解安全系数最大值的非线性优化问题；通过优化算法求得三维边坡安全系数。

本发明步骤S1：根据三维边坡的实际情况，确定该三维边坡的多边形界面，将三维边坡分割成若干个块体，通过多边形界面集合和块体集合定义该三维边坡，具体方法如下：

S1.1：如图2、图3所示，根据三维边坡的实际情况，确定该三维边坡的多边形界面，将三维边坡分割成若干个块体，通过多边形界面集合和块体集合定义该三维边坡V。

该三维边坡V通过构成该三维边坡的多边形界面PO_j集合定义为：

V＝PO₁PO₂…PO_j…PO_p-1PO_p (1)

式中，j为边界面的编号，PO_j为编号为j的多边形界面，p为三维边坡V中多边形界面总数，PO_j的地层岩性、风化程度等特征将通过界面参数(黏聚力c_j和内摩擦角)反映。

该三维边坡V通过组成该三维边坡的块体B_k集合定义为：

V＝B₁B₂…B_k…B_q-1B_q (2)

式中，k为边界面的编号，B_k为编号为k的块体，q为三维边坡V中块体总数，B_k的地层岩性、风化程度等特征将通过块体参数(重度d_k)反映。

S1.2：空间多边形界面PO_j由一系列按顺序排列的顶点来描述，可用下式表示为：

式中i表示多边形界面PO_j中顶点的编号，为编号为j的界面上第i个顶点，h1为多边形PO_j的顶点总数；为表示方便，将多边形界面PO_j最后一个点与第一点设为重合以形成闭合环路，即

多边形界面PO_j的边界为：

式中为多边形PO_j的边界线段。

多边形PO_j的法线向量np_j可由式(5)求得：

当i＝1时

多边形PO_j的边界线段所对应的内法线向量为：

需要指出的是，通过上述公式计算得到的多边形边界线段内法线向量永远指向多边形内部，与顶点旋转方向无关。

根据上述法线向量定义，空间任意一点x位于空间多边形界面PO_j内的充要条件为，如图4所示：

式中，ap_j为多边形界面PO_j内任一定点。

本发明步骤S2：以块体界面PO_j作用力的大小、方向和作用点为变量体系，建立满足极限分析下限定理许可静力场要求的约束方程，包括静力平衡方程、不违反屈服准则的约束不等式，具体步骤如下：

S2.1：根据强度储备的概念，引入描述边坡强度的折减系数F，通过该变量描述边坡安全储备如下式所示；

c_e＝c/F (8)

式中c、为边坡材料实际凝聚力和内摩擦系数，c_e、f_e为边坡材料进行折减后的凝聚力和内摩擦系数，c_e、f_e是对应不同的折减系数F计算而得；

S2.2为描述三维边坡界面力的方向特性，在每个结构面上建立局部坐标系。

该三维边坡V任一界面PO_j(除边坡轮廓面外)的作用力包括法向力和位于该界面的切向力(可用两个分量表示)，其量值分别为N_j、Tm_j和Tl_j。

如图5所示，为描述界面力的方向特性，在每个空间多边形界面PO_j上建立临时局部坐标系，以该面的法线向量np_j为N‘轴正向；以空间多边形任一边界线段的内法线向量为Tm‘轴正向，i可取为多边形PO_j的任一边界线段，通常取i＝1；Tl‘轴正向则可根据N‘轴、Tm‘轴通过右手螺旋法则确定对于包含该界面的块体，还需进一步根据该界面在对应块体内的内法线向量与np_j的关系，最终确定该界面上局部坐标和轴的正向：

并有

式中i表示块体B_k中边界面的局部编号(区别与全局编号，块体B_k内从1～n的连续编号)，为在块体B_k中编号为i的多边形界面，为界面上作用力的三个分量；j为边界面的编号，PO_j为编号为j的多边形界面；和PO_j表示同一个界面，但局部坐标系可能不同。

通过这样的处理，每条界面上的界面力只需要一组变量N_j、Tl_j和Tm_j，就可以准确描述共用界面上作用力与反作用力的关系。

S2.3在三维空间中，多边形界面PO_j上界面力的作用点r_j可用其可空间坐标(r_xj,r_yj,r_zj)表示。根据式(7)定义，为保证作用点始终位于空间多边形PO_j内，必须满足如下条件：

条件1，作用点在空间多边形PO_j所在平面上：

(r_j-ap_j)·np_j＝0 (11)

条件2：作用点位于界面凸多边形区域内：

经上述步骤以块体界面作用力的大小、方向和作用点为变量体系，建立满足极限分析下限定理许可静力场要求的约束方程后，如图7所示，本发明步骤S3：以安全系数最大值为目标函数，以满足极限分析下限定理许可静力场要求的约束方程为约束条件，形成三维边坡稳定分析下限解法的优化算法模型，具体步骤如下：

S3.1：根据极限分析下限定理，建立目标函数；

边坡安全系数f：

f＝max(F) (13)

S3.2：在三维边坡中需要保证每一个块体力系平衡，力系平衡的充分必要条件是力系中各力在三个坐标轴上投影的代数和为零，且各力对三个坐标轴的力矩之和也分别为零。

具体到某一块体B_k而言，力的平衡方程可以表示为：

式中i表示块体B_k中边界面或外部荷载的局部编号(区别与全局编号，块体B_k内从1～n的连续编号)；为在块体B_k中编号为i的多边形界面， 为界面上作用力的三个分量，h2为块体B_k内界面总数；为块体B_k中编号为i的外部荷载，为已知量，h3为作用于B_k的外部荷载总数。lp为平衡方程投影方向的方向导数，

并有：

块体B_k的力矩平衡方程可以表示为：

式中，lm为旋转轴方向向量；为边界多边形作用点到力矩点M_k的力臂矢量；为外部荷载作用点到力矩点M_k的力臂矢量；RG_k＝RG_k-M_k为块体Bk重心RG_k到力矩点M_k的力臂矢量。

在实际应用中，M_k可取为块体重心PG_k，可进一步简化为：

S3.3：不违反屈服准则约束条件：每一个界面PO_j上的作用力应不违反摩尔库仑屈服准则，如下式所示：

c_j、分别为界面PO_j的实际凝聚力和内摩擦角；A_j为界面PO_j的面积。除边坡轮廓面外，上述不等式方程与各界面一一对应。式(18)的几何意义为，切向力T_j(Tm_j，Tl_j)在多边形PO_j上位于以界面力作用点r_j为圆心、半径为的圆内，如图6所示。

为降低屈服准则约束条件的非线性，将上述圆用其内接等边多边形代替，内接等边多边形第i条边对应的内法线向量可通过下式计算得到。

mlⁱ＝(-cosθⁱ,-sinθⁱ)(i＝1,2,…,h4) (19)

式中，h4为内接等边多边形边的条数，θⁱ为mlⁱ与Tm轴负方向的夹角，通过式(20)求得。

θⁱ＝(2i-1)π/h4 (20)

根据式(7)空间多边形界面点集的定义，式(18)可转化为：

(T_j-alⁱ)·mlⁱ≥0(i＝1,2,…,h4) (21)

即

T_j·mlⁱ-alⁱ·mlⁱ≥0(i＝1,2,…,h4) (22)

式中alⁱ为内接等边多边形界面第i条边上的任一定点。

根据内接多边形几何关系可知则式(21)可转化为：

S3.4：为反映岩土体不受拉的特性，采用了式(24)所示的约束条件：

N_j≥0 (24)

本发明步骤S4可以通过输出现有成熟的优化求解算法计算文件，也可以通过编程软件编制优化算法将运算步骤S3形成的优化算法模型，求解三维边坡安全系数f。

图8是利用本发明公开的一种基于极限定理的任意多面楔形体安全系数求解方法的实施例图。

通过C#编程软件，计算边坡安全系数，首先，输入构建三维边坡的各点、各线段、各界面多边形、各块体参数：然后，将推导出的优化算法表达式(1)-(24)编程输入至优化计算软件(如C、C++、C#等)求解安全系数。

图示的椭球形底滑面边坡是三维边坡稳定分析领域经常采用的算例，很多学者给出了相应于不同分析方法的稳定计算结果。计算结果表明，本发明计算得到的安全系数为2.174，介于极限平衡法(Zhang 1988，"Three-dimensional stability analysis ofconcave slopes in plan view."J GeotechEng 114(6):658–671.)和极限分析上限解法(Chen,Z.,H.Mi,F.Zhang and X.Wang(2003)."A simplified method for 3D slopestability analysis."Canadian Geotechnical Journal 40(3):675-683.)之间，比极限平衡法计算得到的安全系数大0.052，比极限分析上限解结果小0.088。如果考虑到三维极限平衡法所作的大量简化假定，本方法方法的计算结果理应更为合理。

本发明提出的以界面PO_j为基础的变量体系，极大简化了约束方程未知量数量，约为以节点为基础的变量体系变形数量的1/10，为处理复杂三维边坡问题提供了基础。本发明根据空间多变形定义将在界面PO_j上屈服条件由“莫尔圆”转化为内接多边形的，使得原有四次方约束方程转化为二次方约束方程，极大降低了约束方程的非线性。本发明公开的基于极限分析下限定理的三维边坡安全系数求解方法，为实际边坡工程的三维稳定性分析提供了新的手段。

Claims (4)

1.一种基于极限分析下限定理的三维边坡安全系数求解方法，其特征在于：它包括如下步骤：

S1：根据三维边坡的实际情况，确定该三维边坡的多边形界面，将三维边坡分割成若干个块体，通过多边形界面集合和块体集合定义该三维边坡；

S2：以块体界面作用力的大小、方向和作用点为变量体系，建立满足极限分析下限定理许可静力场要求的约束方程，该约束方程包括静力平衡方程、不违反屈服准则的约束不等式；

S3：以安全系数最大值为目标函数，以满足极限分析下限定理许可静力场要求的约束方程为约束条件，形成三维边坡稳定分析下限解法的优化算法模型；

S4：通过该优化算法模型求解三维边坡安全系数。

2.根据权利要求1所述的基于极限分析下限定理的三维边坡安全系数求解方法，其特征在于：所述步骤S1具体方法如下：

S1.1：根据三维边坡的实际情况，确定该三维边坡的多边形界面，将三维边坡分割成若干个块体，通过多边形界面集合和块体集合定义该三维边坡V；

该三维边坡V通过构成该三维边坡的多边形界面PO_j集合定义为：

V＝PO₁PO₂…PO_j…PO_p-1PO_p (1)

式中，j为边界面的编号，PO_j为编号为j的多边形界面，p为三维边坡V中多边形界面总数，PO_j的地层岩性、风化程度、结构面分布特征通过界面黏聚力c_j和内摩擦角反映；

该三维边坡V通过组成该三维边坡的块体B_k集合定义为：

V＝B₁B₂…B_k…B_q-1B_q (2)

式中，k为边界面的编号，B_k为编号为k的块体，q为三维边坡V中块体总数，B_k的地层岩性、风化程度等特征将通过块体参数反映；

S1.2：空间多边形界面PO_j由一系列按顺序排列的顶点来描述，可用下式表示为：

式中i表示多边形界面PO_j中顶点的编号，为编号为j的界面上第i个顶点，h1为多边形PO_j的顶点总数；为表示方便，将多边形界面PO_j最后一个点与第一点设为重合以形成闭合环路，即

多边形界面PO_j的边界为：

式中为多边形PO_j的边界线段；

多边形PO_j的法线向量np_j可由式(5)求得：

当i＝1时

多边形PO_j的边界线段所对应的内法线向量为：

需要指出的是，通过上述公式计算得到的多边形边界线段内法线向量永远指向多边形内部，与顶点旋转方向无关；

根据上述法线向量定义，三维边坡空间任意一点x位于空间多边形界面PO_j内的充要条件为：

式中，ap_j为多边形界面PO_j内任一定点。

3.根据权利要求2所述的基于极限分析下限定理的三维边坡安全系数求解方法，其特征在于：所述步骤S2以块体界面PO_j作用力的大小、方向和作用点为变量体系，建立满足极限分析下限定理许可静力场要求的约束方程，包括静力平衡方程、不违反屈服准则的约束不等式，具体步骤如下：

S2.1：根据强度储备的概念，引入描述边坡强度的折减系数F，通过该变量描述边坡安全储备如下式所示；

c_e＝c/F (8)

式中c、为边坡材料实际凝聚力和内摩擦系数，c_e、f_e为边坡材料进行折减后的凝聚力和内摩擦系数，c_e、f_e是对应不同的折减系数F计算而得；

S2.2为描述三维边坡界面力的方向特性，在每个结构面上建立局部坐标系；

该三维边坡V任一界面PO_j(除边坡轮廓面外)的作用力包括法向力和位于该界面的切向力(可用两个分量表示)，其量值分别为N_j、Tm_j和Tl_j；

为描述界面力的方向特性，在每个空间多边形界面PO_j上建立临时局部坐标系，以该面的法线向量np_j为N‘轴正向；以空间多边形任一边界线段的内法线向量为Tm‘轴正向，i可取为多边形PO_j的任一边界线段，通常取i＝1；Tl‘轴正向则可根据N‘轴、Tm‘轴通过右手螺旋法则确定对于包含该界面的块体，还需进一步根据该界面在对应块体内的内法线向量与np_j的关系，最终确定该界面上局部坐标和轴的正向：

并有

式中i表示块体B_k中边界面的局部编号；为在块体B_k中编号为i的多边形界面，为界面上作用力的三个分量；j为边界面的编号，PO_j为编号为j的多边形界面；S2.3在三维空间中，多边形界面PO_j上界面力的作用点r_j可用其可空间坐标(r_xj,r_yj,r_zj)表示，根据式(7)定义，为保证作用点始终位于空间多边形PO_j内，必须满足如下条件：

条件1，作用点在空间多边形PO_j所在平面上：

(r_j-ap_j)·np_j＝0 (11)

条件2：作用点位于界面凸多边形区域内：

4.根据权利要求3所述的基于极限分析下限定理的三维边坡安全系数求解方法，其特征在于：所述步骤S3以安全系数最大值为目标函数，以满足极限分析下限定理许可静力场要求的约束方程为约束条件，形成三维边坡稳定分析下限解法的优化算法模型，具体步骤如下：

S3.1：根据极限分析下限定理，建立目标函数；

边坡安全系数f：

f＝max(F) (13)

S3.2：在三维边坡中需要保证每一个块体力系平衡，力系平衡的充分必要条件是力系中各力在三个坐标轴上投影的代数和为零，且各力对三个坐标轴的力矩之和也分别为零；

具体到某一块体B_k而言，力的平衡方程可以表示为：

式中i表示块体B_k中边界面或外部荷载的局部编号(区别与全局编号，块体B_k内从1～n的连续编号)；为在块体B_k中编号为i的多边形界面， 为界面上作用力的三个分量，h2为块体B_k内界面总数；为块体B_k中编号为i的外部荷载，为已知量，h3为作用于B_k的外部荷载总数，l_p为平衡方程投影方向的方向导数，

并有：

块体B_k的力矩平衡方程可以表示为：

式中，lm为旋转轴方向向量；为边界多边形作用点到力矩点M_k的力臂矢量；为外部荷载作用点到力矩点M_k的力臂矢量；RG_k＝PG_k-M_k为块体Bk重心PG_k到力矩点M_k的力臂矢量，

在实际应用中，M_k可取为块体重心PG_k，可进一步简化为：

S3.3：不违反屈服准则约束条件：每一个界面PO_j上的作用力应不违反摩尔库仑屈服准则，如下式所示：

c_j、分别为界面PO_j的实际凝聚力和内摩擦角；A_j为界面PO_j的面积，

式(18)的几何意义为，切向力T_j(Tm_j，Tl_j)在多边形PO_j上位于以界面力作用点r_j为圆心、半径为的圆内；

为降低屈服准则约束条件的非线性，将上述圆用其内接等边多边形代替，内接等边多边形第i条边对应的内法线向量可通过下式计算得到：

mlⁱ＝(-cosθⁱ,-sinθⁱ)(i＝1,2,…,h4) (19)

式中，h4为内接等边多边形边的条数，θⁱ为mlⁱ与Tm轴负方向的夹角，通过式(20)求得：

θⁱ＝(2i-1)π/h4 (20)

根据式(7)空间多边形界面点集的定义，式(18)可转化为：

(T_j-alⁱ)·mlⁱ≥0(i＝1,2,…,h4) (21)

即

T_j·mlⁱ-alⁱ·mlⁱ≥0(i＝1,2,…,h4) (22)

式中alⁱ为内接等边多边形界面第i条边上的任一定点；

根据内接多边形几何关系可知则式(21)可转化为：

S3.4：为反映岩土体不受拉的特性，采用了式(24)所示的约束条件：

N_j≥0 (24)。