CN112214921B - 一种基于应力矢量的三维边坡稳定性评价方法 - Google Patents

Abstract

本发明涉及一种基于应力矢量的三维边坡稳定性评价方法，首先，将三维边坡宽、长和高所围成的最小立方体离散成若干个立方体搜索网格，并进行编号，各搜索网格内包含一个以上数值计算单元，所有的立方体搜索网格组成了多个条柱；然后，在各搜索网格内确定唯一对应的数值计算单元作为目标数值计算单元；接着，根据各条柱中各搜索网格对应的目标数值计算单元确定各条柱的潜在破坏单元；最后，由所有的潜在破坏单元得到最危险滑移面，同时计算计算三维边坡整体安全系数和局部区域安全系数。本发明可分析复杂条件下的三维边坡稳定性，无需对滑移面形状及滑入、滑出点限制，可求解极限与非极限状态边坡稳定性，可以对复杂应力状况下边坡进行滑移面搜索。

Description

一种基于应力矢量的三维边坡稳定性评价方法

技术领域

本发明属于岩土体边坡稳定性分析及防灾减灾技术领域，涉及一种基于应力矢量的三维边坡稳定性评价方法。

背景技术

边坡稳定性分析及其加固治理问题是一个历史悠久且又充满活力的科学分支，人类与滑坡灾害作斗争的努力一直没有中断过。从1970年开始，对于边坡的三维稳定性研究渐渐展开，现在所形成的三维分析方法多源自二维稳定性分析方法，以三维极限平衡法和有限元强度折减法等传统方法为主要研究手段，但现阶段的三维稳定性分析方法依然存在缺陷，如三维极限平衡法对土体条间力等进行了假定且将滑移面假设为椭球面，与实际情况相差较大。这些方法都存在诸多的假定，尤其是假设了滑移面的形状，极大限制了极限平衡法的应用范围。从实际情况看，滑移面形状可以为任意形状，且边坡土体内应力情况复杂。然而，随着计算机技术的发展，国内外在进行边坡稳定性分析时逐渐开始使用有限元法与计算机相结合的方法，能更好地对滑移面进行确定，进而对边坡稳定性有了更加准确的分析。如在用有限元分析建立应力场的基础上，采用蚂蚁算法搜索出更接近实际情况的滑移面从而对边坡稳定性进行分析，但这种算法有着收敛慢、容易陷入局部最优的缺点。除此以外还有在已知应力场的情况下，使用多种群遗传算法构建适应度函数搜索最危险滑移面，达到边坡稳定性分析的目的；但这类算法有着编程过程复杂，交叉率和变异率等参数通常需要靠经验选择的缺点，不易推广应用。

发明内容

本发明的目的是解决现有技术中存在的问题，本发明提出了一种基于土体应力应变状态下的三维边坡空间滑移面搜索方法，进而建立了基于应力矢量的三维边坡稳定性评价方法。

本发明提供一种基于离散应力应变场的三维边坡滑移面条柱网格式搜索方法，其属于岩土体边坡稳定性分析及防灾减灾技术领域。本发明重新构建条柱式搜索网格，离散三维边坡应力应变场，提出了靶单元(潜在破坏单元)判别准则，创建了一种条柱式网格化搜索模式求解三维边坡空间潜在危险滑移面。该方法根据边坡几何参数建立三维边坡空间坐标系统，将三维边坡的应力应变状态离散为若干个空间单元，依据靶单元的空间坐标及其应力应变数值，可精准生成三维边坡最危险滑移面。本发明考虑三维空间边坡下滑方向具有矢量性，提出任意空间条柱内离散滑移面上的法向量和下滑方向矢量求解方法，从而建立基于应力矢量的三维边坡稳定性分析方法。

为达到上述目的，本发明采用的技术方案如下：

一种基于应力矢量的三维边坡稳定性评价方法，步骤如下：

(1)建立搜索网格；

通过数值分析软件建立边坡三维模型并计算，提取三维模型几何信息和应力应变计算结果，导入按上述算法编制的程序或软件中，设置搜索网格尺寸参数及网格坐标生成方式，建立直角坐标系下三维边坡滑移面搜索网格和模型，具体过程如下：

建立x轴平行于三维边坡宽、y轴平行于三维边坡长且z轴平行于三维边坡高的三维直角坐标系x-y-z，将三维边坡宽、长和高所围成的最小立方体离散成n×m×k个立方体搜索网格，对每个搜索网格进行编号，即完成搜索网格的建立；

n×m×k个立方体搜索网格组成了n×m个条柱，各条柱的底面位于三维边坡宽和长所在平面上，同时各条柱的高等于三维边坡高，m为所有的条柱排列形成的平行于x轴的行数，n为所有的条柱排列形成的平行于y轴的列数，k为各条柱中立方体搜索网格的数量，n、m、k的取值满足：各搜索网格内包含一个以上数值计算单元；

(2)搜索目标数值计算单元；

设搜索网格i的几何中心为Z点，Z点的空间坐标为(x_iz,y_iz,z_iz)；搜索网格i的八个顶点为V_i1、V_i2、V_i3、V_i4、V_i5、V_i6、V_i7、V_i8，八个顶点的空间坐标分别对应为(x_i1,y_i1,z_i1)、(x_i2,y_i2,z_i2)、……、(x_i8,y_i8,z_i8)；搜索网格i的数值计算单元j的四个顶点为A、B、C、D，四个顶点的空间坐标分别对应为(x_ja,y_ja,z_ja)、(x_jb,y_jb,z_jb)、(x_jc,y_jc,z_jc)、(x_jd,y_jd,z_jd)；搜索同时满足以下四个条件的数值计算单元，即为搜索网格i的目标数值计算单元；

①x_{i_min}<x_ja、x_jb、x_jc、x_jd<x_{i_max}；

②y_{i_min}<y_ja、y_jb、y_jc、y_jd<y_{i_max}；

③z_{i_min}<z_ja、z_jb、z_jc、z_jd<z_{i_max}；

④v_{j_abcd}＝v_{j_bcdz}+v_{j_cadz}+v_{j_dabz}+v_{j_bacz}；

其中，x_{i_min}、x_{i_max}、y_{i_min}、y_{i_max}、z_{i_min}、z_{i_max}分别为搜索网格i的八个顶点的空间坐标中的最小x坐标、最大x坐标、最小y坐标、最大y坐标、最小z坐标、最大z坐标；v_{j_abcd}表示搜索网格i的数值计算单元j的体积；v_{j_bcdz}表示由Z点和数值计算单元j的B、C、D顶点围成的四面体体积；v_{j_cadz}表示由Z点和数值计算单元j的C、A、D顶点围成的四面体体积；v_{j_dabz}表示由Z点和数值计算单元j的D、A、B顶点围成的四面体体积；v_{j_bacz}表示由Z点和数值计算单元j的B、A、C顶点围成的四面体体积；

按此方法，依次确定n×m×k个立方体搜索网格对应的目标数值计算单元；

若搜索网格内不包含目标数值计算单元，则该搜索网格为无效网格，在计算过程中按零标记，不参与计算过程；

(3)搜索潜在破坏单元；

搜索每个条柱的k个搜索网格对应的目标数值计算单元中同时满足以下两个条件的目标数值计算单元，该目标数值计算单元对应的搜索网格即为该条柱的潜在破坏单元；

①f≤[f]，

②ε_p≥[ε_p]；

其中，f表示目标数值计算单元的安全系数；[f]表示安全系数允许值；C为土体的粘聚力；σ_n表示目标数值计算单元所在位置处任意滑移面上法线方向正应力；

为土体的摩擦角；τ_d表示目标数值计算单元所在位置处任意滑移面上下滑方向剪应力；ε_p表示目标数值计算单元对应的塑性应变；[ε_p]表示滑移面搜索时允许的塑性应变值，取值大小根据边坡安全控制要求确定，取值范围介于0～1％ε_pmax之间，ε_pmax表示边坡最大塑性应变值；

σ_n和τ_d通过以下公式计算得到：

其中，F_N′表示F_N的转置；F_N表示任意方向平面上正应力矢量；N表示该任意方向平面的法线方向；

表示该任意方向平面的法向量；F_D′表示F_D的转置；F_D表示任意方向平面上剪应力矢量；D表示该任意方向平面的下滑方向；

表示该任意方向平面的下滑方向矢量；σ_pq表示目标数值计算单元对应的应力张量(应力状态的数学表示)，p＝1，2，3，q＝1，2，3，数学上应力为二阶张量，三维空间中需九个分量来确定，σ₁₁、σ₂₂、σ₃₃、τ₂₁、τ₁₂、τ₃₁、τ₁₃、τ₃₂、τ₂₃为应力张量的9个分量；

和

的确定流程如下：

(I)生成滑移面；

将搜索网格i的12条边进行分段，以所有边上分段后形成的断点为集合，取其中不共线的任意三点S_i1、S_i2、S_i3为组合，由S_i1、S_i2、S_i3点形成的面即为搜索网格i的滑移面；

(II)判断S_i1、S_i2、S_i3点是否同时位于搜索网格i所在条柱z轴方向的4条棱上，如果是，则将S_i1、S_i2、S_i3点作为H、I、J点进入步骤(IV)；反之，则进入步骤(III)；

(III)确定S_i1、S_i2、S_i3点所在的平面与搜索网格i所在条柱z轴方向的4条棱的四个交点H、I、J、K；

(IV)计算

和

的计算公式如下：

其中，

表示面HIJ的法向量；

表示点H与点I组成的线段的向量；

表示点H与点J组成的线段的向量；x_H、y_H、z_H分别对应为H点的x坐标、y坐标、z坐标；x_I、y_I、z_I分别对应为I点的x坐标、y坐标、z坐标；x_J、y_J、z_J分别对应为J点的x坐标、y坐标、z坐标；

的计算公式如下：

①当H、I、J点的z坐标相等时，

②当只有H、J点的z坐标相等时，

③当只有H、I点的z坐标相等时，

④当H、I、J点的z坐标均不相等时，

如果存在多个目标数值计算单元满足以上两个条件，则将f值最小的目标数值计算单元对应的搜索网格作为潜在破坏单元；如果不存在目标数值计算单元满足以上两个条件，则以包含目标数值计算单元且最靠近该条柱顶部的搜索网格作为潜在破坏单元；

按此方法，依次确定n×m个条柱的潜在破坏单元；

(4)将n×m个条柱的潜在破坏单元的中心点空间坐标拟合生成空间曲面，形成三维边坡的整体空间最危险滑移面，同时按以下公式计算三维边坡的整体安全系数和局部区域安全系数，并将其与参考值进行比较，定量评价边坡稳定性，参考值为与边坡相关的国家及地方规范规定的允许安全系数值；

式中，K_W表示三维边坡的整体安全系数；d表示潜在破坏单元的总数；k_v表示第v个潜在破坏单元的安全系数；σ_nv表示第v个潜在破坏单元的滑移面上法线方向正应力；τ_dv表示第v个潜在破坏单元的滑移面上下滑方向剪应力；C为土体的粘聚力；

为土体的摩擦角；

式中，K_L表示三维边坡的局部区域安全系数；e表示需要求解的局部区域所包含的潜在破坏单元的总数；k_g表示需要求解的局部区域所包含的第g个潜在破坏单元的安全系数；σ_ng表示需要求解的局部区域所包含的第g个潜在破坏单元的滑移面上法线方向正应力；τ_dg表示需要求解的局部区域所包含的第g个潜在破坏单元的滑移面上下滑方向剪应力；C为土体的粘聚力；

为土体的摩擦角。

作为优选的技术方案：

如上所述的一种基于应力矢量的三维边坡稳定性评价方法，步骤(1)具体如下：

(1.1)根据三维边坡的几何尺寸建立三维直角坐标系x-y-z，三维边坡宽平行于x轴，三维边坡长平行于y轴，三维边坡高平行于z轴；

(1.2)将三维边坡宽、长和高所围成的最小立方体划分为沿x轴依次排列的n个大立方体；

(1.3)将每个大立方体划分为沿y轴依次排列的m个中立方体；

(1.4)将每个中立方体划分为沿z轴依次排列的k个小立方体，至此，三维边坡宽、长和高所围成的最小立方体离散成n×m×k个立方体搜索网格；

(1.5)对每个搜索网格进行编号，即完成搜索网格的建立。

如上所述的一种基于应力矢量的三维边坡稳定性评价方法，所有的划分都为均匀划分。

如上所述的一种基于应力矢量的三维边坡稳定性评价方法，步骤(1)中，各搜索网格内数值计算单元的个数为1～125。

如上所述的一种基于应力矢量的三维边坡稳定性评价方法，步骤(3)中，搜索每个条柱的k个搜索网格对应的目标数值计算单元中同时满足以下两个条件的目标数值计算单元时，按从下往上或从上往下的顺序。

如上所述的一种基于应力矢量的三维边坡稳定性评价方法，步骤(3)中，将搜索网格i的12条边进行分段是指均匀分段。

如上所述的一种基于应力矢量的三维边坡稳定性评价方法，步骤(3)中，搜索网格i的每条边被分为3～10段。

如上所述的一种基于应力矢量的三维边坡稳定性评价方法，步骤(4)中，与边坡相关的国家及地方规范包括《建筑边坡工程技术规范》GB50330-2013和《建筑基坑支护技术规程》JGJ120-2012。

有益效果：

(1)本发明适用于复杂条件下的三维边坡稳定性分析，无需对三维边坡滑移面形状及滑入、滑出点限制，可用求解极限与非极限状态边坡稳定性，可以对复杂应力状况下边坡进行滑移面搜索；

(2)本发明适用计算机编程，算法简捷，计算效率高，便于基于商业数值分析软件进行三维边坡稳定性分析技术的二次开发；

(3)本发明的方法所得空间滑移面位置和安全系数空间分布更具有合理性和准确性，更能体现出三维边坡空间效应，优势明显，它是一项解决三维边坡稳定性分析、控制及评价的核心技术。

附图说明

图1为搜索网格的示意图；

图2为搜索网格i与数值计算单元j的对应关系示意图；

图3为在搜索网格i内确定滑移面的示意图；

图4为搜索网格内滑移面与所在条柱z轴方向棱相交示意图；

图5为任意条柱内滑移面的法向量及下滑方向矢量求解示意图；

图6为三维基坑几何模型平面图(单位：m)；

图7为三维基坑几何模型的a-a断面图(单位：m)；

图8为三维基坑几何模型的b-b断面图(单位：m)；

图9为三维边坡数值计算模型示意图；

图10为位移等高线云图；

图11为第95根条柱潜在破坏单元搜索过程结果图解显示；

图12为第140根条柱潜在破坏单元搜索过程结果图解显示；

图13为基坑三维空间滑移面及安全系数空间分布图；

图14为基坑三维空间局部滑移面及安全系数分布图；

图15为图6的1-1断面处二维滑移面位置及安全系数分布；

图16为图6的2-2断面处二维滑移面位置及安全系数分布；

图17为图6的3-3断面处二维滑移面位置及安全系数分布。

具体实施方式

下面结合具体实施方式，进一步阐述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

对于三维边坡，通常受到地形、周边环境等影响，存在空间效应问题。三维空间内的滑移面会呈现多种可能的形状和特征。为实现三维空间滑移面的搜索，本发明为解决当前三维边坡稳定性分析滑移面搜索时引入过多假设及不考虑土体空间效应的弊端，基于边坡应力应变场提出一种三维边坡滑移面搜索方法及其应力矢量的稳定性评价方法。

一种基于应力矢量的三维边坡稳定性评价方法，步骤如下：

(1)建立搜索网格；

建立x轴平行于三维边坡宽、y轴平行于三维边坡长且z轴平行于三维边坡高的三维直角坐标系x-y-z，将三维边坡宽、长和高所围成的最小立方体离散成n×m×k个大小相等的如图1所示的立方体搜索网格，对每个搜索网格进行编号，以实现搜索过程控制，即完成搜索网格的建立；n×m×k个立方体搜索网格组成了n×m个条柱，各条柱的底面位于三维边坡宽和长所在平面上，同时各条柱的高等于三维边坡高，m为所有的条柱排列形成的平行于x轴的行数，n为所有的条柱排列形成的平行于y轴的列数，k为各条柱中立方体搜索网格的数量，n、m、k的取值满足：各搜索网格内包含一个以上数值计算单元(优选包含1～125个数值计算单元)，此处数值计算单元指建立数值模型时划分的计算单元；具体步骤如下：

(1.1)根据三维边坡的几何尺寸建立三维直角坐标系x-y-z，三维边坡宽平行于x轴，三维边坡长平行于y轴，三维边坡高平行于z轴；

(1.2)将三维边坡宽、长和高所围成的最小立方体均匀划分为沿x轴依次排列的n个大立方体；

(1.3)将每个大立方体均匀划分为沿y轴依次排列的m个中立方体；

(1.4)将每个中立方体均匀划分为沿z轴依次排列的k个小立方体，至此，三维边坡宽、长和高所围成的最小立方体离散成n×m×k个大小相等的立方体搜索网格；

(1.5)对每个搜索网格进行编号，即完成搜索网格的建立；

(2)搜索目标数值计算单元；

如图2所示，设搜索网格i的几何中心为Z点，Z点的空间坐标为(x_iz,y_iz,z_iz)；搜索网格i的八个顶点为V_i1、V_i2、V_i3、V_i4、V_i5、V_i6、V_i7、V_i8，八个顶点的空间坐标分别对应为(x_i1,y_i1,z_i1)、(x_i2,y_i2,z_i2)、……、(x_i8,y_i8,z_i8)；搜索网格i的数值计算单元j的四个顶点为A、B、C、D，四个顶点的空间坐标分别对应为(x_ja,y_ja,z_ja)、(x_jb,y_jb,z_jb)、(x_jc,y_jc,z_jc)、(x_jd,y_jd,z_jd)；建立搜索网格i的搜索范围矩阵S_i，其中S_i＝[x_{i_min},x_{i_max},y_{i_min},y_{i_max},z_{i_min},z_{i_max}]，x_{i_min}、x_{i_max}、y_{i_min}、y_{i_max}、z_{i_min}、z_{i_max}分别为搜索网格i的八个顶点的空间坐标中的最小x坐标、最大x坐标、最小y坐标、最大y坐标、最小z坐标、最大z坐标，搜索同时满足以下四个条件的数值计算单元，即为搜索网格i的目标数值计算单元；

①x_{i_min}<x_ja、x_jb、x_jc、x_jd<x_{i_max}；

②y_{i_min}<y_ja、y_jb、y_jc、y_jd<y_{i_max}；

③z_{i_min}<z_ja、z_jb、z_jc、z_jd<z_{i_max}；

④v_{j_abcd}＝v_{j_bcdz}+v_{j_cadz}+v_{j_dabz}+v_{j_bacz}；

其中，x_{i_min}、x_{i_max}、y_{i_min}、y_{i_max}、z_{i_min}、z_{i_max}分别为搜索网格i的八个顶点的空间坐标中的最小x坐标、最大x坐标、最小y坐标、最大y坐标、最小z坐标、最大z坐标；v_{j_abcd}表示搜索网格i的数值计算单元j的体积；v_{j_bcdz}表示由Z点和数值计算单元j的B、C、D顶点围成的四面体体积；v_{j_cadz}表示由Z点和数值计算单元j的C、A、D顶点围成的四面体体积；v_{j_dabz}表示由Z点和数值计算单元j的D、A、B顶点围成的四面体体积；v_{j_bacz}表示由Z点和数值计算单元j的B、A、C顶点围成的四面体体积；(注：以上四面体体积均为标量，当运用矢量计算体积时，若体积为负值，须取绝对值)

按此方法，依次确定n×m×k个立方体搜索网格对应的目标数值计算单元；

若搜索网格内不包含目标数值计算单元，则该搜索网格为无效网格，在计算过程中按零标记，不参与计算过程；

(3)搜索潜在破坏单元；

按从下往上或从上往下的顺序搜索每个条柱的k个搜索网格对应的目标数值计算单元中同时满足以下两个条件的目标数值计算单元，该目标数值计算单元对应的搜索网格即为该条柱的潜在破坏单元；

①f≤[f]，

②ε_p≥[ε_p]；

其中，f表示目标数值计算单元的安全系数；[f]表示安全系数允许值，可根据相关规范取值(安全系数允许值为国家规范给定的数值)；C为土体的粘聚力；σ_n表示目标数值计算单元所在位置处任意滑移面上法线方向正应力；

为土体的摩擦角；τ_d表示目标数值计算单元所在位置处任意滑移面上下滑方向剪应力；ε_p表示目标数值计算单元对应的塑性应变；[ε_p]表示滑移面搜索时允许的塑性应变值，取值大小根据边坡安全控制要求确定，取值范围介于0～1％ε_pmax之间，ε_pmax表示边坡最大塑性应变值；

ε_p、ε_pmax均为事先经数值计算求得的应变(或现场监测取得)；C、

为土体固有物理参数，通过试验取得；

σ_n和τ_d通过以下公式计算得到：

其中，F_N′表示F_N的转置；F_N表示任意方向平面上正应力矢量；N表示该任意方向平面的法线方向；

表示该任意方向平面的法向量；F_D′表示F_D的转置；F_D表示任意方向平面上剪应力矢量；D表示该任意方向平面的下滑方向；

表示该任意方向平面的下滑方向矢量；σ_pq表示目标数值计算单元对应的应力张量(应力状态的数学表示)，p＝1，2，3，q＝1，2，3，数学上应力为二阶张量，三维空间中需九个分量来确定，σ₁₁、σ₂₂、σ₃₃、τ₂₁、τ₁₂、τ₃₁、τ₁₃、τ₃₂、τ₂₃为应力张量的9个分量；

在任意空间条柱内，均存在整体滑移面与立柱相交生成的立柱范围内局部滑移面，当条柱网格划分较小时，可假定网格内滑移面近似平面，用H、I、J、K四点分别表示滑移面与条柱z轴方向的4条棱的交点，其中任意三个点即可确定

和

以下取H、I、J三点所确定的平面为例说明条柱内滑移面的法向量和下滑方向矢量的求解过程；

和

的确定流程如下：

(I)生成滑移面；

搜索网格内不同的滑移面，对应求得不同的安全系数，每个搜索网格内任意不共线的三点均可确定一个滑动面(平面)，从理论上讲，这样的滑移面可以有无数个，为实现最危险滑移面的搜索，本发明给出一种立方体搜索网格内任意滑移面形成方法，具体为：将搜索网格i的12条边进行均匀分段(搜索网格i的每条边被分为3～10段，图3示出了搜索网格i的每条边被分为4段的情况)，以所有边上分段后形成的断点为集合，取其中不共线的任意三点S_i1、S_i2、S_i3为组合，如图3所示，由S_i1、S_i2、S_i3点形成的面即为搜索网格i的滑移面，每条边分段数越多，形成的点越多，生成的滑移面也就越多，计算精度也相对更为准确；

按照以上滑移面形成方法，每一个搜索网格均可生成若干个滑移面，这些滑移面可按上述方法分别求出对应的安全系数，在计算求解过程中，记录安全系数最小值及其对应的滑移面坐标，即可求解出任意搜索网格i的最小安全系数及最危险滑移面；

(II)判断S_i1、S_i2、S_i3点是否同时位于搜索网格i所在条柱z轴方向的4条棱上，如果是，则将S_i1、S_i2、S_i3点作为H、I、J点进入步骤(IV)；反之，则进入步骤(III)；

(III)如图4所示，确定S_i1、S_i2、S_i3点所在的平面与搜索网格i所在条柱z轴方向的4条棱的四个交点H、I、J、K；

(IV)计算

和

的计算公式如下：

其中，

表示面HIJ的法向量；

表示点H与点I组成的线段的向量；

表示点H与点J组成的线段的向量；x_H、y_H、z_H分别对应为H点的x坐标、y坐标、z坐标；x_I、y_I、z_I分别对应为I点的x坐标、y坐标、z坐标；x_J、y_J、z_J分别对应为J点的x坐标、y坐标、z坐标；

的计算公式如下：

①当H、I、J点的z坐标相等时，如图5所示，此时面HIJ为垂直于z轴的平面，在自重应力作用下，此时认为下滑方向矢量为零，即

②当只有H、J点的z坐标相等时，如图5所示，此时面HIJ仅为沿x轴方向具有坡度的斜平面，此时的下滑方向矢量

与向量

一致，即

③当只有H、I点的z坐标相等时，如图5所示，此时面HIJ仅为沿y轴方向具有坡度的斜平面，此时的下滑方向矢量

与向量

一致，即

④当H、I、J点的z坐标均不相等时，如图5所示，此时面HIJ是为沿x轴、y轴方向同时具有坡度的斜平面，下滑方向矢量则为②+③两种情况的叠加，根据矢量和运算，

即为条柱内任意面HIJ的下滑方向矢量，即

具体推导过程如下：

如图5所示，首先，在平面HIJK内，过J点做垂直于HI的线段JI’，再过I点沿y轴做水平线垂直相交于条柱棱上点J’处，过H点沿y轴做水平线垂直相交于条柱线上点K’，此时构成面HIJ’K’，其特征符合以上第②条，其下滑方向矢量则为

在面HIJK内，当以面HIJ’K’为参考平面时，由JI’垂直于HI，则面HIJK的特征则符合以上第③条，面HIJK下滑方向矢量则为

(等于

)；

根据矢量叠加原理，④＝②+③，面HIJK在条柱内的下滑方向矢量：

如上所述，根据以上各向量间几何关系，在已知HIJK四点坐标的情况下，通过向量运算，即求得面HIJK的法向量及下滑方向矢量；

如果存在多个目标数值计算单元满足以上两个条件，则将f值最小的目标数值计算单元对应的搜索网格作为潜在破坏单元；如果不存在目标数值计算单元满足以上两个条件，则以包含目标数值计算单元且最靠近该条柱顶部的搜索网格作为潜在破坏单元；

按此方法，依次确定n×m个条柱的潜在破坏单元；

(4)将n×m个条柱的潜在破坏单元的中心点空间坐标拟合生成空间曲面，形成三维边坡的整体空间最危险滑移面(也可以根据需求设定局部区域坐标范围，生成局部空间滑移面)，同时按以下公式计算三维边坡的整体安全系数和局部区域安全系数，并将其与参考值进行比较，定量评价边坡稳定性，参考值为与边坡相关的国家及地方规范(包括《建筑边坡工程技术规范》GB50330-2013和《建筑基坑支护技术规程》JGJ120-2012)规定的允许安全系数值；

式中，K_W表示三维边坡的整体安全系数；d表示潜在破坏单元的总数；k_v表示第v个潜在破坏单元的安全系数；σ_nv表示第v个潜在破坏单元的滑移面上法线方向正应力；τ_dv表示第v个潜在破坏单元的滑移面上下滑方向剪应力；C为土体的粘聚力；

为土体的摩擦角；

式中，K_L表示三维边坡的局部区域安全系数；e表示需要求解的局部区域所包含的潜在破坏单元的总数；k_g表示需要求解的局部区域所包含的第g个潜在破坏单元的安全系数；σ_ng表示需要求解的局部区域所包含的第g个潜在破坏单元的滑移面上法线方向正应力；τ_dg表示需要求解的局部区域所包含的第g个潜在破坏单元的滑移面上下滑方向剪应力；C为土体的粘聚力；

为土体的摩擦角；

此外，需注意的是，上述已经解释过潜在破坏单元为特定的目标数值计算单元对应的搜索网格，搜索网格的应力和应变等于搜索网格中心点对应的目标数值计算单元的应力和应变，本方法认为每个搜索网格内所有点的应力、应变值均相同，所以用搜索网格中心点应力、应变表示，所以说当网格划分越细时，精度就更高。

本发明的方法考虑了三维空间边坡下滑方向的矢量性，建立基于应力矢量的三维边坡稳定性分析方法，无需对三维边坡滑移面形状及滑入、滑出点限制，可用求解极限与非极限状态边坡稳定性，可以对复杂应力状况下边坡进行滑移面搜索。该方法除得到整体安全系数值以外，还可以得到局部区域安全系数值，可为边坡分区开挖及局部加固设计提供科学的计算依据，相对于传统方法更能体现出三维边坡空间效应，有更强的适用性和计算优势；

现结合具体案例进行说明，具体如下：

(i)建立三维基坑几何模型，其平面图如图6所示，其a-a断面如图7所示，其b-b断面如图8所示；

(ii)计算参数与模型；

基坑土体为均质土层，本构模型采用Mohr Coulomb Plasticity，土层参数如表1所示：

表1

(以上表中的参数都是数值计算需要的，C、

是在本方法中也用到的，E、γ、μ、Ψ只有在数值计算时用到，相当于前期计算，本方法是后期计算)

该模型采用Abaqus有限元软件计算，荷载仅考虑重力荷载，边界条件为基坑模型底面约束Z轴方向位移为零，基坑模型侧面约束与其垂直方向(X轴或Y轴)位移为零，模型单元形状采用四面体，单元类型为三维实体单元C3D10M:A 10-node modified quadratictetrahedron，方程求解器采用直接法和非对称矩阵存储，该有限元网格自动划分为152168个单元，如图9所示；

如图9所示，利用商业有限元软件建立的三维边坡数值计算模型，荷载只考虑自重，其边界约束条件为：上边界面自由无约束条件，前后左右边界面进行水平约束，下底边界面进行水平及竖向约束；采用非线性静力分析，计算得到三维边坡的应力场和应变场，经计算求解后，其中位移等高线云图如图10所示；

分别对其它所得应力应变场进行提取并存储，以下不再一一列出；

(iii)计算结果对比分析；

基于本发明方法，通过编制计算，共划分270根条柱，2710个搜索网格，经过具体实施步骤(1)、(2)、(3)(即建立搜索网格——搜索目标数值计算单元——搜索潜在破坏单元)，得到任意条柱内潜在破坏单元，提取其中第95、140根条柱潜在破坏单元搜索过程图解显示，如图11、图12所示；

如上所示，本发明方法可根据工程需要，找出任意边坡内部区域潜在破坏位置，并能找出最危险的滑面及其对应安全系数；

经过实施步骤(4)(即生成最危险滑移面，计算三维边坡的整体安全系数和局部区域安全系数)，最终生成基坑三维空间滑移面及安全系数空间分布，如图13、图14所示；

以上是通过本发明方法得到的基坑整体三维空间滑移面及基坑局部三维空间滑移面，滑移面形状为任意，土体应力状态任意，无相关假设条件，在求解滑移面位置和形状的同时也得到了基坑安全系数大小及空间分布，这种形式滑移面非常有利于工程设计与技术人员可直观、定量的判断出基坑最危险的位置，这是现有一些边坡稳定性评价方法不能实现的功能；

在得到基坑三维滑移面的同时，也能在指定断面位置生成平面条件下的二维滑移面，如图6中所示的1-1、2-2、3-3断面，其中，1-1位于基坑长边中部区域，2-2位于基坑长边端部基坑角部区域，3-3位于基坑短边中部区域，1-1、2-2、3-3断面除了所在位置不同之外，其它设计参数均相同；

利用本方法分别生成三个断面位置的二维滑移面及安全系数分布，如图15、图16、图17所示；其中虚线曲线表示用传统圆弧条分法计算得到的滑移面，实线曲线表示用本发明方法计算得到的滑移面；通过对比分析发现：(a)传统圆弧条分法在不同位置的断面处生成的滑移面位置均相同，安全系数大小也相同，且只能得到整体的安全系数K_W＝0.828；(b)本发明方法在不同位置的断面处的滑移面位置不同，整体安全系数值分别为：K_W1-1＝0.872，K_W2-2＝1.063，K_W3-3＝0.956，而且也同时得到安全系数在滑移面上分布；(c)1-1断面位置，传统的圆弧条分法与本发明方法安全系数相对接近，但在2-2及3-3断面处安全系数相差较大，这是由于2-2断面位于角部，3-3位于基坑短边区域，这两个位置相对于1-1断面三空间效应明显，尤其是2-2基坑角部空间效应最为显著，这点工程实际状况一致；

综合来看，以上3点体现了基坑三维边坡的空间效应的影响，以传统的圆弧条分法为代表的极限平衡方法不能考虑基坑空间效应的影响，而本发明方法基于三维应力状态进行边坡稳定性分析充分体现了空间效应的影响。对于三维边坡，采用本发明方法不仅无须对滑移面形状、土体应力状态做出假设，而且能考虑三维空间效应对边坡稳定性的影响，最终计算结果更符合三维边坡实际受力状况，且有利于边坡工程领域的设计和技术人员应用。

Claims (8)

1.一种基于应力矢量的三维边坡稳定性评价方法，其特征在于，步骤如下：

(1)建立搜索网格；

建立x轴平行于三维边坡宽、y轴平行于三维边坡长且z轴平行于三维边坡高的三维直角坐标系x-y-z，将三维边坡宽、长和高所围成的最小立方体离散成n×m×k个立方体搜索网格，对每个搜索网格进行编号，即完成搜索网格的建立；

n×m×k个立方体搜索网格组成了n×m个条柱，各条柱的底面位于三维边坡宽和长所在平面上，同时各条柱的高等于三维边坡高，m为所有的条柱排列形成的平行于x轴的行数，n为所有的条柱排列形成的平行于y轴的列数，k为各条柱中立方体搜索网格的数量，n、m、k的取值满足：各搜索网格内包含一个以上数值计算单元；

(2)搜索目标数值计算单元；

设搜索网格i的几何中心为Z点，Z点的空间坐标为(x_iz,y_iz,z_iz)；搜索网格i的八个顶点为V_i1、V_i2、V_i3、V_i4、V_i5、V_i6、V_i7、V_i8，八个顶点的空间坐标分别对应为(x_i1,y_i1,z_i1)、(x_i2,y_i2,z_i2)、……、(x_i8,y_i8,z_i8)；搜索网格i的数值计算单元j的四个顶点为A、B、C、D，四个顶点的空间坐标分别对应为(x_ja,y_ja,z_ja)、(x_jb,y_jb,z_jb)、(x_jc,y_jc,z_jc)、(x_jd,y_jd,z_jd)；搜索同时满足以下四个条件的数值计算单元，即为搜索网格i的目标数值计算单元；

①x_{i_min}<x_ja、x_jb、x_jc、x_jd<x_{i_max}；

②y_{i_min}<y_ja、y_jb、y_jc、y_jd<y_{i_max}；

③z_{i_min}<z_ja、z_jb、z_jc、z_jd<z_{i_max}；

④v_{j_abcd}＝v_{j_bcdz}+v_{j_cadz}+v_{j_dabz}+v_{j_bacz}；

其中，x_{i_min}、x_{i_max}、y_{i_min}、y_{i_max}、z_{i_min}、z_{i_max}分别为搜索网格i的八个顶点的空间坐标中的最小x坐标、最大x坐标、最小y坐标、最大y坐标、最小z坐标、最大z坐标；v_{j_abcd}表示搜索网格i的数值计算单元j的体积；v_{j_bcdz}表示由Z点和数值计算单元j的B、C、D顶点围成的四面体体积；v_{j_cadz}表示由Z点和数值计算单元j的C、A、D顶点围成的四面体体积；v_{j_dabz}表示由Z点和数值计算单元j的D、A、B顶点围成的四面体体积；v_{j_bacz}表示由Z点和数值计算单元j的B、A、C顶点围成的四面体体积；

按此方法，依次确定n×m×k个立方体搜索网格对应的目标数值计算单元；

若搜索网格内不包含目标数值计算单元，则该搜索网格为无效网格，在计算过程中按零标记，不参与计算过程；

(3)搜索潜在破坏单元；

搜索每个条柱的k个搜索网格对应的目标数值计算单元中同时满足以下两个条件的目标数值计算单元，该目标数值计算单元对应的搜索网格即为该条柱的潜在破坏单元；

①f≤[f],

②ε_p≥[ε_p]；

其中，f表示目标数值计算单元的安全系数；[f]表示安全系数允许值；C为土体的粘聚力；σ_n表示目标数值计算单元所在位置处任意滑移面上法线方向正应力；

为土体的摩擦角；τ_d表示目标数值计算单元所在位置处任意滑移面上下滑方向剪应力；ε_p表示目标数值计算单元对应的塑性应变；[ε_p]表示滑移面搜索时允许的塑性应变值，取值大小根据边坡安全控制要求确定，取值范围介于0～1％ε_pmax之间，ε_pmax表示边坡最大塑性应变值；

σ_n和τ_d通过以下公式计算得到：

其中，F_N′表示F_N的转置；F_N表示任意方向平面上正应力矢量；N表示该任意方向平面的法线方向；

表示该任意方向平面的法向量；F_D′表示F_D的转置；F_D表示任意方向平面上剪应力矢量；D表示该任意方向平面的下滑方向；

表示该任意方向平面的下滑方向矢量；σ_pq表示目标数值计算单元对应的应力张量，p＝1,2,3，q＝1,2,3，σ₁₁、σ₂₂、σ₃₃、τ₂₁、τ₁₂、τ₃₁、τ₁₃、τ₃₂、τ₂₃为应力张量的9个分量；

和

的确定流程如下：

(I)生成滑移面；

将搜索网格i的12条边进行分段，以所有边上分段后形成的断点为集合，取其中不共线的任意三点S_i1、S_i2、S_i3为组合，由S_i1、S_i2、S_i3点形成的面即为搜索网格i的滑移面；

(II)判断S_i1、S_i2、S_i3点是否同时位于搜索网格i所在条柱z轴方向的4条棱上，如果是，则将S_i1、S_i2、S_i3点作为H、I、J点进入步骤(IV)；反之，则进入步骤(III)；

(III)确定S_i1、S_i2、S_i3点所在的平面与搜索网格i所在条柱z轴方向的4条棱的四个交点H、I、J、K；

(IV)计算

和

的计算公式如下：

其中，

表示面HIJ的法向量；

表示点H与点I组成的线段的向量；

表示点H与点J组成的线段的向量；x_H、y_H、z_H分别对应为H点的x坐标、y坐标、z坐标；z_I、y_I、z_I分别对应为I点的x坐标、y坐标、z坐标；x_J、y_J、z_J分别对应为J点的x坐标、y坐标、z坐标；

的计算公式如下：

①当H、I、J点的z坐标相等时，

②当只有H、J点的z坐标相等时，

③当只有H、I点的z坐标相等时，

④当H、I、J点的z坐标均不相等时，

如果存在多个目标数值计算单元满足以上两个条件，则将f值最小的目标数值计算单元对应的搜索网格作为潜在破坏单元；如果不存在目标数值计算单元满足以上两个条件，则以包含目标数值计算单元且最靠近该条柱顶部的搜索网格作为潜在破坏单元；

按此方法，依次确定n×m个条柱的潜在破坏单元；

(4)将n×m个条柱的潜在破坏单元的中心点空间坐标拟合生成空间曲面，形成三维边坡的整体空间最危险滑移面，同时按以下公式计算三维边坡的整体安全系数和局部区域安全系数，并将其与参考值进行比较，定量评价边坡稳定性，参考值为与边坡相关的国家及地方规范规定的允许安全系数值；

式中，K_W表示三维边坡的整体安全系数；d表示潜在破坏单元的总数；k_v表示第v个潜在破坏单元的安全系数；σ_nv表示第v个潜在破坏单元的滑移面上法线方向正应力；τ_dv表示第v个潜在破坏单元的滑移面上下滑方向剪应力；C为土体的粘聚力；

为土体的摩擦角；

式中，K_L表示三维边坡的局部区域安全系数；e表示需要求解的局部区域所包含的潜在破坏单元的总数；k_g表示需要求解的局部区域所包含的第g个潜在破坏单元的安全系数；σ_ng表示需要求解的局部区域所包含的第g个潜在破坏单元的滑移面上法线方向正应力；τ_dg表示需要求解的局部区域所包含的第g个潜在破坏单元的滑移面上下滑方向剪应力；C为土体的粘聚力；

为土体的摩擦角。

2.根据权利要求1所述的一种基于应力矢量的三维边坡稳定性评价方法，其特征在于，步骤(1)具体如下：

(1.1)根据三维边坡的几何尺寸建立三维直角坐标系x-y-z，三维边坡宽平行于x轴，三维边坡长平行于y轴，三维边坡高平行于z轴；

(1.2)将三维边坡宽、长和高所围成的最小立方体划分为沿x轴依次排列的n个大立方体；

(1.3)将每个大立方体划分为沿y轴依次排列的m个中立方体；

(1.4)将每个中立方体划分为沿z轴依次排列的k个小立方体，至此，三维边坡宽、长和高所围成的最小立方体离散成n×m×k个立方体搜索网格；

(1.5)对每个搜索网格进行编号，即完成搜索网格的建立。

3.根据权利要求2所述的一种基于应力矢量的三维边坡稳定性评价方法，其特征在于，所有的划分都为均匀划分。

4.根据权利要求1所述的一种基于应力矢量的三维边坡稳定性评价方法，其特征在于，步骤(1)中，各搜索网格内数值计算单元的个数为1～125。

5.根据权利要求1所述的一种基于应力矢量的三维边坡稳定性评价方法，其特征在于，步骤(3)中，搜索每个条柱的k个搜索网格对应的目标数值计算单元中同时满足以下两个条件的目标数值计算单元时，按从下往上或从上往下的顺序。

6.根据权利要求1所述的一种基于应力矢量的三维边坡稳定性评价方法，其特征在于，步骤(3)中，将搜索网格i的12条边进行分段是指均匀分段。

7.根据权利要求1或6所述的一种基于应力矢量的三维边坡稳定性评价方法，其特征在于，步骤(3)中，搜索网格i的每条边被分为3～10段。

8.根据权利要求1所述的一种基于应力矢量的三维边坡稳定性评价方法，其特征在于，步骤(4)中，与边坡相关的国家及地方规范包括《建筑边坡工程技术规范》GB50330-2013和《建筑基坑支护技术规程》JGJ120-2012。

Images